Basic Model

We consider a diploid, randomly mating population in a seasonally fluctuating environment. While asymmetry in various model parameters will be explored later, we start with a fully symmetric model having a yearly cycle with g generations of winter followed by g generations of summer. The genome consists of L unlinked loci with two alleles each: one summer-favored and one winter-favored allele. For a given multilocus genotype, let ns and nw be the number of loci homozygous for the summer and winter allele, respectively, and nhet the number of heterozygous loci, with ns+nw+nhet=L.

In the basic model, loci are interchangeable in their effects (see Stochastic Simulations for a more general model), and the fitness of a multilocus genotype can be computed as a function of ns, nw and nhet. In the simplest case, fitness depends only on ns+0.5⋅nhet in summer and nw+0.5⋅nhet in winter, i.e., half the number of currently favored allele copies. To allow for dominance effects, we generalize this simple scenario and assume that fitness in summer depends on the summer scorezs:=ns+ds⋅nhet[1]and fitness in winter depends on the winter scorezw:=nw+dw⋅nhet.[2]The parameters ds and dw quantify the dominance of the currently favored allele in summer and winter, not with respect to fitness, but with respect to the seasonal scores zs and zw. Because we are interested in whether temporally fluctuating selection can maintain polymorphism in the absence of other stabilizing mechanisms, we only consider values of ds and dw between 0 and 1, and do not allow values >1, which would correspond to standard heterozygote advantage.

The relationship between the seasonal score z (z=zs in summer and z=zw in winter), and fitness, w, is given by a monotonically increasing fitness function w(z). This function specifies the strength of selection and accounts for epistasis. With discrete generations, the allele-frequency dynamics at a focal locus are driven by the relative fitnesses of the three possible genotypes at that locus—for example, the ratio of the fitness of homozygotes and heterozygotes. We say that there is no epistasis if these ratios and thus the strength of selection are independent of the number of other loci and their contributions to z. This is the case when fitness is multiplicative across loci:w=∏i=1Lwi⇔ln(w)=∑i=1Lln(wi),[3]where wi is the fitness value at locus i. In our model, this is achieved by setting w(z)=exp(z) because then Eq. 3 is fulfilled with wi=exp(1) if locus i is homozygous for the currently favored allele, wi=exp(ds) or wi=exp(dw) if it is heterozygous, and wi=exp(0)= 1 if it is homozygous for the currently disfavored allele. With v(z):=ln(w(z)), the multiplicative model is characterized by v″(z)= 0. We therefore use the second derivative of the logarithm of fitness v″(z) as a measure of epistasis (see ref. 39 for a similar definition of epistasis). Under positive or synergistic epistasis (v″(z)>0), the logarithm of fitness increases faster than linearly with z, and thus selection at a focal locus increases in strength with increasing contribution of the other loci to z. By contrast, under negative or diminishing-returns epistasis (v″(z)<0), selection at a focal locus becomes weaker with increasing contribution of other loci to z. We focus on two classes of fitness functions (Fig. 1). The first is of the formw(z)=(1+z)y[4]with a positive parameter y. Although for y> 1 in Eq. 4, fitness increases faster than linearly with increasing z (Fig. 1A), transformation to the logarithmic scale v(z)=y⋅ln(1+z) reveals that epistasis is negative for all y (Fig. 1B). Epistasis (v″) becomes more negative with increasing y, but v′ and thus the selective advantage of having an additional favored allele still increases with y for all z. Thus, larger y values correspond to overall stronger per-locus selection.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Examples for fitness functions generated by Eq. 4 or Eq. 5 with various parameters. In A, fitness is shown on a linear scale, and in B, on a logarithmic scale.

The second class of fitness functions is of the formw(z)=exp(zq)[5]with q≥1. This function reduces to the multiplicative model with q=1 (cyan lines in Fig. 1) and has positive epistasis with q>1 (e.g., magenta lines in Fig. 1).

In summary, fitness is computed in two steps. The first maps the multilocus genotype onto a seasonal score z to which loci contribute additively, essentially a generalized counter of the number of favored alleles (Eqs. 1 and 2), and the second maps z to fitness (Fig. 1). This two-step process disentangles dominance (step 1) from selection strength and epistasis (step 2).

In our model, genotypes with many summer alleles have a high summer score but a low winter score and vice versa, a form of antagonistic pleiotropy. Previous theoretical studies suggest that antagonistic pleiotropy is most likely to maintain polymorphism if for each trait affected by a locus the respective beneficial allele is dominant (40, 41). Such “reversal of dominance” also facilitates the maintenance of polymorphism in single-locus models for temporally fluctuating selection (42). Hypothesizing that reversal of dominance would also help to maintain polymorphism under multilocus temporally fluctuating selection, we assume that ds and dw in Eq. 1 and 2 take the same value, d, which means that dominance switches between seasons (see Stochastic Simulations for a more general model). For d<0.5, we have “deleterious reversal of dominance” (41), and the currently favored allele is always recessive, whereas for d>0.5, we have “beneficial reversal of dominance” (41), and the currently favored allele is always dominant. If d= 0.5, the seasonal score z is additive, not just between loci, but also within loci. Importantly, the value of d also determines the relative fitness of heterozygous intermediates, multilocus genotypes with the same number of summer and winter alleles and at least some heterozygous loci, compared with homozygous intermediates, which also have the same number of summer and winter alleles, but are fully homozygous (Fig. 2). For d<0.5, heterozygous intermediates have a lower seasonal score, z, and therefore a lower fitness in both seasons than homozygous intermediates. For d=0.5, heterozygous and homozygous intermediates have the same score and fitness. Finally, for d>0.5, heterozygous intermediates have a higher score and fitness. Interestingly, because d measures dominance not at the scale of fitness but at the scale of the seasonal score, z, beneficial reversal of dominance for fitness is neither sufficient nor necessary for d>0.5 (SI Appendix, Fig. S1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Values of the seasonal score, z, as a function of the dominance parameter, d, for two example four-locus genotypes: a heterozygous intermediate (blue line) and a homozygous intermediate (black solid line).

Multiple mechanisms could underlie a seasonal reversal of dominance. For example, metabolic control theory suggests that deleterious mutations affecting multistep metabolic pathways are generally recessive (43). If selection is fluctuating such that each allele is favored during one season and deleterious in the other season, we might thus expect a beneficial reversal of dominance. Alternatively, changes in dominance could be mediated by seasonal changes in gene expression. But even without changes in the genotype–phenotype map, seasonal changes in dominance are possible. In the example scenario in Fig. 3, the additive seasonal score, z, is a composite phenotype, a weighted average of two (also additive) traits—for example, starvation tolerance and fecundity—and there is antagonistic pleiotropy. Although the allelic effects on the two traits remain constant throughout the year, d>0.5 because the relative importance of the traits changes between seasons. This scenario requires changes (though not necessarily a reversal; SI Appendix, Fig. S2A) in dominance with respect to the pleiotropic effects of a locus on the two traits. For example, the winter allele (blue) in Fig. 3 produces higher starvation tolerance and is dominant for this trait, whereas it leads to smaller fecundity and is recessive there. One specific way in which such changes in dominance across pleiotropic effects can arise is via branched enzyme pathways with saturation or feedbacks (44).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Potential mechanistic underpinning for beneficial reversal of dominance. There is antagonistic pleiotropy for two traits, and the seasonal scores for winter and summer are computed as weighted averages of traits 1 and 2, with the relative importance of the two traits switching between seasons. The dashed line indicates the average of the two homozygote traits. If the heterozygotes are closer to the fitter homozygote with respect to both traits 1 and 2, there is a beneficial reversal of dominance at the level of the seasonal score, z. See SI Appendix, Fig. S2 for alternative scenarios.

Deterministic Analysis

In this section, we assume that population size is so large that genetic drift does not play a role. We also assume that mutations are rare enough that the allele-frequency dynamics will equilibrate before a new mutation arises at one of the L loci. This simple deterministic framework allows us to develop an intuitive understanding of the conditions for stable polymorphisms for various genotype-to-fitness maps. The intuitions developed here will then be checked and extended with stochastic simulations in the next section.

We will first confirm that the conditions under which seasonally fluctuating selection can maintain polymorphism are restrictive when contributions to the seasonal score z are additive within loci (d=ds=dw= 0.5 in Eqs. 1 and 2). Then, zs+zw=L for all possible genotypes, and the mean z over time is z∗=L/2. The long-term success of a genotype depends on its geometric mean fitness, or, equivalently, on the arithmetic mean of the logarithm of fitness, v(z). Jensen’s inequality or simple geometric considerations (Fig. 4) tells us that the arithmetic mean of v(zs) and v(zw) for a given genotype will be smaller than or equal to v(z∗) if v″<0 everywhere (Fig. 4A), equal to v(z∗) if v″=0 (Fig. 4B), and larger than or equal to v(z∗) if v″>0 (Fig. 4C).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

The logarithm of summer fitness (red) and winter fitness (blue) and the average logarithm of fitness (gray) as a function of a genotype’s summer score, zs. Assuming d=0.5, the winter score is zw=L−zs, leading to the mirror symmetry around L/2. If the fitness function is concave on the logarithmic scale, intermediate types have the highest average log-fitness (A); if the fitness function is log-linear, then all types have the same average log-fitness (B); and if the fitness function is convex on a logarithmic scale, extreme types have the highest average log-fitness and thus the highest geometric mean fitness (C).

The interannual allele-frequency dynamics (e.g., from summer to summer or from winter to winter) with multiplicative fitness (v″= 0) and d= 0.5 are thus neutral. No balancing selection emerges. With positive epistasis (v″>0), extreme types with either only summer or only winter alleles have the highest geometric mean fitness. Therefore, the population ends up in a state where all loci are fixed for the summer allele or all for the winter allele. With negative epistasis (v″< 0), the genotypes with the highest geometric mean fitness are those with the same number of summer and winter alleles and thus zs=zw. There are always some genotypes heterozygous at one or more loci that fulfill this condition (heterozygous intermediates; Fig. 2). For an even number of loci, zs=zw is also true for genotypes homozygous for the summer allele at half of the loci and homozygous for the winter allele at the other half (homozygous intermediates). When one of the homozygous intermediates fixes in the population, it cannot be invaded by any mutant starting at small frequency (under the assumptions of the deterministic model; see SI Appendix, section S1 for a detailed proof), and all polymorphism is eliminated. For an odd number of loci, homozygous intermediates do not exist, and some polymorphism may be maintained, at least at one locus. This case, which we will examine in more detail below, appears to be the only way in which seasonally fluctuating selection can maintain polymorphism under additivity (d=0.5).

Next, we explore whether deviations from additivity (d≠0.5) can facilitate multilocus polymorphism. Under multiplicative selection, i.e., without epistasis, the conditions for polymorphism at one locus are not affected by the dynamics at other loci. Thus, given the fitness values for individual loci (exp(d) for heterozygotes, exp(1) and exp(0) for currently favored and disfavored homozygotes, respectively, see text below Eq. 3), we can conclude that polymorphism is possible ifexp(d)2>exp(1)⋅exp(0)⇔d>0.5.[6]That is, there must be a reversal of dominance with respect to z, such that at any time the currently favored allele is dominant.

Now, we explore whether such a beneficial reversal of dominance can also maintain polymorphism in the presence of epistasis. In each case, a necessary condition for polymorphism is that a population fixed for the fittest possible fully homozygous genotype can be invaded by mutants. As we have seen above, with synergistic epistasis (v″>0), there are two fully homozygous genotypes with maximum fitness, the one with the summer allele at all loci and the one with the winter allele at all loci. In both cases, the resident type has score L in one season and score 0 in the other season, whereas mutants differing in one position have scores L−1+d and d. For mutants to invade, we thus needv(d)+v(L−1+d)>v(L)+v(0).[7]Using our example class of fitness functions with positive epistasis, Eq. 5 with q>1, we thus obtain the conditiondq+(L−1+d)q>Lq.[8]The critical value of d, dcrit, at which extreme types become invasible, satisfies(dcritL)q+(1+dcrit−1L)q=1.[9]For q=2, dcrit takes values 0.707, 0.954, and 0.995, with 1, 10, and 100 loci, respectively. For q>1 in general, dcrit approaches 1 as the number of loci increases. To see this, note first that the condition in Eq. 9 is always fulfilled for d=1 and thus dcrit≤1. Thus, as the number of loci, L, goes to infinity, (dcrit−1)/L becomes small, and we can approximate the second term on the left-hand side of Eq. 9 by a Taylor expansion around 1 to obtain(dcritL)q+1+q⋅dcrit−1L+O(1L2)=1.[10]Multiplying both sides by L and letting L go to infinity, we can conclude that dcrit is approximately 1 if the number of loci is large. Thus, for fitness functions of the type in Eq. 5 with positive epistasis, seasonally fluctuating selection can, in principle, maintain polymorphism at many loci, but the respective favored allele would have to be almost completely dominant, requiring large seasonal changes in dominance.

With diminishing-returns epistasis (v″<0), the fittest possible fully homozygous genotype carries the summer allele at half of the loci and the winter allele at the other half of the loci (assuming an even number of loci). Thus, a necessary condition for polymorphism is that this homozygous intermediate type can be invaded by mutants. The resident type has score L/2 in both seasons whereas mutants differing in one position have scores L/2+d and L/2−1+d. Thus, the resulting necessary condition for polymorphism isw(L/2+d)⋅w(L/2−1+d)>w(L/2)2.[11]Again, this condition is always fulfilled for d=1. For fitness functions of the form Eq. 4 with any exponent y, the critical dominance coefficient dcrit at which homozygous intermediates become invasible satisfies(1+L2+dcrit)(1+L2−1+dcrit)=(1+L2)2.[12]This quadratic equation has a negative solution, which is not relevant for our model, and a positive solutiondcrit=12(−1−L+(2+L)1+1(2+L)2).[13]From Eq. 13, dcrit decreases as L increases and approaches 0.5 as L goes to infinity. The intuition here is that the second derivative of the logarithm of fitness v″(z)=−y(1+z)−2 decreases with increasing z. Therefore, for large L, epistasis around the intermediate type with z=L/2 is weak, and the conditions for polymorphism approach those without epistasis. In other words, with increasing L, selection against temporal variation around the intermediate type becomes weaker, and a smaller change in dominance is sufficient to overcome it.

Our results so far suggest that for a broad class of fitness functions, seasonally fluctuating selection can maintain polymorphism if in both seasons the respective favored allele is sufficiently dominant. We call this mechanism “segregation lift” because it is based on a positive aspect of two alleles segregating at the same locus, as opposed to the negative aspect of segregation load. However, the preceding analysis does not tell us whether polymorphism will be maintained at all loci, or just one or a few of them. Also, it is still unclear how efficient segregation lift is at maintaining multilocus polymorphism in finite populations with genetic drift and recurrent mutations and whether genetic load is a problem. To address these questions, we now turn to stochastic simulations.

Stochastic Simulations

We use Wright–Fisher type individual-based forward simulations (see SI Appendix, section S2 for details). That is, for every individual in a generation independently, two individuals are sampled as parents in proportion to their fitnesses. We focus on diminishing-returns fitness functions of type Eq. 4 both because diminishing-returns epistasis appears to be more common and plausible (e.g., refs. 21⇓–23) and because the above theoretical arguments suggest that it is more conducive to multilocus polymorphism than synergistic epistasis. Specifically, the critical dominance parameter, dcrit, for diminishing-returns epistasis in Eq. 13 is generally smaller than the one for synergistic epistasis in Eq. 9. In addition, we run some simulations for the multiplicative model.

Additional parameters in the stochastic simulations are the symmetric mutation probability μ per allele copy per generation and the population size N. We generally keep population size constant, but also run supplementary simulations with seasonal changes in population size. Table 1 gives an overview of the model parameters and the ranges explored. In most simulated scenarios, selection and dominance effects are strong relative to mutation [w′(z) and d−0.5 are much larger than the mutation rate μ]. Although natural populations are often larger and mutation rates smaller than the values used here, many population genetic processes depend only on the product Nμ (e.g., ref. 45). Thus, large populations with small mutation rates may be well approximated by computationally more manageable smaller populations with larger mutation rates.

In addition to the basic model, we design a “capped” model to assess the relevance of genetic load. In this model, each individual can be drawn at most 10 times as a parent of individuals in the next generation, i.e., contribute at most 10 gametes. Once an individual has reached that number, its fitness is set to 0 so that it cannot be drawn again. To better understand the role of offspring-number capping, we also run supplementary simulations with a cap of three, the smallest possible value that still allows for differences in offspring number between individuals in the population.

From the simulation output, we estimate an “effective strength of balancing selection” (Materials and Methods and SI Appendix, section S2), which tells us whether and how fast a rare allele increases in frequency over a full yearly cycle. As expected from the above theoretical arguments, additive contributions within loci (d= 0.5) are not conducive to multilocus polymorphism (Fig. 5). For even numbers of loci, i.e., situations where both homozygous and heterozygous intermediates exist (Fig. 2), the effective strength of balancing selection estimated from the simulations is negative, indicating that rare alleles tend to become even rarer. For small odd numbers of loci, the effective strength of balancing selection is positive, but only one or two loci at a time fluctuate at intermediate frequency (SI Appendix, Fig. S4). As the number of loci increases, the effective strength of balancing selection eventually becomes negative even for odd numbers, decreases overall in absolute value, and finally approaches zero (effective neutrality) from below (Fig. 5). This behavior is independent of the exponent, y, of the fitness function Eq. 4. Also, as expected, effective balancing selection (Materials and Methods) emerges if the dominance parameter, d, is larger than a certain critical value, which decreases with the number of loci and is only weakly influenced by mutation rate (Fig. 6).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Effective strength of balancing selection (be in Eq. 14 in Materials and Methods) in the additive case (d=0.5) as a function of the number of loci. Solid lines indicate means and dashed lines indicate means ± two standard errors. Simulations are always run for successive odd and even numbers. N=1,000,g=15,μ=0.0001.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Critical value of the dominance parameter, dcrit, such that the effective strength of balancing selection (be in Eq. 14 in Materials and Methods), is positive (stable polymorphism) if d>dcrit, and negative (unstable polymorphism) if d<dcrit. Symbols represent means across replicates, and lines represent averages ± two standard errors. Since the pattern for odd numbers of loci is more complex (Fig. 5), only even values for the number of loci are included here. N=1,000,y=2,g=15.

From now on, we will focus on scenarios with large numbers of loci. For the case of 100 loci, Fig. 7 shows example allele-frequency trajectories for three different dominance parameters, d. For small d, each locus is almost fixed either for the summer or winter allele. For large d, all loci fluctuate at intermediate frequency. The critical dominance parameter, dcrit, with 100 loci is close to 0.5, independently of the exponent of the fitness function, y (Fig. 8A). For d<0.5, i.e., if the currently favored allele is recessive, the effective strength of balancing selection is negative and polymorphism is unstable (Fig. 8A). As d increases beyond 0.5, i.e., as the currently favored allele becomes more dominant, effective balancing selection becomes stronger (Fig. 8A). Both the stabilizing and destabilizing effects increase with increasing exponent y (Fig. 8A).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 7.

Three examples of allele-frequency trajectories for N= 1,000,L= 100,g= 15,y= 4,μ= 10−4, and d= 0.15 (A), d= 0.5 (B), and d= 0.65 (C). Only 10 randomly selected loci (shown in different colors) out of 100 loci are shown for 5 years (150 generations) in the middle of the simulation run (years 301–305).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 8.

Influence of the dominance parameter d on effective strength of balancing selection (A; be, Eq. 14, Materials and Methods), retardation factor (B), magnitude of fluctuations (C; se, Eq. 15, Materials and Methods), and predictability of fluctuations (D). Symbols indicate averages across replicates for the uncapped vs. capped model variant (often overlapping), and solid vs. dashed lines in A, C, and D indicate the respective means ± two standard errors. Lines in B simply connect maximum-likelihood estimates obtained jointly from all replicates. N= 1,000,L= 100,g= 15,μ= 10−4. See SI Appendix, Fig. S5 for more detailed information on the distribution and frequency dependence of seasonal allele frequency changes. The vertical gray lines are at d= 0.5. Bal. sel., balancing selection; exp., expected.

A tendency for rare alleles to increase in frequency does not guarantee that the average lifetime of a polymorphism is larger than under neutrality (42, 46). This is particularly interesting for fluctuating selection regimes with positive autocorrelation where alleles regularly go through periods of low frequency (42). We therefore compute the so-called retardation factor (46), the average lifetime of a polymorphism in the selection scenario relative to the average lifetime under neutrality (see SI Appendix, section S2 for detailed methodology). The results for 100 loci are consistent with those for the effective strength of balancing selection: For d>0.5, polymorphism under segregation lift is lost more slowly than under neutrality (Fig. 8B).

To quantify seasonal fluctuations, we compute an effective selection coefficient (Materials and Methods and SI Appendix, section S2). We also compute the predictability of fluctuations as the proportion of seasons over which the allele frequency changes in the expected direction, e.g., where the summer-favored allele increases over a summer season. Both the effective selection coefficient and the predictability of fluctuations have a maximum at intermediate values of d and increase with increasing exponent y of the fitness function Eq. 4 (Fig. 8 C and D). For even higher values of d, fluctuations are not as strong, presumably because heterozygotes are fitter, and therefore more copies of the currently disfavored allele enter the next generation. Also, effective strength of balancing selection, effective selection coefficient, and predictability of fluctuations increase with the number of generations per season (SI Appendix, Fig. S6).

With an offspring-number cap of 10, the results for the capped model generally match the results for the uncapped model in all respects, especially for d>0.5 (Fig. 8). It appears that the capping mechanism only rarely takes effect because of relatively narrow distributions of the seasonal score, z, within a generation (SI Appendix, Fig. S7), and consequent relatively low variance in fitness (SI Appendix, Fig. S8). For example, for d= 0.7 and y= 0.5, the fittest individual in the population was on average 1.2 times fitter than the least fit individual and 4.9 times fitter for y=4 (see SI Appendix, section S4 for a supporting heuristic analysis). When the offspring-number cap is set to three, substantial quantitative differences between uncapped and capped simulations are seen (SI Appendix, Fig. S9). However, 0.5 remains the critical dominance parameter. This result also holds for the multiplicative model (SI Appendix, Fig. S10), but with otherwise larger differences between capped and uncapped model versions, even with an offspring-number cap of 10. Also, the retardation factor for the multiplicative model is often below one, even when the effective strength of balancing selection is positive. The reason appears to be that there is larger variance in fitness under the multiplicative model (SI Appendix, Fig. S11) and that fluctuations are sometimes so large that alleles go to fixation (SI Appendix, Fig. S12). Both of these effects are weakened by offspring number capping, so that the capped multiplicative model behaves more similarly to the diminishing-returns model (SI Appendix, Fig. S10).

Additional simulations for the diminishing-returns model suggest that the finding of stable multilocus polymorphism for d>dcrit≈ 0.5 still holds under various forms of asymmetry, e.g., when one season has more generations than the other (SI Appendix, Fig. S13); when the exponent of the fitness function, y, differs between summer and winter (SI Appendix, Fig. S14); and when there are seasonal changes in population size (SI Appendix, Fig. S15). To explore the effects of asymmetry in the dominance parameters, we fixed the winter dominance parameter, dw, at 0.4 and varied the summer dominance parameter, ds. Stable polymorphism then arises for ds>0.6 (SI Appendix, Fig. S16), suggesting that stable polymorphism requires an arithmetic mean dominance parameter >0.5. Compared with the diminishing-returns model, the multiplicative model seems less robust to asymmetry (SI Appendix, Fig. S17).

With an increasing number of loci under the diminishing-returns model and with d>0.5, the strength of balancing selection, the retardation factor, and the magnitude and predictability of fluctuations all decrease (Fig. 9). Population size hardly influences effective strength of balancing selection and effective selection coefficient, measures which are based on average allele-frequency changes (Fig. 9 A and C), but large populations maintain polymorphism for longer (Fig. 9B) and have more predictable allele-frequency fluctuations (Fig. 9D). In small populations, polymorphism can even be lost slightly faster than under neutrality (Fig. 9B).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 9.

Influence of population size and the number of seasonally selected loci on effective strength of balancing selection (A; be in Eq. 14, Materials and Methods), retardation factor (B), magnitude of fluctuations (C; se in Eq. 15, Materials and Methods), and predictability of fluctuations (D). Symbols indicate averages across replicates and lines in A, C, and D indicate means ± two standard errors (in C and D, standard errors are too small to be visible). Lines in B simply connect maximum-likelihood estimates obtained jointly from all replicates. Note that in B, some points are missing because the rate of loss of polymorphism was too small to be quantified. d=0.7,y=4,g=15,μ=10−4. Bal. sel., balancing selection; exp., expected.

Finally, we consider a generalized model where parameters vary across loci and may be asymmetric between seasons. Independently, for each locus l, we draw four parameters: Summer effect size Δs,l and winter effect size Δw,l are drawn from a log-normal distribution. For this, we draw a pair of parameters from a bivariate normal distribution with mean 0, SD 1, and correlation coefficient 0.9, and then apply the exponential function to each of them. Summer and winter dominance parameters, ds,l and dw,l, are drawn independently from a uniform distribution on [0,1]. Seasonal scores are then computed as z=∑l=1Lcl, where the contribution cl of locus l in summer is 0 for winter–winter homozygotes, ds,lΔs,l for heterozygotes, and Δs,l for summer–summer homozygotes. Winter contributions are computed analogously. Because all effect sizes Δ are positive, all loci exhibit a trade-off between summer and winter effects. We use y=4 here because it led to the most stable polymorphism in the basic model.

The results indicate that polymorphisms with different parameters can be maintained in the same population, with their allele frequencies fluctuating on various trajectories (Fig. 10A). With a sufficiently high total number of loci, hundreds of stable polymorphisms (positive expected frequency change of a rare allele; see SI Appendix, section S2 for details) can be maintained in populations of biologically plausible size (Fig. 10B). The number of loci classified as stable depends only weakly on population size. However, only a small proportion of the polymorphisms classified as stable also exhibit detectable allele-frequency fluctuations, defined as changes in the expected direction by at least 5% in at least half of the seasons (Fig. 10 C and D). The number of detectable polymorphisms is highest at an intermediate total number of loci and increases with population size (Fig. 10D). Detectable polymorphisms tend to have larger summer and winter effect sizes than polymorphisms that are only stable (Fig. 10E). Compared with unstable polymorphisms, stable polymorphisms are more balanced in their summer and winter effect sizes (two-sample t-test on |ln(Δs,l/Δw,l)|, p<2.2⋅10−16; Fig. 10E; see also SI Appendix, Fig. S18). Many stable polymorphisms have asymmetric dominance parameters, but for almost all of them, detectable or not, the average dominance parameter across seasons is >0.5 (Fig. 10F).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 10.

Stability of polymorphism and detectability of allele-frequency fluctuations when parameters vary across loci and seasons. (A) Snapshot of allele-frequency trajectories for stable polymorphisms in one simulation run. (B) Average number of stable polymorphisms as a function of the total number of loci for different population sizes. (C and D) As in A and B, but only for polymorphisms that are also detectable. (E) Winter effect size, Δl,w, vs. summer effect size, Δl,s, for stable and detectable and only stable polymorphisms. The plot shows pooled results over ten simulation runs with independently drawn parameters. Oval isoclines indicate the shape of the original sampling distribution, with 75% of the sampling probability mass inside the outermost isocline. (F) Corresponding dominance parameters (see also SI Appendix, Fig. S19). The dominance parameters were originally drawn from a uniform distribution on the unit square. Parameters: y=4,g=10,μ=10−4 and in A, C, E, and F N=10,000,L=100.