Effect of migration and environmental heterogeneity on the maintenance of quantitative genetic variation: a simulation study

The paradox of high genetic variation observed in traits under stabilizing selection is a long‐standing problem in evolutionary theory, as mutation rates appear too low to explain observed levels of standing genetic variation under classic models of mutation–selection balance. Spatially or temporally heterogeneous environments can maintain more standing genetic variation within populations than homogeneous environments, but it is unclear whether such conditions can resolve the above discrepancy between theory and observation. Here, we use individual‐based simulations to explore the effect of various types of environmental heterogeneity on the maintenance of genetic variation (VA) for a quantitative trait under stabilizing selection. We find that VA is maximized at intermediate migration rates in spatially heterogeneous environments and that the observed patterns are robust to changes in population size. Spatial environmental heterogeneity increased variation by as much as 10‐fold over mutation–selection balance alone, whereas pure temporal environmental heterogeneity increased variance by only 45% at max. Our results show that some combinations of spatial heterogeneity and migration can maintain considerably more variation than mutation–selection balance, potentially reconciling the discrepancy between theoretical predictions and empirical observations. However, given the narrow regions of parameter space required for this effect, this is unlikely to provide a general explanation for the maintenance of variation. Nonetheless, our results suggest that habitat fragmentation may affect the maintenance of VA and thereby reduce the adaptive capacity of populations.


Introduction
As genetic variation is the fundamental basis upon which evolution acts, it is important to understand how variation is maintained in order to provide a foundation for answering various questions in biology and related fields, such as the 'missing heritability' problem in medical genetics (Maher, 2008), conservation of biodiversity (Cook & Sgr o, 2017) and population potential to respond to change (Houle, 1992). And yet, the relative importance of factors that influence variation and the mechanism(s) under which it is maintained are not well understood (Barton & Turelli, 1989;Mackay et al., 2009). Quantitative traits commonly experience stabilizing selection (Kingsolver et al., 2001), which in theory should erode genetic variation. However, high levels of standing variation and heritability of quantitative traits are consistently observed in nature (Johnson & Barton, 2005). This paradoxa high degree of genetic variation maintained in the face of stabilizing selectionremains a long-standing, unsolved problem in evolutionary biology and quantitative genetics theory.
The most widely studied explanation for this paradox is mutation-selection balance (henceforth referred to as MSB), the appeal of which lies in its intuitive logic: mutation, as the ultimate source of genetic variation, provides enough input to offset the eroding effect of selection, leading to a state of equilibrium. Under such models, stabilizing selection is assumed according to a Gaussian fitness function with parameter V S setting the . Multiple MSB models have been proposed, most notably the continuum-of-alleles model from Kimura & Crow (1964) of which two main approximations have been put forth: the Gaussian approximation [Kimura, 1965; later expanded by Lande (1976)] and the house-of-cards approximation (Turelli, 1984;B€ urger et al., 1989). The continuum-of-alleles model makes the basic assumption that for a locus with a continuous distribution of possible alleles, mutation can result in a new allele with an effect different from the pre-existing one. In general, the two approximations differ in the way each handles this assumption. For example, for an arbitrary diploid locus i with an allele of effect x i and a mutation with effect y i , under the Gaussian approximation the mutated allele would take on a new value which is conditional on the previous state (x i + y i ). This results in a predicted equilibrium genetic variance V G ðGÞ ffi 2n ffiffiffiffiffiffiffiffiffiffiffiffi ffi la 2 V S p , where n is the total number of loci that compose the quantitative trait, l the per locus mutation rate and a 2 the variance of the distribution of mutational effects. Under the houseof-cards approximation, x i may take on any value independent of the previous value, and mutation results in the replacement of x i by y i . In this case, the predicted variance at equilibrium, V G ðHCÞ ffi 4nlV S , is independent of the distribution of mutation effects. In contrast with the continuum-of-alleles model and its approximations, the diallelic model (Bulmer, 1972;Barton, 1986) assumes only two possible values at a locus with equal forward and backward mutation rates, such that mutation causes a flip between states (e.g. x i = {0,1}). In this case, the resulting equilibrium variance has been shown to be generally comparable to the house-ofcards approximation V G (HC) (Johnson & Barton, 2005).
Although these models have been extensively studied, conflicting evidence over their applicability in relation to realistic biological parameters has left debate open (Barton & Turelli, 1989;Johnson & Barton, 2005;Zhang & Hill, 2005). The assumptions and requirements of these models are unrealistic: primarily, they require extreme mutation rates or too many loci per trait (Johnson & Barton, 2005). This can be most simply illustrated by considering two empirical observations: heritabilities, h 2 , for the majority of traits are relatively high (h 2 = 0.2~0.6; Mousseau & Roff, 1987), and stabilizing selection is typically relatively strong in nature, as evidenced by the median value of estimates of reported quadratic selection gradient (c = À0.1; Kingsolver et al., 2001). The value c = À0.1 implies that the ratio of selection to phenotypic variation is V S /V P = 5 (Kingsolver et al., 2001;Johnson & Barton, 2005), which can be rearranged to account for typical estimates of heritability, yielding an expected empirical range of values for variation maintained in a natural population of approximately 0.04 < V A /V S < 0.12. Putting this in the context of Turelli's (1984) V G (HC) model under the assumption that V G % V A , then 0.01 < nl < 0.03. If, for example, there are 100 loci underlying a given trait, then this would require per-locus mutation rates on the order of~10 À4 , which is 10 to 100 times higher than most estimates from quantitative genetic studies, typically thought to be in the range of 10 À6 to 10 À5 (Barton & Turelli, 1989;B€ urger, 2000). Alternatively, given instead 1000 loci underlying a trait, then a per-locus mutation rate on the order of 10 À5 satisfies the bounds, and this discrepancy is not an issue. However, this alternative scenario is also not realistic, as pleiotropy puts limits on the number of traits that could be independently under selection, so very few traits could maintain variation in this way (Johnson & Barton, 2005).
More complex extensions of these models have been investigated (B€ urger et al., 1989;Zhang et al., , 2004B€ urger, 2014), namely by extending to multiple traits, such that 'apparent stabilizing selection' is generated through pleiotropic effects. Pure pleiotropy and joint effects models have been studied, as well as other extensions to include factors such as dominance or balancing selection. Nonetheless, these models have yet to provide a sufficient explanation for the patterns of variation maintained as suggested by empirical data (Johnson & Barton, 2005; but see Zhang & Hill, 2005). Another extension of mutation-selection balance models is to explore the interaction between spatial structure and genetic redundancy, which occurs when many different genotypes can yield the same phenotype (this is assumed to be operating when stabilizing selection occurs). When there is high genetic redundancy, drift can cause divergence among populations in the combinations of alleles that yield an optimal phenotype under uniform stabilizing selection, and migration can cause modest increases in variance, although the magnitude of this effect is also moderate and unlikely to provide a general explanation (Goldstein & Holsinger 1992;Phillips, 1996).
The above models all assume stabilizing selection to a constant environment, yet environments varying in space and time can also affect the maintenance of variation. If environmental heterogeneity can maintain sufficient differences in allele frequencies within a subdivided population, and migrant individuals introduce novel variants that can be maintained for some time, then migration can result in an increase in genetic variation within a population. For evolution in a quantitative trait, structured populations with limited amounts of gene flow have the potential to increase within-population variance (Slatkin, 1978;Lythgoe, 1997;Barton, 1999;Tufto, 2000;Spichtig & Kawecki, 2004;Huisman & Tufto, 2012;B€ urger, 2014), and a temporally fluctuating environment has been shown in some cases to increase variance under particular regions of parameter space (Kondrashov & Yampolsky, 1996;B€ urger & Gimelfarb, 2002;Gulisija & Kim, 2015). However, these studies have not explicitly framed results in terms of the relative increase in variance over mutation-selection to address the problem of the discrepancy between MSB predictions and empirical observations, and the effect(s) of time and space have not been investigated simultaneously for quantitative traits. Results about the maintenance of variation in spatial and temporally heterogeneous environments can also be derived from single-and two-locus population genetic models (Levins, 1964;Felsenstein, 1976;Hedrick, 1986Hedrick, , 2006Nagylaki & Lou, 2008), but such models do not give much insight into quantitative traits, due to the assumptions involved in extrapolating to many loci under phenotypic selection (as per Spichtig & Kawecki, 2004).
Here, we assess how different kinds of environmental heterogeneity affect the maintenance of genetic variation. We simulate spatial and temporal heterogeneity both independently and simultaneously to explore how the maintenance of variation is affected by variation in migration rate, population size, mutation rate, strength of selection and pattern of environmental variability (Table 1). In all cases, we explicitly focus on comparing heterogeneous vs. homogeneous environments to represent the increase in variance relative to that expected under mutation-selection balance. This provides some indication of how much more variance can be maintained under heterogeneous environments, which can be compared to the discrepancy between the mutationselection balance predictions and empirical observations described above.

Simulation set-up
A diploid Wright-Fisher two-patch scenario was modelled using the stochastic, individual-based simulation program Nemo (v2.3.45; Guillaume & Rougemont, 2006) under various parameter combinations of population size N (1000, 10 000), strength of selection V S (2, 5, 10), mutation rate l (10 À6 , 10 À5 , 10 À4 ) and backwards migration m (rate ranging from 10 À5 to 0.1 per generation) between two equal-sized patches with some local optima h. See Table 1 for detailed reference of each parameter and range of values investigated.
In each simulation, relative fitness is determined by Gaussian stabilizing selection acting on a single quantitative trait of 50 unlinked loci, each with an independent allelic value a (n = 50; r = 0.5). Nemo defines individual fitness of a quantitative trait according to the equation: where z is the phenotypic value of the individual, equal to Σ(a i,1 + a i,2 ) for all i diploid loci, h the local optimum and w 2 the strength of selection on the phenotype (V S = w 2 + V E ). As no effect of environmental variance was included (V E = 0), hereafter we use V S to represent the strength of selection on the genotype to simplify comparison to theoretical models. At initialization of a simulation, each locus is assigned an allelic value (a i ) drawn randomly from an exponential distribution (mean = 0.05). Subsequent mutation causes a change of (À1) Á a i , such that locus i has only two possible states, +a i or Àa i for the entire simulation (this follows a pseudo-diallelic model with a house-of-cards mutation scheme; note that the diallelic model precludes the gradual evolution of genetic architecture that tends to occur over long periods of time under a continuous mutation model; Yeaman & Whitlock, 2011). A different set of mutation effect sizes is drawn for each simulation replicate within a given parameter set, and the same set of effect sizes was then used for all subsequent parameter sets, to minimize the effect of stochastic outlier alleles on differences among parameter sets. Loci have additive effects on the phenotype, with no dominance or epistasis and free recombination between adjacent loci (r = 0.5). A total of six scenarios were run, each differing in the relation of local optima between patches: 1 Homogeneous patches (h 1 = h 2 = 1). The control set is used for comparison to the other scenarios; the local optimum is the same in each patch such that any differences that arise are due to drift rather than selection. 2 Pure spatial heterogeneity (h 1 = +1; h 2 = À1). This set of simulations introduced spatial environmental heterogeneity with opposite local optima between each patch. 3 Pure temporal heterogeneity (h t,1 = h t,2 ). To reflect a population inhabiting a changing environment, temporal fluctuations were modelled using an oscillating optimum defined by a sine wave centred about zero with amplitude of 1; hðtÞ ¼ sinð 2pt P Þ. Simulations were run with varying periodicity (P = 4, 16, 40, 100,  (2) and (3). Three subsets were used to represent differing degrees of spatial variation: (A, h t,1 = Àh t,2 ) opposing optima reflected about the horizontal axis; (B, h t,1 = h t,2 + 2) in phase temporal fluctuations with one patch oscillating about À1 and the other about +1; and (C, h t,1 = h t+(P/4),2 ) a partial phase shift between optima functions, such that h 2 lags behind h 1 (see Fig. S1 for illustration). All simulations were run for an initial period of 40 000 generations under these conditions, allowing the trait value to stabilize at an approximate equilibrium (assessed visually for the absence of change in trajectory of the trait value over time). Following selection and migration, a census of the population was taken every t generations (t = 100, or t = 5, for temporal simulations with P ≤ 100 generations) for the mean trait value and genetic variation (in terms of V A , the additive genetic variation) present in each patch. These data were averaged over the stable period of the last 10 000 generations (1000 for temporal simulations with P ≤ 100) for 40 replicates. For consistency, a 40950 matrix of allele values was generated such that each replicate had a different complement of mutation effect sizes, but that this complement was held constant across all simulation scenarios. Under the simulations that included temporal heterogeneity, preliminary results showed a cycling pattern emerge for both the trait value and V A over time; to use an estimate as representative as possible, the V A maintained within a cycle was sampled at 20 evenly spaced time points within a cycle, and the mean was calculated over the last ten cycles of the simulation. To obtain a measure of variation maintained within the simulated population that is comparable to empirical evidence, the ratio V A /V S was calculated from the final mean V A estimate of each simulation set. All data processing and analysis were performed with R v3.3 (R Development Core Team, 2016).

Results
Pure homogeneity, Set 1 (h 1 = h 2 = 1) We first describe patterns in the homogeneous environment, which acts as a control demonstrating the effect of finite population MSB processes in a two-patch model, to study the effects of mutation rate (l), selection strength (V S ) and population size (N), on the maintenance of variation (Fig. 1a). When compared to general MSB models (Turelli, 1984;B€ urger et al., 1989), the observed genetic variance at equilibrium is on the same order of magnitude as model predictions (See Fig. S2). As expected from quantitative genetic theory (Falconer & Mackay, 1996), variance was somewhat higher in larger populations, although the difference is minimal under most parameters (Fig. S3). Also following expectations, V A increases with l, with the effect of mutation rate approximately linear and independent of selection. Similarly, V A increases with V S (selection weakens), which qualitatively matches predictions of Turelli (1984) andB€ urger et al. (1989), although some quantitative deviations occur from the predictions, likely due to differences in the mutation model (see Figs S2.1 and S2.2 for simulated results under both a diallelic and continuous mutation model). Although from Fig. 1 there appears to be an interaction between selection and mutation rate, the relative effect of a 10fold increase in mutation (from l = 10 À5 to l = 10 À4 ) is similar for V S = 2 vs. V S = 10 (V A changes by a factor of 8.87 vs. 8.59 for V S = 2 vs. V S = 10; m = 0.1).
Notably, these patterns are qualitatively independent of N ( Fig. S3), and so the results below are presented for the subset of N = 1000 only. There is minimal response to rate of migration between the two patches, except for a small increase at low-to-intermediate m, which is consistent with predictions of Goldstein & Holsinger (1992), and arises due to an interaction between effects of genetic redundancy and genetic drift (hereafter, referred to as the 'redundancy-drift effect'; see Discussion).
Pure spatial heterogeneity, Set 2 (h 1 = +1; h 2 = À1) With spatial environmental heterogeneity, two key differences arise when compared to the control case (Set 1): an overall large increase in the amount of V A maintained, and the appearance of a threshold at high migration, above which the divergence among populations and amount of V A maintained both decrease [which is qualitatively consistent with the critical migration threshold predicted by Yeaman & Otto (2011)]. Although the migration rate at which V A peaks varies as a function of V S , the maximum V A reached under these parameters (ffi 0.32) is similar among levels of V S . Given the number of loci, mutation effect sizes and mutation rates that we used, spatial environmental heterogeneity generates up to a~10-fold increase in V A compared to the homogeneous case (for the same parameters, V A [homogeneous] = 0.033), as shown in Fig. 1b. This demonstrates that the redundancy-drift effect described previously by Goldstein & Holsinger (1992) is small relative to the effect of migration and spatial heterogeneity, at least under these conditions. Additionally, the effect of mutation becomes relatively less pronounced under spatial heterogeneity, as most of the variation is maintained by migration in this case. Whereas a 10-fold change in l (from l = 10 À5 to l = 10 À4 ) results in an increase by a factor of 8.87 in Fig. 1a, this increase only changes V A by a factor of 1.09 under spatial heterogeneity (V S = 2, m = 0.1).
Pure temporal heterogeneity, Set 3 (h t,1 = h t,2 ) In contrast to the case with spatial heterogeneity, under a temporally fluctuating optimum there is in general little change in V A in comparison with Set 1, the homogeneous case. As shown in Fig. 2a, a similar qualitative pattern is observed in response to migration rate, with an increase at low-to-intermediate m and rapid drop with m > 0.001, akin to the redundancy-drift effect observed in the homogeneous case. The resulting maximum peak in V A however is approximately 45% greater than the homogeneous case (V S = 5, l = 10 À4 , m = 2.5 9 10 À4 ), suggesting that temporal heterogeneity can only marginally increase the importance of this effect. Different periods (P) for the temporal oscillation of the environment were simulated, and complex interactions were found ( Fig. S4.1). Setting P ≤ 100 resulted in a response to migration similar to the homogeneous case and an overall increase in V A with l, as expected. However, the behaviour for P = 100 deviates under very strong selection (V S = 2), with the redundancydrift effect causing a peak in V A at lower migration rates than for shorter periodicities (see is minimal response to migration rate (i.e. no effect of redundancy-drift) and very little increase in variance when l = 10 À6 relative to simulations with shorter periods, but much more increases in variance under long periods at l = 10 À4 (see Fig. S4.1; V A under P = 1000 is 1.3-fold~2.5-fold higher than P = 100 under l = 10 À4 ). The degree of increase in V A under temporal heterogeneity relative to homogeneity therefore depends on mutation, period and selection, but is typically considerably less than under spatial heterogeneity with intermediate migration.
Combination of spatial and temporal heterogeneity, Set 4 For spatial and temporal heterogeneity combined (h t,1 = h t,2 + 2), the amount of V A and patterns that arise are qualitatively similar to those in the case of pure spatial heterogeneity. Similar threshold behaviour is observed, and similar levels of V A are maintained, peaking between migration rates of 0.01~0.1, again depending on V S (Fig. 2b). Under this form of combined spatial and temporal heterogeneity, there is a max 31% increase in V A compared to pure spatial heterogeneity (at V S = 5, l = 10 À4 , m = 0.001), which is very small when compared to the 9.75-fold difference in V A for the pure spatial heterogeneity vs. homogeneity under the same parameters V S and l.
In contrast, the complex patterns involving periodicity described above that were observed in the purely temporal set do not appear to carry over under combined spatial and temporal heterogeneity, likely because they are swamped out by the relatively larger effect of migration, which is relatively consistent across different period lengths (see Supplementary Materials, S4.2). There is an overall increase in V A maintained, which appears more similar to the case of pure spatial heterogeneity than pure temporal heterogeneity or homogeneity. This result is perhaps unexpected, as spatial and temporal effects of heterogeneity might be expected to work synergistically. Figure 3 shows a comparison among the six different patterns of environmental variation for a subset of the total parameter combinations used. For a given parameter set, there is a general trend of increased V A maintained with an increase in heterogeneity, particularly with the degree of spatial heterogeneity between patches. Temporal heterogeneity on its own does little to increase V A compared to the homogeneous set. When temporal heterogeneity is combined with spatial heterogeneity, the increase in V A is most pronounced in sets where the temporal fluctuations in the two patches are most out of phase ( Fig. 3; see Supplementary Materials, S1 for relation between h t,1 and h t,2 ). Of the three combinations investigated, subset 4Bwhere h t,1 6 ¼ h t,2 for any tmaintains a considerable higher level of V A compared to subset 4A, which in turn maintained more V A than in 4C, which was most in phase between patches. This further supports the emphasis on spatial heterogeneity over temporal heterogeneity as a stronger force for maintaining variation.

Comparing sets
To explore whether environmental heterogeneity could yield more realistic ratios of V A /V S (i.e. 0.04 < V A /V S < 0.12), we compared values for this ratio under homogeneous vs. heterogeneous simulation scenarios with all other parameters held constant (Fig. 4). Without any spatial heterogeneity, none of the parameter combinations result in V A that falls within this range for expected typical values (see Supplementary Materials, S5 and S6, for simulation sets that did not fall into range). However, with spatial heterogeneity, the combination of high migration (m ≥ 0.01) and strong selection (V S = 2 or 5) does result in V A within the typical range (i.e. an increase of~109 relative to MSB). For the combination of spatial and temporal heterogeneity, there is little change for most parameter combinations compared to spatial heterogeneity alone (Fig. 4b); for those where there is a noticeable deviation, whether the parameter set falls in or out of the target range remains unchanged (with an exception at V S = 2, m = 0.001).
Given migration between two patches with different optima, the genetic load observed due to deviation from the optimum (Fig. 5a) and due to genetic variance (Fig. 5b) increases considerably when m ≥ 0.01 (load, L ¼ WmaxÀ w Wmax ). This behaviour is expected, as the migrant individuals are introducing maladapted alleles, which increases both the variance and deviation from the optimum, contributing to load. As optimal local adaptation occurs when at least some alleles are strongly differentiated, this result is also reflected in the differences in allele frequency between patches, which decrease with increasing migration (See Fig. S7). For the homogeneous case, the observed load is minimal (L [homogeneous] ( L[heterogeneous]), and on the order of the genomic mutation rate, as expected from classical theory on mutation load (B€ urger, 2000; Fig. S2B). We note that these results were qualitatively indistinguishable for load calculated with W max = 1, because genetic variance is sufficiently high such that at least some nearly optimal genotypes are typically present in the population (See Fig. S8).

Discussion
The results of these simulations show that environmental heterogeneity can have a large impact on the maintenance of additive genetic variation (V A ) under a variety of parameter values. The direction and magnitude of this impact depend upon the rate of migration between subpopulations, with intermediate rates of dispersal providing the largest effect. This effect is also dependent on the type of environmental heterogeneityspatial, temporal or a combination of both all result in different magnitude and pattern of response to migration. Spatial heterogeneity has a much larger effect than temporal, and the effect of combined heterogeneity increases with the degree of spatial difference when temporal oscillations in each patch are more out of phase.

The redundancy-drift effect
For a range of migration rates (dependent on V S and l), spatial environmental heterogeneity substantially increases the variation a population can maintain. However, migration alone, even between two environmentally similar patches, can increase the variation maintained at very low rates of migration. This result is due to the interaction of genetic redundancy and the stochastic force of genetic drift, as described by Goldstein & Holsinger (1992). With enough population structure (i.e. limited gene flow), the genotypes existing in each patch are expected to undergo different histories of mutation and drift, resulting in relatively independent fluctuations in allele frequency within each patch, and possibly local fixation or loss. Over time, this results in a different pool of genotypes present in each patch, even if the mean phenotype is equivalent. Migration between patches with different allele frequencies will then increase the variance within patches; movement of individuals is adapted to a given environment with one complement of alleles that result in a similarly fit phenotype, such that it can be maintained in the new population; then, this flux of allelic variants/combinations should increase the variation maintained in the population. This effect is mainly dependent on two factors: (i) the number of possible redundant genotypes, and (ii) the degree of population structurepatches must be sufficiently connected such that new genotypes can be introduced by migration at a high enough rate to affect genetic variance, but not overly connected such that drift is no longer acting (relatively) independently in each patch. This genetic redundancy-drift effect is evident in the control case, where there is an observablealbeit minimalresponse to low migration, which persists regardless of population size, rate of mutation or strength of selection. This effect is somewhat increased under temporal heterogeneity, resulting in larger peaks over the same range of low dispersal as the homogeneous case, particularly under high mutation and strong selection (l = 0.0001, V S = 2). The increase of this effect may be due to the temporally variable selection within each patch causing different alleles to change in frequency more rapidly than migration homogenizes differences, resulting in a larger pool of differentiated loci than would be generated by drift alone.

Environmental heterogeneity
Nonetheless, this redundancy-drift effect is small relative to the effect of spatial heterogeneity. In this case, selection acts to pull each subpopulation towards its optima, resulting in different allele frequencies within in each patch. With migration, the alleles that would otherwise be purged by selection are persistently reintroduced into the other patch, increasing the V A maintained in the population. However, there is a threshold behaviour to this systema critical migration rate beyond which polymorphism can no longer be maintained (Felsenstein, 1976;Yeaman & Otto, 2011). Past this point, the allele frequencies of each patch are similar enough that migrants are much less likely to reintroduce novel/differentiated alleles, and the effect of migration on V A is reduced. Most populations experience some form of temporally variable selection, be it seasonal fluctuations or some otherwise unstable environment. Due to this, there has been much interest around fluctuating selection and its effect(s) on the maintenance of variation (Kondrashov & Yampolsky, 1996;B€ urger & Gimelfarb, 2002;Siepielski et al., 2009;Morrissey & Hadfield, 2012;Gulisija & Kim, 2015). Difficulties arise when attempting to draw comparisons across such studies due to a variety of different assumptions, such as the way in which the moving optimum is modelled and the nature of the underlying genetic architecture. A simple sinusoidal wave is commonly used, but varying parameters of this function (amplitude, centre, periodicity) appear to have complex interactions with each other as well as other parameters such as mutation rate and genetic redundancy. B€ urger & Gimelfarb (2002) found that in a diallelic multilocus model with recurrent mutation, the relative genetic variance increased with period length (typically ≥ 24 generations), but the magnitude of change and pattern of response to period depended on the amplitude of the oscillations. They observed the highest levels of genetic variance when the amplitude was set such that the optimum cycled between the most extreme possible genotypes, and that greater amplitudes, where the optimum fell beyond the maximum possible genotypic values, resulted in a decline in V A at high period lengths (≥ 52 generations). However, given the small number of loci (2-6) and scaling of allelic effects, this represents a scenario with highly nonredundant genetics. By contrast, Kondrashov & Yampolsky (1996) explored the maintenance of variation under temporal heterogeneity with high genetic redundancy and found that a considerable increase in variance (~3 orders of magnitude) was observed when the amplitude of the fluctuations of the optimum exceeded the width of the fitness function (amplitude, d > √[S 2 + V E ]). Although the assumptions they made about allele effect size and number of loci meant that the local optimum never exceeded the maximum/ minimum possible phenotypes, the large amplitude implied by the above inequality results in unnaturally large differences in fitness for a given phenotype over time. For example, an individual with a phenotype equal to +d would suffer a fitness cost of 86% when the optimum = Àd, when S = d and V E = 0. Similarly, for the standard set of parameters used in their simulations where substantial variance was maintained, the fitness cost for the same difference in phenotype exceeded 99%. By comparison, the fitness cost for the same difference in phenotype in our simulations was 63% (under strongest selection, V S = 2; this fitness cost is reduced to 18% under weaker selection, V S = 10). Common garden experiments typically show an overall magnitude of local adaptation of~45% (Hereford, 2009), so the selection regimes used by Kondrashov & Yampolsky (1996) constitute much stronger selection than normally observed. Thus, although they found that temporal heterogeneity could maintain much more variance than observed in our simulations, this was likely a consequence of drastic changes in fitness causing the underlying alleles to cycle rapidly in frequency, with increases in V A as the alleles reached intermediate frequency. Thus, the capacity for temporal heterogeneity to maintain substantial amounts of variance seems to depend upon either very extreme environmental fluctuations or a narrow range of rapid fluctuations with strict assumptions on the amount of genetic redundancy available. As empirical data become available for quantitative analyses of a population undergoing temporal variation, this will provide a realistic standard set of parameter values to be used in theory studies, allowing for a more quantitative assessment of the likely effects of temporal fluctuations and further evaluation of the relative importance of the above models.

Spatial vs. temporal vs. combined heterogeneity
Interestingly, it appears that spatial heterogeneity has a much stronger effect than temporal, both in isolation and when combined. This is evidenced by the difference in magnitude of V A maintained in Sets 2 vs. 3 (pure spatial and pure temporal heterogeneity, respectively), as well as the comparisons among the three scenarios of combined spatial and temporal heterogeneity (Set 4 A, B, C). The combination of spatial and temporal heterogeneity only marginally increased the variation maintained by spatial heterogeneity alone, and only in the case of complete spatial heterogeneity (Set 4B). Interestingly, among the three combined sets, V A increased with the spatial component (in general, , which was robust to changes in l and V S . Thus, our results indicate that spatial heterogeneity is more important for the maintenance of variation than a temporally fluctuating environment, at least for the genetic architectures and assumptions about fitness used here. Additionally, from the ratio V A /V S (Figs 4 and S5), temporal fluctuations alone fail to result in any simulation maintaining sufficient variation such that it is within the range 0.04 ≤ V A /V S ≤ 0.12, whereas the scenarios incorporating spatial heterogeneity do result in a shift towards this range. This indicates that it is insufficient on its own as a mechanism to explain the discrepancy between mutation-selection balance models and empirical observations. population should fall within the range 0.04 ≤ V A / V S ≤ 0.12, which would require either an unrealistically large number of loci or very high mutation rates under models of mutation-selection balance, as described in Introduction. To explore whether heterogeneity could result in more realistic ratios of V A /V S , we compared simulations run under homogeneous vs. heterogeneous conditions and examined the degree to which heterogeneity increased this ratio. The shift towards this window for simulations with spatial heterogeneity suggests that this is a viable mechanism for the maintenance of variation, but is limited to populations experiencing strong selection and high migration rates. As migration-selection balance is a spatially segregated form of negative frequency dependence, it seems likely that other mechanisms causing frequency dependence within a single population could also maintain considerable variation, such as sexually antagonistic selection (Turelli & Barton, 2004;Barson et al., 2015), nontransitive competition among genotypes (Lively & Sinervo, 1996) or fine-grained environmental heterogeneity (Sokolowski et al., 2007). This could be a particularly potent effect if strong selection acts on a few loci and indirectly protects linked loci from the expected erosion of variance under simple mutation-selection balance (e.g. Patten et al., 2010).
The results presented in Fig. 4 are intended to illustrate how the maintenance of realistic levels of V A is possible under environmental heterogeneity, rather than to represent the full span of parameter combinations that result in realistic behaviour. Little is known about the true values for many parameters involved, and changes in some values, such as the number of loci, are likely to alter the quantitative results. Rather, these results illustrate that for a given set of assumptions about the genetic basis of a trait, considerably more variation can be maintained by migration-mutation selection balance than mutation-selection balance alone. Moreover, the magnitude of this change can be sufficient to explain empirical observations in some cases, indicating environmental heterogeneity is a viable mechanism to consider for the maintenance of variation. However, we stress that given the sensitivity of this result to the balance between selection and migration and the narrowness of the region that maintains significantly more variance, it seems unlikely that environmental heterogeneity will generally resolve the 'paradox of variation'.
Although these results are qualitatively quite clear, comparison to established analytical models of mutation-selection balance is somewhat complicated by the differences in mutation scheme. The majority of theory and analytic models assume either diallelic loci of equal effect or quantitative traits with a normally distributed continuum of alleles. Our simulations attempt to combine the two varying allelic effects across loci (a i ), but with a single |a| per locus, and an exponential distribution of possible values a i . We chose this approach to reflect a (potentially) more biologically realistic scheme of mutation. For example, evidence suggests that the distribution of allelic effects at a locus is more likely to follow a gamma rather than a normal distribution, with a small number of loci of large effect and many loci of small effect (Orr, 2003). However, this leads to a less straightforward parameterization of V m or a 2 , which causes challenges when making comparisons to other analytical predictions (e.g. Supplementary Materials, Comparison of simulations to theory). Additionally, we note that altering the shape parameter of the gamma distribution did not substantially affect the maintenance of variation in our model, suggesting that the consequences of these assumptions are relatively small.
Given the relatively narrow range of parameter space over which we observe a significant effect of spatial heterogeneity on variance, it may often be difficult to evaluate whether natural populations commonly fall within this range, given the hurdles involved in quantifying selection and migration in nature (see Moore et al., 2007 for a comprehensive study). Evidence of migration maintaining quantitative genetic variation in spatially heterogeneous environments was found in lodgepole pine (Yeaman & Jarvis, 2006), but genomic signatures of this effect at individual loci may be difficult to detect. Presumably the increase in V A we observed was driven by increased heterozygosity at the underlying loci. However, given the difficulties involved in detecting signatures of local adaptation at individual loci using genome scans (De Mita et al., 2013;Lotterhos & Whitlock, 2015;Hoban et al., 2016), detecting the secondary effect of migration causing an increase in heterozygosity would likely be still more difficult. We did not explore the effect of migration-selection on heterozygosity at individual loci here, as a rigorous treatment would require extensive exploration of the effect of variation in the underlying genetic architecture and distribution of mutation effect sizes.

Implications
We have shown that spatial environmental heterogeneity and migration have the potential to increase genetic variance substantially, although temporal heterogeneity has a much more modest effect. As local adaptation is common (Hereford, 2009), this may be an important factor affecting evolvability. It has been noted that conservation management lacks proper consideration of evolutionary theory, implementation of which can improve overall efforts and long-term outcome (Cook & Sgr o, 2017). The described effect of migration suggests that the ecology underlying the maintenance of variation may be important to consider when planning conservation efforts with a focus on genetics.
For example, a major criterion for the IUCN (International Union for Conservation of Nature) Red List 3 1 ( 2 0 1 8 ) 1 3 8 6 -1 3 9 Benson et al., 2011). The genetic consequences of population size have been emphasized due to [entirely legitimate] concerns about reduced diversity, which is often used as a proxy measure for the 'health' of a population of concern in conservation biology. Here, we showed that a 10-fold change in population size had less effect on V A than a reduction in migration, which suggests that migration may also be an important factor to consider when assessing and managing the genetic health of populations of concern. Indeed, the related concepts of gene flow and connectivity have been previously suggested for consideration and implementation in conservation management, albeit with some contention, as too much migration can result in a reduction in variability (Weeks et al., 2011;Urban et al., 2012;Cook & Sgr o, 2017). Longterm maintenance of evolvability may depend in some cases on an intact and interconnected environment, rather than endogenous generation of variation (mutation-selection), which would depend less on environment. Therefore, a move towards improving connectivity between subpopulations and protection against habitat fragmentation may be a pertinent consideration for the maintenance of quantitative genetic variation, in addition to the well-recognized problems of reducing inbreeding and mitigating demographic fluctuations.
In conclusion, our results demonstrate a clear effect of migration on the maintenance of variation in all four investigated scenarios and highlight the potential for environmental heterogeneity to substantially increase V A . However, it seems unlikely that a single mechanism best explains the maintenance of variation that we see in nature; rather, the many concepts put forth over the past decades may be viable explanations in only certain scenarios, just as has been demonstrated here with environmental heterogeneity. Further investigation into potential mechanisms, particularly into those scenarios under which they succeed and fail, and how biologically relevant successful scenarios may be, is needed to get a complete picture of the fundamental process that is the maintenance of variation.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1 The change in local optima (h) over time for each patch in the three sub-scenarios (A, B, C) of simulations with combined spatial and temporal heterogeneity Figure S2.1 Comparison of genetic variance observed in simulations (filled) to predictions of mutationselection balance models (the house-of-cards approximation & the Stochastic house-of-cards) Figure S2.2. Comparison of genetic variance (V A ) maintained under a continuous mutation scheme and diallelic mutation scheme (as used for the main results) under homogeneous scenario (h 1 = h 2 ) and spatially heterogeneous patches (h 1 6 ¼ h 2 ) Figure S3 Effect of migration on the maintenance of genetic variation (V A ) for homogeneous patches (h 1 = h 2 ) with a population size of 1000 (solid line) or 10 000 (dashed line) for different mutation rates (l) and strength of selection (V S ) Figure S4.1 Effect of migration on the maintenance of genetic variation (V A ) for temporally heterogeneous patches Figure S4.2 Effect of migration on the maintenance of genetic variation (V A ) for patches with combined spatial and temporal (h t,1 = h t,2 + 2) heterogeneity for various periodicities Figure S5 V A /V S for simulations with temporal heterogeneity compared to simulations with homogeneous patches, for various rates of migration, mutation, and selection strength Figure S6.1 V A /V S for simulations with combined temporal and partial spatial heterogeneity (set 4A; h t,1 = Àh t,2 ) compared to simulations with homogeneous patches, for various rates of migration, mutation, and selection strength Figure S6.2 V A /V S for simulations with combined temporal and partial spatial heterogeneity (set 4C; h t,1 = h t+(P/4),2 ) compared to simulations with homogeneous patches, for various rates of migration, mutation, and selection strength Figure S7 Genetic Load due to (A) deviation from the optimum (B) genetic variance, assuming Wmax = 1 (i.e., the theoretical optimal phenotype) observed for spatial heterogeneity (solid) and homogeneous patches (broken) given mutation rate l = {10 À5 , 10 À4 }; V S = 5; N = 1000 Figure S8 Difference in frequencies of the '+' allele between patch 1 and patch 2 across the range of effect sizes for increasing rate of migration (Left to Right: m = 0.0001, 0.001, 0.01) Figure S9 Effect of ordering of lifecycle events on the observed genetic variance