Environmental fluctuations and the maintenance of genetic diversity in age or stage-structured populations

Bulletin of Mathematical Biology, Vol. 0

1996 Society

58, No. 1, pp. 103-127,

1996

Elsetier

Inc.

Science

for Mathematical 0092-8240/96

$15.00

Biology + 0.00

0092-8240(95)00309-E

ENVIRONMENTAL FLUCTUATIONS AND THE MAINTENANCE OF GENETIC DIVERSITY IN AGE OR STAGE-STRUCTURED POPULATIONS n STEPHEN

ELLNER

Biomathematics Graduate Program, Department of Statistics, North Carolina State University, Raleigh, NC 2769543203, U.S.A. (Email:

[email protected])

The ability of random fluctuations in selection to maintain genetic diversity is greatly increased when generations overlap. This result has been derived previously using genetic models with very special assumptions about the population age structure. Here we explore its robustness in more realistic population models, with very general age structure or physiological structure. For a range of genetic models (haploid, diploid, single and multilocus) we find that the condition for maintaining genetic diversity generalizes almost without change. Genetic diversity is maintained by selection if a product of the form (generation overlap) x (selection intensity) X (variability in the selection regime) is sufficiently large, where the generation overlap is measured in units of Fisher’s reproductive value. This conclusion is based on a local evolutionary stability analysis, which differs from the standard “protected polymorphism” criterion for the maintenance of genetic diversity. Simulation results match the predictions from the local stability analysis, but not those from the protected polymorphism criterion. The condition obtained here for maintaining genetic diversity requires fitness fluctuations that are substantial but well within the range observed in many studies of natural populations.

1. Introduction. The potential impact on evolution of random fluctuations in the selection regime has received considerable theoretical attention, from the seminal analysis by Haldane and Jayakar (1963) up to recent times (Gillespie, 1991). One general conclusion from numerous studies by mathematical population geneticists (see e.g. Hedrick, 1986; Karlin, 1988; Turelli, 1988; Barton and Turelli, 1989; Frank and Slatkin, 1990) has been that temporal fluctuations in selection can contribute significantly to the maintenance of genetic diversity in a population only under very restrictive conditions, unless the fluctuations act to induce heterozygote advantage (as in Gillespie’s SAS-CFF models (Gillespie, 1991). However, this conclusion is based on the analysis of models with the conventional assumption of discrete, non-overlapping generations. Recently it was shown (Ellner and 103

104

S. ELLNER

Hairston, 1994; Ellner and Sasaki, 1995) that with overlapping generations, fluctuating selection acting on a quantitative trait can maintain genetic diversity as long as the fluctuations are sufficiently large. Specifically, Ellner and Hairston (1994; hereafter, EH) considered stabilizing selection on a scalar quantitative trait, with random year-to-year fluctuations in the optimum trait value, and selection affecting the survival of newborn offspring. For variance of the optimum below a threshold value, the evolutionarily stable population is monomorphic, while for variance above the threshold the population is genetically polymorphic. The value of the threshold depends on the form of selection. An important and frequently used example is the Gaussian selection function, in which fitness is proportional to exp( - i, where z is the offspring’s trait value and 8, is the optimal trait value in year t. The “selection variance” a: is an inverse measure of the strength of selection-selection is intense if uw2 is small. Under Gaussian selection, genetic variance is maintained if yue2/u: > 1 (8 = generation overlap, u02= Var(8,)). This result applies to haploid as well as diploid inheritance, which is significant because it demonstrates that the mechanism does not depend on any form of heterozygote advantage. EH’s results are analogous to results for ecological models of interspecific competition, which showed that environmental fluctuations can promote coexistence via the “storage effect” when generations overlap (Chesson and Warner, 1981) and more generally (see Chesson, 1994, and references therein). However, whereas the ecological models consider competition between a finite number of unchanging species, the EH model addresses natural selection acting on a population containing a continuum of possible values for the trait of interest (e.g. continuous variation in height, weight, timing of reproduction, etc.). This difference in assumptions is reflected in different stability criteria, and the necessary condition for stable coexistence of alleles in the genetic models is more stringent that those for stable coexistence of species in the ecological models (see EH for discussion). The goal of this paper is to show that EH’s results continue to hold under more realistic models of overlapping generations. The EH model distinguishes only between new offspring and adults. New offspring are subject to selection, while adults experience a constant mortality rate that is dependent on age, time and genotype; newborns become adults as soon as they reach age 1. Clearly the EH model is a very special case of overlapping generations, and it would be a poor description of many if not most age-structured populations (such as humans). Before the theory can (or should) be tested empirically, it must be determined to what extent the results depend on the special form of generation overlap.

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

105

Here I consider a population in which survival and fecundity may depend on an individual’s age or stage, where “stage” is a classification based on an attribute, such as size, that changes over an individual’s lifetime. The population structure and transition rules are arbitrary, apart from an irreducibility assumption that amounts to saying that there are no disconnected sub-populations. The unexpected conclusion is that this increase in generality does not complicate the results. Indeed for Gaussian stabilizing selection the condition for maintaining genetic variance is not changed at all, except that y, the amount of generation overlap, requires a more general definition that covers both structured and unstructured models. The mathematical analysis is a local analysis in which the perturbations in allele frequency and the deviations between mutant and wild-type alleles are both assumed to be small. Small-fluctuation expansions, originally derived in mathematical demography (Tuljapurkar, 1989, 19901, are used to find the rate of increase of a rare allele’s frequency. The “small fluctuations” are in the fitness of the rare allele relative to population mean fitness. Local and global stability conditions need not coincide exactly, even in the unstructured model, but they are nearly identical for biologically reasonable parameter values, and local stability analysis gives very accurate predictions of the model’s long-term dynamics (Ellner and Sasaki, 1995). To demonstrate that the local stability analysis remains reliable for the structured model, simulation results are presented here for a model with three stage classes, that is a rough model for the freshwater copepod populations that originally motivated EH and this study. The assumptions and model are described in section 2, and section 3 discusses the stability criterion. Sections 4 and 5 contain the mathematical analysis, and the simulation results are discussed in section 6.

2. Population Model and Genotype Frequency Dynamics. I consider a discrete-time model for a structured population with overlapping generations. Structured means that the population is divided into a finite number of stages according to some attribute, which could be a life-cycle stage or any other attribute that varies over an individual’s lifetime, such as age, size, habitat, location, etc. The model also allows for the possibility that there are several types of newborns that would be assigned to different stages at birth (e.g. seedlings of different sizes). Population size (number of individuals) is assumed to be constant, and moreover I assume that the number of individuals in each stage is constant. These assumptions hold asymptotically in many structured models with density-dependent recruitment, and transients can be ignored for the analysis here. Thus in each time step of the model, a constant fraction of

106

S. ELLNER

the population dies and is replaced by newborn offspring. The population turnover is determined by the number of individuals in each stage and the stage-specific mortality rate. These assumptions are analogous to the “soft selection” model for selection in a spatially subdivided population, in which each sub-population contributes a fixed fraction of the individuals in the next generation. I assume that selection on the trait of interest acts only on newborns. Selection is assumed to be stabilizing with a fluctuating optimum. Specifically, the relative fitness in year t of a newborn with phenotype Z (a scalar quantity) is given by k(Z - 19,>,where 0, is the (random) optimum value of the phenotype in year t, and the fitness function k is a smooth symmetric function with k(O) = 1 and k(lZl) monotonically decreasing. In cases where newborns are distributed among several stage classes, the same relative fitnesses are assumed to hold within each class. Because the population sizes are constant, the state variables for the model are ~$2, t), the fraction of individuals in class i who have genotype z in year t. For simplicity, consider first the case of asexual haploid single-locus inheritance. For this case we can assume there is a one-to-one correspondence between phenotype (measured trait value) and genotype, and that the genotype index z is equal to the phenotype. I also assume for exposition that the population contains at most a countable set of phenotypes, but a continuous distribution of types is easily handled (replace sums by integrals in the obvious way). Let: ni = fraction of the total population in class i (i = 0,1,2,. . . , K). Fij = number of class-i offspring (prior to selection) produced by a class-j individual. w(z, t> = k(z-&,), the relative fitness of genotype z offspring in year t. Pij = annual transition rate from class j to class i (i.e. of the nj individuals “now” in class j, Pijnj will be alive and in class i “next year”). Sij = fraction of class-i individuals derived from survival from last year of class-j individuals, Sij = Pijnj/ni. Ri = fraction of class-i individuals derived from new recruitment, Ri = 11 cj sij. Recall that all 12,F, P and S values are assumed to be constant over time. I assume, for simplicity, that all the life cycle defined by the P and F values is irreducible and primitive, which is generally the case in demographic models if post-reproductive life stages are ignored (Caswell, 1989).

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

107

For any class i into which offspring do not recruit, gene frequencies “now” are simply a weighted average of the gene frequencies “last year” in the classes that feed into i: x,(z,t)

= C Siixj(t)

for i = 1,2 )...) K.

For a class i into which there is recruitment, are type 2 in year t is

(1)

the fraction of recruits that

(2)

Let +ij denote the fractional contribution recruits entering class i prior to selection,

4ij =nj&j

I

of class j to the pool of new

C IZkFik>

(3)

k

and let Fiji(t) denote the average fitness of offspring in class i in year t,

‘itt)

=

C C Y

4ijxjtY,

t)ktY

-

0,).

(4)

i

Substituting (3) and (4) into (21, the fraction of recruits that are type z is

x,(t + 1, z) can then be computed as a weighted average of the frequencies of type z in new recruits and in survivors. Let Ri be the fraction of class i derived from recruitment, Ri = 1 - Cj Sij, and let Bij be the fraction of class i derived from recruitment from parents in class j, Bij = Ri4ij. We thus have, for any class i, Xi(Z,

t +

1) = C SijXj( Z, t) + C Bijxj( zz, t)k( 2 - e,)/Wi( t) i i

(5)

108

S. ELLNER

Equation (5) can be expressed as the matrix equation

(6)

where q(z, t) = k(Z - e,>/iq(t>.

(7)

Equation (6) is not actually a linear system because the xi(z, t) appear in the formula for tl;i (equation (4)) and therefore in ri(z, t). Iterations similar to (6) can be derived similarly for more general models, with the details depending on the genetic system (haploid vs diploid, number of loci) and sexual system (monoecious vs dioecious, random vs assortative mating, etc.). The general principle is that genotype frequencies in newborns are generated by mating and reproduction, modified by selection and finally combined with the frequency of survivors from previous rounds of recruitment. A one-locus diploid model is considered below. 3. Evolutionary Genetic Stability. The goal of the analysis here is to identify the conditions under which a genetically monomorphic population is evolutionarily stable vs those in which any monomorphic population is evolutionarily unstable and therefore genetic diversity is maintained (Cohen and Levin 1991). A formal definition of evolutionary stability for genetic models was given by Eshel and Feldman (1982) under the name evolutionary genetic stab&@ (EGS). The EGS stability criterion requires that any rare mutant allele is at a selective disadvantage relative to “wild-type” alleles

FLUCTUATING

SELECTION:

STRUCTURED

POPULATIONS

109

established in the population. The biological assumption underlying the use of the EGS criterion is that the population is subject to recurrent mutation at loci affecting the trait of interest, and as a result the only alleles maintained in the long run are those which have a selective advantage over possible mutants. Thus, before applying the EGS criterion some biological constraints must be placed on the set of possible mutants (e.g. to disallow Drosophila mutants with the flying speed of a hawk). Here we assume that mutations only affect the value of the single trait Z under selection. The EGS criterion is conceptually very close to the evolutionarily stable strategy (ESS) solution concept from evolutionary game theory, but it takes account of the genetic dynamics rather than simply assuming that high-fitness phenotypes will increase in frequency. Consequently, genetic parameters which are ignored in a conventional ESS analysis (mode of inheritance, recombination rate, etc.) may strongly affect the results. For example, in the unstructured version of the present model, Ellner and Sasaki (1995) found that the number of loci controlling the trait affected the pattern of bifurcations in the genotype distribution as the variance in the selection regime is increased. However, the EGS stability criterion ignores genetic drift and assumes that evolution of the trait is not mutation-limited. To formally define evolutionary stability, we use the invasibility criterion for the increase vs decrease of a rare invading allele (see e.g. Turelli, 1978; Chesson and Ellner, 1989; Metz et al., 1992). We assume that for any initial set of alleles, the model has a unique stationary distribution p* for the joint distribution of genotype frequencies and optimum phenotype 8, and the stationary dynamics are ergodic. For any such stationary distribution, let &( p*> denote the set of alleles present in the population (ti is well defined, since ergodicity implies that an allele must either be present at all times or absent at all times). An ESS is a set Y of alleles maintained in a polymorphism (i.e. P’=&‘( ~“1 for some stationary distribution), such that any allele z PP has a negative boundary growth rate p(z). The boundary growth rate for an allele z is defined to be its long-run rate of increase or decrease when introduced as a rare invader:

p(z) = ,$n,

f

E ln{~,+,(z)/~,(z)},

(8)

where P,(z) is the frequency of z in year t and the expectation is taken over p”. The invasibility criterion generalizes to a continuum-of-alleles setting simply by introducing the type-z mutants as a delta-function perturbation to the population state, and then (8) is calculated for the growth or decay

110

S. ELLNER

over time in the magnitude of that component in the population distribution. If the population already contains type-z alleles, (8) is computed for a new sub-population of rare mutants with the same phenotypic effect and P,(z) is the frequency of z-alleles descended from the mutants. Invasibility analysis has the virtue of working so far on the type of model considered here, i.e. it correctly predicts the results of simulations for both single-locus and multi-locus models (Elmer and Hairston, 1994; Sasaki and Ellner, 1995; Ellner and Sasaki, 1995). This is no mean feat, since multi-locus quantitative trait models have been extremely difficult to analyze even under constant selection with discrete generations (e.g. Barton and Turelli, 1991; Nagylaki, 1993). However, the invasibility approach is still a heuristic linear stability analysis which lacks a completely rigorous mathematical justification. It has been proved valid for a general class of unstructured models (Chesson and Ellner, 1989; Ellner, 1989), but not for the structured models considered here, so conclusions based on an invasibility analysis must be checked by simulations. 4. Stability of Monomorphic Population: Haploid One-locus Model. In this section I consider a monomorphic population of type-y individuals, with type z introduced as a rare invader. For /z -yl small, I derive expansions of the boundary growth rate for z invading y and use the result to derive the conditions for evolutionary stability of genetic monomorphism. The main conclusion is that for the Gaussian selection model, the ESS is genetically polymorphic whenever ys2/crz > 1, where y is the generation overlap when individuals are weighted by their reproductive value as calculated from the transition matrix for a monomorphic population. For non-Gaussian selection, it is still the case that higher generation overlap and stronger selection destabilize a monomorphic population and promote the maintenance of genetic variability. Assume to begin with that selection is Gaussian. Without loss of generality, the measurement unit for the trait can be chosen so that ai = 1, E( 0,) = 0, and henceforth this is always assumed to hold for Gaussian selection model. To calculate the boundary growth rate for z invading y, (6) is linearized about zero frequency for type z. The result is easily seen to be (6) with z deleted from all sums defining the Ei, leaving only type y in the sum. The linearized dynamics for the frequency of type z are then given by (6) with Ti(z, t, = k(z - 8,)/k(y

- e,)

(9)

for all i. The matrices in (6) are therefore independent and identically distributed, so the boundary growth rate is given by the dominant Lyapunov

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

111

exponent of the random matrix product process (6). Moreover, Taylor expansion shows that ~$2, t) = 1 + O(z -y) for Iz -yI small, and the rest of the matrix is constant over time and independent of z and y, so the matrices in (6) have “small fluctuations.” Thus the small-fluctuations expansion for the dominant Lyapunov exponent (e.g. Tuljapurkar, 1989, 1990) can be used to approximate the boundary growth rate, as follows. Let A(t) denote the matrix in equation (6) for z invading y and let p(z, t> = Y(Z,t) - 1. We can write A(t) =A, + ( ,@B + ( P(V)

(10)

- ( ,@)B,

where ( p) = EP(z, t), A, is the matrix that results from setting Y(Z,t) = 1 in A(t) and B is the matrix whose ijth element is Bij. The first step is to approximate the dominant Lyapunov exponent of the mean matrix A, + ( p >B. Since A, + ( /3)B is constant, its dominant exponent is simply the logarithm of its dominant eigenvalue. A, is non-negative and has A = 1 as an eigenvalue with corresponding right eigenvector e = (1, 1,. . . , 1); therefore, 1 is the dominant eigenvalue of A, (Horn and Johnson, 1985, Corollary 8.1.30). Let u = (u,, uI, u2,. . . , u,) be a dominant left eigenvector of A,,, normalized so that its entries sum to 1. (All entries in the dominant left eigenvector are positive because the life cycle is irreducible, so this normalization is possible). The well-known formula for eigenvalue sensitivities in terms of left and right eigenvectors (see e.g. Caswell, 1989, equation 6.8) gives dh/dAij = ui for A, and therefore

where dh Bij = C ~,Bij = C ~iRi~ij = C ~iRi. i j dAij i i,i i,i

@=C

(11)

The second step in approximating p(A(t)) is to use the small-fluctuations expansion to approximate the difference between p( A, + ( p >B) and p(A(t)). Suppose first that the wild-type phenotype y # 0. The leading term in the expansion is proportional to the variances and covariances of the perturbations to matrix entries. The random perturbations are ( p(z, t) ( P))B which is O(z -y), so we have p(A) = p(A, + ( /3)B) + O((z -yj2). On the other hand, exp( - (0 - ~)*/2) P(zyt’

=

exp( -(B-y)*/2)

- 1=

(2 -ym

-y)

+ O((z

-y,“)

112

S. ELLNER

and hence ( p > = -y(z -y) + O((z -yj2>, which implies that p(Ao + ( /3 )B) differs from p( A,) by an amount of O(z - y). Thus to leading order in (z -y) we have p(A(t)) = P(A, + ( /?)B). B is non-negative and has at least one positive entry (i.e. some class has non-zero fecundity); hence p(A) > 1 for ( p ) > 0 sufficiently small. This occurs for any z having the same sign as y but slightly smaller magnitude. Therefore, a population monomorphic at y # 0 can always be invaded by a nearby type and cannot ever be an ESS. Now consider a population monomorphic at y = 0, invaded by a nearby

type t. In this case, Taylor expansion of (9) gives ( p) = z”( a,2 - 1)/2 + 0(23> Var( p(z, t>> = z2cq2

(12)

and, therefore, by (ll), p(A, + ( P)B)

= @z2(q2

- 1)/2 + 0(z3).

The leading order perturbation term in the small-fluctuations can be expressed as (Tuljapurkar, 1989, equation 4.2.7)

(13) expansion

(14) where C(q, kl) is the covariance between the (i, j) and (k, 0 entries in A(t), A is the dominant eigenvalue of the mean matrix A,, + ( /3 >B and Aij is the sensitivity (partial derivative) of A with respect to the changes in the (i, j)th entry of the mean matrix. In the present case, using (12),

C(ij,kl) =BijBklVar(

P(z, t)) =BijBklz2q2_

Therefore, to compute (14) to 0(z2>, A and the A, only need to be computed to O(1). Thus we can use 1 as the value of A in (141, and for the Aij we can use the corresponding values for A,, namely Aij = ui. To 0(z2), therefore, (14) equals

-

7

& r,l,k,l

Bij~iBkluk

=

-@;2u’ .

(15)

FLUCTUATING

SELECTION:

STRUCTURED

POPULATIONS

113

Combining (13) and (15) we have that &4(t))

=A((1

- @‘)a,2 - 1) + Ok31

(16)

for type z invading a monomorphic population at phenotype 0. When (16) is positive, any monomorphic population can be invaded, and therefore genetic variance must be maintained by selection. The condition for maintaining genetic variance is therefore (1 - @‘)$ > 1,

(174

in scaled units such that aw2= 1, and (1 - @)ae”/a,2 > 1,

(17b)

in unscaled units, for the Gaussian selection function. It is illuminating to express (1 - @,) in terms of the stage-specific survivorships and reproductive values. Let P and F be the survival and fecundity matrices for a monomorphic population (that is, the matrices whose (i, j)th entries are Pij and Fij, respectively). The state-specific reproductive values for a monomorphic population are given by the dominant left eigenvalue of (P + F). Recall that Pij = (ni/ni)Sij and note that by definition Ri = C, nkFik/ni, SO that by equation (31, Bij = Ri~ij = (r~~/rz~)&~. Therefore, A,=D-‘(P+F)D

and

S=D-‘PD,

where D is the diagonal matrix with diagonal entries (n,, IZ,,. . ., n,). Because this represents a change of basis, the dominant left eigenvectors of A, and (P + F) (u and U, respectively) are related by u = uD. The generation overlap, when individuals are weighted according to their reproductive value, is given by def

y = uPn/un.

The denominator is the total reproductive value of the population, and the numerator is the total reproductive value of individuals that survive to the next year. We then have U(DSD-‘)n Y= =

u(DD-‘)n

=

(zD)(D-‘n)

use = C uisij = C i,i

=1-Q.

use

(uD)S(D-‘n)

i

ui

= ue =

C Ui(l i

-Ri)

114

S. ELLNER

The condition for maintaining becomes

genetic variance, equation (17b), therefore

yq2/u,2

> 1.

(17c)

This is, as claimed, identical to the condition derived for the unstructured model with Gaussian selection, apart from the more general definition of y in terms of reproductive value. For non-Gaussian selection, the analysis above is easy to repeat using Taylor expansion of (9) to find moments of p to the necessary order of accuracy. It is convenient to set k(x) = exp(--b * g(x)>, where g is a function with g(0) = 0 and g(l) = 1, and b therefore is a measure of the selection intensity. Taylor expanding (9) in z about t = y then gives ( p) = b(z - y)Eg’(y - 0,) + O((y - z>~), so a population monomorphic at phenotype y can be an ESS only if Eg '(y - 0,) = 0. If y” is any such “candidate ESS,” then the higher-order terms in the Taylor expansion can be used as above to obtain

p(AW) A ; Hz

-y*,"(ybV-B),

(18)

where Q= Var(g’(y*

- @,>),

g =E(g”(y*

exactly as in EH. The condition for maintaining fore that ybb+

- e,)),

(19)

genetic variance is there-

(20)

for all candidates ESSs y*. Note that the condition defining y* involves only the fitness function g(x), and not b or any of the transition or birth rates, so changes in those rates will only affect y in equation (20). Because y and Q are both positive, the effect of any increases in the generation overlap y or the strength of selection b will be to de-stabilize any candidate ESS y*. 5. Diploid Sexual Model. In this section, I show that the results of the previous section apply almost verbatim to structured populations with diploid sexual reproduction and additive allele effects. The intuitive reason for the carryover of results is that when a rare allele z invades a monomorphic population of type y; the population initially consists of rare (z, y) heterozygotes and common (y, y> homozygotes (i.e. all terms corresponding

FLUCTUATING

SELECTION: STRUCTURED

POPULATIONS

115

to (z, t) homozygote terms drop out in the linearization). Therefore the invading allele again faces a single type of competitor and increases according to the fitness of (z, y) relative to (y, y). Thus the boundary growth rate for allele z should be the same as that for z + y invading 2y in the haploid model of the previous section. The remainder of the section verifies that this heuristic argument holds up. Sexual reproduction adds only one significant complication to the model: offspring genotype frequencies depend on the mating pattern. A general assumption, which covers a variety of mating and sexual systems (see Charlesworth, 1980, Chapter 31, is the following: a fied fiactionp.ij,,of all class-m offspring (prior to selection) is derived j?om the union of randomly chosen stage-i male gametes, with randomly chosen stage-j female gametes. The

distributions pij,m summarize the net effect of assortative mating, stage-dependent mating success in males and stage-dependent reproductive success in females. Note that even if all individuals are hermaphroditic (e.g. many flowering plants which produce both pollen and ovules), we need the ~ij,(n to account for possible stage-dependent variation in reproductive rates via male vs female gametes from an individual. This model can also accommodate facultative clonal reproduction (e.g. fission of one large individual into two smaller individuals) as long as clonal “offspring” are not subject to selection on the trait of interest, simply by including clonal reproduction in the transition rates Pij. A second potential complication is that the sexes may have different survival and transition probabilities Sij. This would not occur if all individuals are hermaphroditic or if the “adults” are actually resting stages such as dormant eggs or seeds in which sexual differences are not expressed, but in general there is no reason why the sexes should have the same vital rates. It turns out that such gender differences do not cause any difficulties, but for the sake of exposition assume to begin with that males and females have identical vital rates. Then because selection is not gender biased, it is eventually the case that males and females (if there are such) have identical stage-dependent genotype frequencies, We can therefore use xj(z, t) to denote the frequency of allele z in stage j individuals, both male and female. Let (y,y’) denote an ordered genotype in which y and y’ are the paternally and maternally contributed alleles, respectively. The frequency of (y, y’) in unselected offspring is then (suppressing for now the offspring class subscript m and the time index t)

C i,j

/JijXitY)xjtY’),

(21)

116

S. ELLNER

which is then modified by selection to

C

+ty’)/W,

/4j(_Y)xj(Y')W(Y

i,i

(23)

where

C

K=

+_Y’>*

/JAijxi(Y)xj(Y’)w(Y

(24)

i,i,y,y’

The frequency of an allele z in selected offspring is then

C

tXi(Z) C

~ij{

+Y> +

Xj(Y>w(z

+Xj(Z)

Y+z

i,i

X

C Xi(Y>w(Z +Y>

)

/’

+Xi(Z)xj(z)w(2z)

Y#= =

=

=

C

Pij{

i,i

+xi(z)C

+Y) + +Xj(Z) C Xi(Y)W(Z

xj(Y)w(Z

Y

( Pij

Cl

+

+Y))/ic

Y

Pji)

Xi(Z> C

Li

2

I

C’i”)( i

C j,y

$ijxj(Y)W(Z

xj(Y)w(Z

+Y)/W

Y

(25)

+y)/E},

where *ij = ( CLij+ /.~ji)ji)J2. Thus the one-locus diploid model with equal vital rates for males and females takes the form of the haploid model (6) with R,{C,, ~ij,mXj(y, ~)w(z +Y, t)/if,(t)I in place of B,,~,(z, t). However, for allele z invading a population monomorphic for allele y, the linearized dynamics for z reduce to (6) with Btni =

RmC

$ij, m

i

q&G t> = w(z +y, t)/w(2y,

t).

(26)

Thus the boundary growth rate of z invading y in this diploid model is the boundary growth rate of z + y invading 2y in the haploid model of the previous section with the Bs defined as above. The calculations for the generation overlap y are also formally unchanged once the single-sex

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

117

definitions of relative contributions to different classes of newborns are appropriately modified for the two-sex case (Aj = Ck I,$.~,~,which gives Bij = c#+~R~ as in the one-sex case, and Fjj = (n,/nj)Bij>. Consequently the threshold conditions (17) and (20) for maintaining genetic variance still apply under diploid sexual reproduction. If the sexes have different vital rates, the situation is only slightly more complicated. The prior histories of the males and females currently in a given class will not have the same probability distribution, so we cannot assume equal allele frequencies in males and females. Thus the number of state variables must double: Let classes 1,2,. . . , K be females, while K + 1,K+2,..., 2K are males. Then xj(t, t) can still unambigously denote the frequency of allele 2 in class j. Let pij,m continue to denote the fraction of offspring in class m (prior to selection) that result from the union of a class-i male and class-j female gamete. With these definitions, many of the ps and many entries in the transition matrix will be 0 unless sex change is common, but it is easily seen that the calculations in this section remain valid without any modification whatsoever. The conclusions therefore continue to hold even when male and female vital rates are not equal. The only complication is that y will involve both male and female “reproductive values,” defined as above by the dominant left eigenvalue of the transition matrix for a genetically monomorphic population. 6. Simulations. The results derived above all carry the caveat that they are based on a local stability analysis. The analysis is “local” in two distinct and independent senses: (i) Stability is checked only against small changes in allele frequency (rare mutants). (ii) Invading mutants are assumed to have only small phenotypic differences from the wild-type allele(s) established in the population. Sense (i) justifies a linear stability analysis and sense (ii) justifies the small-fluctuations expansion for the boundary growth rate. For both senses, the analysis does not indicate the quantitative meaning of “small.” Simulation studies are necessary to determine whether the local stability analysis accurately captures the behaviour of the system. For the simulations I used a three-stage model based roughly on the population which motivated this work- the freshwater copepod Diuptomus sanguineus in Bullhead Pond, Rhode Island (see Fig. 1) The trait under selection is the timing of diapause-more precisely, the date in the Spring when a female switches from producing eggs that hatch within a few days, to making eggs that will remain in diapause in the pond sediment over the

118

S. ELLNER

New eggs

-

u

JI

v!’ h2

n2

(1-h2)(1-d)

AdUltS

(I-h2)d

n3

h3

s(l-h3) (1-h3)(1-8) v

Figure 1. Life cycle for the simulations, based on the vertically structured egg bank of the freshwater copepod Diaptomus sanguineus. Arrows indicate possible transitions and symbols labeling arrows are the transition rates. Stages 1, 2 and 3 represent eggs sitting above the sediment (on the surface of the pond sediment or on the leaves of submerged macrophytes), the “shallow” egg bank and the “deep” egg bank, respectively. Newly produced eggs spend their first hatching season above the sediment; those that do not hatch (a fraction 1 - h, of the total) become covered by a thin layer of sediments and therefore have a reduced hatching probability h,. Eggs in the deep egg bank are buried so deeply in the sediment that they can only hatch if the sediment is disturbed (e.g. by fish building a nest>. The “deep” egg bank is depleted by gradual mortality of older eggs. The number of new eggs each year is assumed to be constant, and hatched individuals from all stages make equal per-capita contributions to the pool of new eggs prior to selection.

summer, when predation eliminates all non-dormant individuals from the population. These diapausing eggs may remain viable in the pond sediments for many years, resulting in high generation overlap. For more information on this trait and the selection regime, see Hairston (1988), Hairston and Dillon (19901, and DeStasio (1989). The three-stage model was implemented with haploid asexual inheritance as in Section 4 and also with diploid sexual inheritance as in section 5. In the diploid model the vital rates for males and females were assumed to be equal, because Diuptomus eggs are not observed to express any gender differences. New offspring were assumed to be produced by random mating between the adults that develop from hatched eggs. The random mating assumption implies that the diploid genotypes of offspring prior to selection will be in the Hardy-Weinberg proportions determined by the allele frequencies in the hatched adult population. Consequently the diploid

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

119

model can take as its state variables the allele frequencies in the egg bank, rather than the genotype frequencies. The haploid and diploid models simulated allele frequency dynamics for 51 alleles evenly spaced between symmetric limits z,, and -z,,,. In the haploid model simulations I used z,, = 0.2 for a check of “local” stability (in the second sense above) and z,, = 2 for a check of “global“ stability. The fitness function was Gaussian with CT,+? = 1, so z,,, = 0.2 vs 2 are “small” vs “large” in terms of their effects on fitness. For the diploid simulations I used z,, = 0.1 and 1, so that the range of possible phenotypes was the same as in the haploid model simulations. Populations were initiated with allele frequencies proportional to exp( -(z/z,,>*) in all stages. Thus the simulations considered very large perturbations of allele frequencies from a monomorphic state. In order to explore the validity of the linear stability analysis, I used a range of parameter values, including values very different from those believed to hold for the real population of D. sanguineus. Specifically, the values of h,, h,, h,, d and s were chosen independently from uniform distributions on the intervals [0.1,0.8] for hi and d and [0.1.0.99] for s. Values of the optimum phenotype 0, in each generation were generated by uniform random draws from a discrete distribution at the points ti = xi

exp( CYX,~),wherexi=

-1+i/15,i=O,l,2

,...,

30.

For (Y= 0 this is a discrete uniform distribution, while for larger values of (Y the distribution becomes increasingly leptokurtic; at (Y= 1 the kurtosis is nearly that of a Gaussian distribution. The value of (Y was chosen uniformly on the interval [0,2] for each simulation. The value of y was calculated for each set of parameters using a formula (derived in the Appendix) that holds for the life cycle in Fig. 1. The simulation results for the haploid model (Fig. 2) show that the local stability analysis is correct even for large perturbations in allele frequencies as long as all possible mutations have small phenotypic effects (left column, 2 max= 0.2): for parameters below the threshold (171, the population is apparently converging to monomorphism, while above the threshold, genetic variance is maintained. As in the unstructured model (Ellner and Sasaki, 19951, the genetic variance consisted of two alleles with equal but opposite effects. If mutations can have large phenotypic effects (Fig. 2 right column, z max= 2.01, there is a small set of parameters for which the monomorphic equilibrium is locally stable but globally unstable. The occurrence of such parameters is shown in Fig. 2 by the high levels of genetic variance maintained in a few of the below-threshold simulations. In contrast there

120

S. ELLNER

Std. Dev. of Allele Distribution Figure 2. Simulation results for one-locus haploid model with the life cycle shown in Fig. 1. Each panel is a histogram showing the frequency distribution of the genotypic standard deviation (i.e. the standard deviation of the allele frequency distribution), based on simulations with 200 randomly chosen values for the parameters h,, h,, h,, d, s and (Y. For each set of parameter values, eight simulations were run of 100,000 generations of selection, and the standard deviation of the allele effects (2) was calculated for the average allele distribution of hatched individuals over the final 5000 generations. The eight simulations involved four different levels of the variance of fluctuations (E = ~-~-y/u: = 0.85,0.95, 1.05, 1.15) and two different sets of alleles: The left column results are for 51 evenly spaced alleles between - 0.2 and + 0.2; right column results are for 51 evenly spaced alleles between -2 and +2. According to the local stability analysis, the population should converge to monomorphism (standard deviation of allele distribution + 0) for E < 1 and to polymorphism (standard deviation of allele distribution > 0) for E > 1.

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

121

were no above-threshold simulations in which appreciable genetic variance was not maintained, indicating that the threshold (17) is a sufficient condition for maintaining genetic variance. The simulation results for the diploid model (Fig. 3) follow exactly the same patterns. One frequently used approach to the maintenance of genetic diversity by selection is to identify the conditions which allow a protected polymorphism

100

f

100

&=.95

:i

[

“8:

,g20‘O

,B ‘0

0

20 0 E

100

100

60

&=I.05

P :i

2

;&f;Gfs

‘O

H Y

20 0

E=l.os

60 40 20 0 f

;EtE;&$

100 ,

Std Dev of Allele Distribution

Std Dev of Allele Distribution

Figure 3. Simulation results as in Figure 2, for the one-locus diploid model with the life cycle shown in Fig. 1. In this case the left column results are for 51 evenly spaced alleles between -0.1 and +O.l; right column results are for 51 evenly spaced between - 1 and + 1. The possible phenotypes therefore ranged between + 0.2 and k 2 in the left and right columns, respectively. Local stability analysis predicts that the bifurcation to polymorphism occurs at E = 1.

122

S. ELLNER

to occur (e.g. Gillespie, 1991; Haldane and Jayakar, 1963; Hedrick, 1986; Karlin and Liberman, 1974). A protected polymorphism is a set of alleles which are maintained by selection in the absence of mutation: they must be “stable against each other,” but not necessarily against “stable against invading mutants.” Protected polymorphisms can occur in the present model whenever uo2> 0, which is true in all the simulations shown in Figs. 2 and 3. (Specifically, alleles z and -z comprise a protected polymorphism for z sufficiently small. Carrying out to order z2 the small fluctuations expansion of section 4 gives boundary growth rate 2y(l - y ) u,~z~ + 0(t3> for z invading -z and also for -z invading z by symmetry. Thus each allele will increase when rare, which is the defining property of a protected polymorphism.) Thus, the protected polymorphism criterion does not correctly predict the circumstances under which genetic variation is maintained when evolution is not mutation-limited.

7. Discussion. The main biological conclusion of this paper is that fluctuating selection can be an effective factor in maintaining genetic diversity in populations with overlapping generations for very general forms of generation overlap which include most any reasonable combinations of age, stage, sex and habitat structure. The threshold conditions for the maintenance of genetic diversity are no more complicated than in the unstructured model, and a single general condition applies across a broad range of genetic systems. For the Gaussian selection model, the criterion for maintaining genetic variance is given by equation (17). Because y I 1 by definition, the maintenance of genetic diversity requires ue2 > aw2,so that (for example) the mean optimum phenotype must often experience fluctuations in relative fitness of 50% between one year and the next. It is common for population genetics theory to assume that selection is much weaker than this. However, Endler (1986, Fig. 7.1) found that published estimates of selection coefficients on polymorphic traits in natural populations were distributed almost evenly over the range 0 (no selection) to 1 (100% effective selection against the inferior type), implying that stronger selection than (17) requires is common. Moreover, Hairston et al. (1995), looking specifically at published data on fluctuations in recruitment, found temporal fluctuations in the recruitment success of dormant propagules (seeds or eggs) by factors of 13.3 in desert/island annual plants, 1150 in chalk grassland annuals and biennials, 614 in chaparral perennials, 1150 in freshwater zooplankton and 31,600 in terrestrial insects. Among long-lived adults, recruitment has been found to vary among years by factors of 4 in desert perennials, 333 in temperate forest perennials, 591 in marine invertebrates, 706 in freshwater fish, 38 in

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

123

terrestrial vertebrates and 2200 in birds (Hairston et al., 1995). These data demonstrate that the required level of variance in selection can readily occur in traits related to the timing of reproduction. The fitness of a desert annual seed that correctly “chooses” to germinate in a year when there is enough rain for reproduction to succeed or correctly “chooses” to remain dormant in a year that is too dry would be orders of magnitude larger than that of a seed that acted incorrectly. It is also likely that (17) is strict relative to the conditions for maintaining genetic diversity in more general models. Three assumptions of the present model eliminate sources of variation in fitness that can potentially contribute to maintaining genetic diversity: (i) The population size of each stage or age class is constant. (ii) Adult vital rates are constant. (iii) There is no selection on traits affecting other stages of the life cycle that are genetically correlated with the trait of interest. Preliminary simulations indicate that violation of assumption (iii> can have a drastic quantitative effect. We (Ellner and Hairston, 1995) simulated a haploid two-locus unstructured model for a population with a diapausing egg bank. One locus controlled the trait subject to fluctuating selection, while the other controlled the hatching probability for diapausing eggs. The traits become genetically correlated as selection generates linkage disequilibrium, so they do not evolve independently. Genetic variation in both traits could be maintained in this model for a,* below the threshold in equation (17). Violation of assumption (i) had very little effect in simulations where population fluctuations are independent of the fluctuations in selection (Ellner and Sasaki, 1995). However, if extreme values of 0, result in high total recruitment, the threshold value of a,* for maintaining genetic diversity can be decreased. In D. sanguineus, years in which the optimum switch date is very late are years when fish predation comes late in the life cycle and is relatively weak (Hairston, 1988) and therefore egg production is higher, all else being equal. Endogenously generated population fluctuations due to overcompensatory density dependence in total recruitment (e.g. the Ricker model) can also decrease the threshold for maintaining genetic variance (D. Baba’i, in preparation). Age-structured population models with randomly fluctuating vital rates have been extensively studied over the past 20 years (see Tuljapurkar, 1989, 1990, for reviews), including some genetic models with genotype affecting the distributions of all age-specific vital rates. However, little has been said in those studies about the maintenance of genetic diversity, except to note that age structure does not preclude the maintenance of polymorphisms by

124

S. ELLNER

heterozygote advantage (Orzack and Tuljapurkar, 1989; Orzack, 1993). Templeton and Levin (1979) considered the effect of generation overlap in a seed pool on the response to fluctuating selection, but they incorrectly concluded (based on heurestic arguments) that a long-lived seed pool would reduce the impact of random fluctuations. The results here are closer in spirit to life history theory for fluctuating environments, which predicts that phenotypic heterogeneity can be an advantageous strategy for “bet hedging” against an unpredictable environment (Bull, 1987; Seger and Brockmann, 1987; Ludwig and Levin, 1991; Yoshimura and Clark, 1993). Taken together, this paper and its predecessors (Ellner and Hairston, 1994; Sasaki and Ellner, 1995; Ellner and Sasaki, 1995) demonstrate that an EGS local stability analysis can be a powerful approach for studying the long-term effects of natural selection in models for quantitative trait evolution. There are two main alternative approaches (Nagylaki 1993). The first approach uses explicit multi-locus genetic models, but “the great complexity and difficulty of the problem are reduced by imposing restrictions on the fitness pattern” (Nagylaki, 1993, p. 627) such as symmetries. The second, which has been heavily applied since its introduction by Lande in the 1970s (Lande, 1975, 1976, 1977, 1982), focuses on phenotypic trait distributions rather than genotype frequencies. The models are derived by approximating the generation-by-generation changes in trait distributions under various assumptions, such as weak selection of Gaussian trait distributions. This approach extends the classical methods of quantitative genetics for predicting the progress in trait means under artifical selection (e.g. Falconer, 1981). Recent papers have derived approximations based rigorously on explicit underlying genetic models, with exact results for some special cases (e.g. Barton and Turelli, 1991; Turelli and Barton, 1990; Nagylaki, 1993). Comparing this paper with the last three cited, it is apparent that the EGS approach is much simpler to use, even though the model here is in many ways more general (overlapping generations, strong selection, fluctuating selection). The drawbacks of the EGS approach are that it says nothing about the transient dynamics far from the equilibrium and it ignores effects of random mutation and drift. However, for many purposes all that matters is the long-term effect of selection, e.g. Barton and Turelli’s (1991) analysis of sexual selection. In that case, an EGS analysis may provide a simpler and more direct route to the desired results than an analysis based on global approximation of the transient dynamics. I thank Michael Turelli, Nelson G. Hairston, Jr. and an anonymous referee for useful comments, Akira Sasaki for graciously instructing me in population genetics while pretending to be my post-doctoral student and an

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

125

anonymous referee on a different paper who raised the issue of general age structure. This work was supported by NSF grant DEB911-8894 to SE and N. G. Hairston, Jr. Computing facilities were partially supported by an NSF SCREMS grant to the Department of Statistics at North Carolina State University. APPENDIX Calculation of y for the Three-stage Model of Fig. 1. The simplest way to calculate y is via the formula y = 1 - 0 = 1 - XuiRi = 1 - u,, which holds in this case because R, = 1, R, = R, = 0. To find u, we need to calculate B + S and fmd its dominant left eigenvector u. For the life cycle shown in Fig. 1, the first row of S is all OSbecause all eggs in stage 1 are newborns. We also have S,, = S,, = 0, and by definition S,, = Pz2 and S,, = Pj3. Therefore, S,, = 1 - S,, and S,, = 1 - S,,, giving 0

S= [ l-(l-h&l-d) 0

0

0

(l-h&l-d)

0

1 - s(l -h,)

SC1 - h3)

1 .

The nonzero entries in B are all in the first row, B,j = hjnj/C h,n,, j = 1,2,3. The relative values of the nj are obtained by solving the steady-state conditions n2 = (1 - h,)n, + (1 - h,)(l n3 = (1 - h,)dn,

+s(l

- d)n,

- h,)n,

With B and S in hand, a left eigenvector (l,q,,qJ of (B + S), corresponding to the dominant eigenvalue A = 1, is easily obtained from the first and third columns of (B + S) because of the zero entry in each. Then by definition u1 = l/(1 + q2 + q3), SO y = (q2 + qJ/(l + q2 + 43).

LITERATURE Barton, N. H. and M. Turelli. 1989. Evolutionary quantitative genetics: how little do we know? Ann. Reu. Genetics 23: 337-370. Barton, N. H. and M. Turelli 1991. Natural and sexual selection on many loci. Genetics 127, 229-255. Bull, J. J. 1987. Evolution of phenotypic variance. Evolution 41, 313-315. Caswell, H. 1989. Matrir Population Models, Sunderland MA: Sinauer Associates, Inc. Charlesworth, B. 1980. Evolution in Age-Structured Populations. Cambridge UK: Cambridge University Press. Chesson, P. L. 1994. Multispecies competition in varying environments. Theor. Population Biol. 45, 227-276. Chesson, P. L. and S. Ellner. 1989. Invasibility and stochastic boundedness in monotonic competition models. J. Math. Biol. 27, 117-138. Chesson, P. L. and R. R. Warner. 1981. Environmental variability promotes coexistence in lottery competitive systems. Am. Naturalist 117, 923-943.

126

S. ELLNER

Cohen, D. and S. A. Levin. 1991. Dispersal in patchy environments: the effects of temporal and spatial structure. Theor. Population Biol. 39, 63-99. De Stasio, B. T. 1989. The seed bank of a freshwater crustacean: copepodology for the plant ecologist. Ecology 70, 1377-1389. Ellner, S. 1989. Convergence to stationary distributions in two-species stochastic competition models. J. Math. Biol. 27, 451-462. Ellner, S. and N. G. Hairston, Jr. 1994. Role of overlapping generations in maintaining genetic variation in a fluctuating environment. Am. Naturalist 143, 403-417. Ellner, S. and A. Sasaki. 1995. Patterns of genetic polymorphism maintained by fluctuating selection with overlapping generations. Theor. Population Biol. (in press). Endler, J. A. 1986. Natural Selection in the Wild. Princeton University Press, Princeton, NJ. Eshel, I. and M. W. Feldman. 1982. On evolutionary genetic stability of the sex ratio. Theor. Population Biol. 21, 430-439.

Falconer, D. S. 1981. Introduction to Quantitative Genetics, 2nd ed. Harlow, UK: Longman Scientific&Technical. Frank, S. A. and M. Slatkin. 1990. Evolution in a variable environment. Am. Naturalist 136, 244-260.

Gillespie, J. H. 1991. The Causes of Molecular Evolution. Oxford, UK,: Oxford University Press. Hairston, N. G., Jr. 1988. Interannual variation in seasonal predation: its origin and ecological importance. Limnology and Oceanography 33, 1245-1253. Hairston, N. G. Jr. and T. A. Dillon. 1990. Fluctuating selection and response in a population of freshwater copepods. Evolution 44, 1796-1805. Hairston, N. G., Jr., S. Elmer and C. M. Kearns. 1995. Overlapping generations: the storage effect and the maintenance of biotic diversity. In Population Dynamics in Ecological Space and Time, 0. E. Rhodes, R. K. Chesser and M. H. Smith (Eds). Chicago: University of Chicago Press (in press). Haldane, J. B. S. and S. D. Jayakar. 1963. Polymorphism due to selection in varying directions. Genetics 58, 237-242. Hedrick, P. W. 1986. Genetic polymorphism in heterogeneous environments: a decade later. Ann. Rev. Ecology and Systematics 17, 535-566.

Horn, R. A. and C. R. Johnson. 1985. Matrix Analysis. Cambridge, U.K.: Cambridge University Press. Karlin, S. 1988. Non-Gaussian phenotypic models of quantitative traits. pp. 123-144 In The Second International Conference on Quantitative Genetics, E. J. Eisen, M. M. Goodman, G. Namkoong and B. S. Weir (Eds). Boston: Sinauer. Karlin, S. and U. Liberman. 1974. Random temporal variation in selection intensities: case of large population size. Theor. Population Biol. 6, 355-382. Lande, R. 1975. The maintenance of genetic variability by mutation in a polygenic character with linked loci. Genetical Res. 26, 221-234. Lande, R. 1976. Natural selection and random genetic drift in phenotypic evolution. Evolution 31, 314-334.

Lande, R. 1977. The influence of the mating system on the maintenance of genetic variability in polygenic characters. Genetics 86, 485-498. Lande, R. 1982. A quantitative genetic theory of life history evolution. Ecology 63,607-615. Ludwig, D. and S. A. Levin. 1991. Evolutionary stability of plant communities and the maintenance of multiple dispersal types. Theor. Population Biol. 40, 285-307. Metz, J. A. J., R. M. Nisbet and S. A. H. Geritz. 1992. How should we define “fitness” for general ecological scenarios? Trends in Ecology and Evolution 7, 198-202. Nagylaki, T. 1993. Evolution of multilocus systems under weak selection. Genetics 134, 627-647.

Orzack, S. H. 1993. Life history evolution and population dynamics in variable environments: some insights from stochastic demography. In Adaptation in Stochastic Environments, Lecture Notes in Biomathematics, J. Yoshimura and C. W. Clark (Eds) Vol. 98, pp. 63-104. Berlin: Springer-Verlag.

FLUCTUATING SELECTION: STRUCTURED POPULATIONS

127

Orzack, S. H. and S. Tuljapurkar. 1989. Population dynamics in variable environments. VII. The demography and evolution of iteroparity. Am. Naturalist 133, 901-923. Sasaki, A. and S. Ellner. 1995. The evolutionarily stable phenotype distribution in a random environment. Evolution (in press). Seger, J. and H. J. Brockmann. 1987. What is bet-hedging? &ford Sumeys in Evolutionary Biology 4, 182-211.

Templeton,

A. R. and D. A. Levin. 1979. Evolutionary

of seed pools. Am.

consequences

Naturalist 114, 232-249.

TuIjapurkar,

S. 1989. An uncertain

life: demography

in random

environments.

Population Biol. 35, 227-294. Tuljapurkar, S. 1990. Population Dynamics in Variable Enrironments. Biomathematics, Vol. 85. Berlin: Springer-Verlag.

Theor.

Lecture Notes in

Turelli, M. 1978. A reexamination of stability in randomly varying versus deterministic environments with comments on the theory of limiting similarity. Theor. Population Biol. 13, 244-267.

Turelli, M. 1988. Population genetic models for polygenic variation and evolution. In E. J. Eisen, M. M. Goodman, G. Namkoong and B. S. Weir (Eds), in The Second International Conference on Quantitatiue Genetics, pp. 601-618. Boston: Sinauer. Turelli, M. and N. H. Barton. 1990. Dynamics of polygenic characters under selection. Theor. Population Biol. 38, 1-57.

Yoshimura, J. and C. W. Clark (Eds) 1993. Adaptation in Stochastic Environments. Lecture Notes in Biomathematics, Vol. 98. Berlin: Springer-Verlag.

Received Received 3 January 1995 Revised version accepted 30 May 1995