Polygenic adaptation: From sweeps to subtle frequency shifts

Evolutionary theory has produced two conflicting paradigms for the adaptation of a polygenic trait. While population genetics views adaptation as a sequence of selective sweeps at single loci underlying the trait, quantitative genetics posits a collective response, where phenotypic adaptation results from subtle allele frequency shifts at many loci. Yet, a synthesis of these views is largely missing and the population genetic factors that favor each scenario are not well understood. Here, we study the architecture of adaptation of a binary polygenic trait (such as resistance) with negative epistasis among the loci of its basis. The genetic structure of this trait allows for a full range of potential architectures of adaptation, ranging from sweeps to small frequency shifts. By combining computer simulations and a newly devised analytical framework based on Yule branching processes, we gain a detailed understanding of the adaptation dynamics for this trait. Our key analytical result is an expression for the joint distribution of mutant alleles at the end of the adaptive phase. This distribution characterizes the polygenic pattern of adaptation at the underlying genotype when phenotypic adaptation has been accomplished. We find that a single compound parameter, the population-scaled background mutation rate Θbg, explains the main differences among these patterns. For a focal locus, Θbg measures the mutation rate at all redundant loci in its genetic background that offer alternative ways for adaptation. For adaptation starting from mutation-selection-drift balance, we observe different patterns in three parameter regions. Adaptation proceeds by sweeps for small Θbg ≲ 0.1, while small polygenic allele frequency shifts require large Θbg ≳ 100. In the large intermediate regime, we observe a heterogeneous pattern of partial sweeps at several interacting loci.


Author summary
It is still an open question how complex traits adapt to new selection pressures. 1 While population genetics champions the search for selective sweeps, quantitative 2 genetics proclaims adaptation via small concerted frequency shifts. To date the 3 empirical evidence of clear sweep signals is more scarce than expected, while 4 subtle shifts remain notoriously hard to detect. In the current study we develop 5 a theoretical framework to predict the expected adaptive architecture of a trait, 6 depending on parameters such as mutation rate, effective population size, size of 7 the trait basis, and the available genetic variability at the onset of selection. For 8 a population in mutation-selection-drift balance we find that adaptation proceeds 9 via complete or partial sweeps for a large set of parameter values. We predict 10 adaptation by small frequency shifts for two main cases. First, for traits with a 11 large mutational target size and high levels of genetic redundancy among loci, and genomic basis of polygenic traits. Consequently, following conceptual work by The fitness for individuals carrying 0, 1, 2, 3 . . . mutations (y-axis) are given for the complete redundancy (a) and relaxed redundancy (b) model of fitness effects, respectively. Grey balls show the fitness of ancestral wild-type individuals (without mutations). Colored balls represent individuals carrying at least one mutation, for time points t < 0 before the environmental change in blue and for t ≥ 0 in red. 161 For the models described above, we use Wright-Fisher simulations for a haploid, at locus i. We start our simulations with a population that is monomorphic for with the derivation of all results is provided in the Mathematical Appendix. 211 During the stochastic phase, we model the origin and spread of mutant copies 212 as a so-called Yule pure birth process following Etheridge et al. (2006) and Hermisson 213 and Pfaffelhuber (2008). The idea of this approach is that we only need to keep 214 track of mutations that found "immortal lineages", i.e. derived alleles that still 215 have surviving offspring at the time of observation (see Fig 2

Simulation model
Extended results including standing genetic variation and time-dependent fitness 226 are given in the Appendix. Assume now that there are currently {k 1 , . . . k L }, 0 ≤ 227 k j N e mutant lineages at the L loci. Then the probability that the next event in 228 the Yule process is either a birth (split) or a new mutation at locus i is 229 k i · p split + p mut,i L j=1 (k j · p split + p mut,j ) Importantly, all these transition probabilities among states of the Yule process 230 are constant in time and independent of the mutant fitness s b , which cancels in 231 the ratio of the rates. As the number of lineages at all loci increases, their joint 232 distribution (across replicate realizations of the Yule process) approaches a limit. 233 In particular, as shown in the Appendix, the joint distribution of frequency ratios 234 x i := k i /k 1 in the limit k 1 → ∞ is given by an inverted Dirichlet distribution Beta function and Γ(z) is the Gamma function. Note that Eq (5) depends only on 238 the locus mutation rates, but not on selection strength. 239 After the initial stochastic phase, when mutant lineages have established and 240 evaded stochastic loss, the dynamics can be adequately described by deterministic 241 selection equations. For allele frequencies p i at locus i, assuming linkage equilibrium, 242 we obtain (consult the Mathematical Appendix M.1 for detailed derivations) whereW andZ are population mean fitness and mean trait value. For the mutant 244 frequency ratios x i = p i /p 1 , we obtain 245ẋ i = d dt We thus conclude that the frequency ratios x i do not change during the deterministic 246 phase. In particular, this means that Eq (5)   The adaptive process is partitioned into two phases. The initial, stochastic phase describes the establishment of mutant alleles. Ignoring epistasis during this phase, it can be described by a Yule process (panel a), with two types of events (yellow box). Either a new mutation occurs and establishes with rate Θ l · s b or an existing mutant line splits into two daughter lines at rate s b . Mutations and splits can occur in parallel at all loci of the polygenic basis, (here 2 loci, shown in green and blue). Yellow and red stars at the blue locus indicate establishment of two recurrent mutations at this locus. When mutants have grown to larger frequencies, the adaptive process enters its second, deterministic phase, where drift can be ignored (panel b). During the deterministic phase, the trajectories of mutations at different loci constrain each other due to epistasis. We refer to the locus ending up at the highest frequency as the major locus (here in blue) and to all others as minor loci (here one in green).

266
While the joint distribution of allele frequencies provides comprehensive information 267 of the adaptive architecture, low-dimensional summary statistics of this distribution 268 are needed to describe and classify distinct types of polygenic adaptation. To 269 this end, we order loci according to their contribution to the adaptive response. 270 In particular, we call the locus with the largest allele frequency at the stopping  Wollstein and Stephan (2014). 279 Concerning our nomenclature, note, that the major and minor loci do not differ 280 in their effect size, as they are completely redundant. Still, the major locus is the 281 one with the largest contribution to the adaptive response and would yield the 282 strongest association in a GWAS case-control study. 283 In the following, we analyze adaptive trait architectures in three steps. In   291 For our biological question concerning the type of polygenic adaptation, the ratio 292 of allele frequency changes of minor over major loci is particularly useful. With 293 "sweeps at few loci", we expect large differences among loci, resulting in ratios 294 that deviate strongly from 1. In contrast, with "subtle shifts at many loci", allele 295 frequency shifts across loci should be similar, and ratios should range close to 296 1. Our theory (explained above) predicts that these ratios are the outcome of 297 the stochastic phase, yet their distribution is preserved during the deterministic 298 phase. They are thus independent of the precise time of observation. For our 299 results in this section, we assume that the mutation rate at all L loci is equal, 300 Θ i ≡ Θ l , for all 1 ≤ i ≤ L. This corresponds to the symmetric case that is most 301 favorable for a "small shift" scenario. 302 Consider first the case of L = 2 loci. There is then a single allele frequency 303 ratio "minor over major locus", which we denote by x. For two loci, the joint 304 distribution of frequency ratios from Eq (5) reduces to a beta-prime distribution. 305 Conditioning on the case that the first locus is the major locus (probability 1/2 for 306 the symmetric model), we obtain for 0 ≤ x ≤ 1, For small locus mutation rates Θ l , the order of the loci is largely determined 328 by the order at which mutations establish at these loci. I.e., the locus where 329 the first mutation establishes ends up as the major locus and the first minor 330 locus is usually the second locus where a mutation establishes. The distribution 331 of the allele frequency ratio x is primarily determined by the distribution of the 332 waiting time for this second mutation after establishment of the first mutation at 333 the major locus. In the 2-locus model, this time will be exponentially distributed, 334 with parameter 1/Θ l . In the L-locus model, however, where L − 1 loci with total 335 mutation rate Θ l (L − 1) compete for being the "first minor", the parameter for the 336 waiting-time distribution reduces to 1/(Θ l (L−1)). We thus see from this argument 337 that the decisive parameter is the cumulative background mutation rate

Expected allele frequency ratio
at all minor loci in the background of the major locus. In Fig 3 (orange dots) we 339 show simulations of a L = 10 locus model with an appropriately rescaled locus 340 mutation rate Θ l → Θ l /9, such that the background rate Θ bg is the same as for the 341 2-locus model. We see that the analytical prediction based on the 2-locus model 342 provides a good fit for the 10-locus model. A more detailed discussion of this type 343 of approximation is given in Appendix A.4. We contrast the expected allele frequency ratios of the first minor locus (with the second largest frequency) over the major locus (with the largest frequency) for 2 loci (blue dots) and for 10 loci (orange dots) with analytical predictions (Appendix, Eq M.16, black curve). E[x] is shown as a function of Θ bg (= Θ l for the 2-locus case). Panels correspond to different strengths of positive selection (s b , rows) and levels of SGV (no SGV, strongly deleterious s d = 0.1, weakly deleterious s d = 0.001, columns). We find that neither factor alters the expected ratio. We do not obtain results for all parameters, as we condition on adaptation from ancestral alleles, such that simulation runs are discarded if sampling conditions are met before the environmental change. Results for 10 000 replicates, standard errors < 0.005 (smaller than symbols).

345
While the distribution of allele frequency ratios, Eqs (5) and (9) 100 loci. Panels in the same row correspond to equal background mutation 357 rates Θ bg = (L − 1)Θ l , but note that the locus mutation rates Θ l are not equal.  across traits with different L is almost invariant. For large Θ bg , they converge 362 for sufficiently large L (e.g. for Θ bg = 10, going from L = 10 to L = 50 and to 363 L = 100). In particular, the background mutation rate Θ bg determines the shape 364 of the major-locus distribution (red in the Figure) for large p → 1 − f w = 0.95 (the 365 maximum possible frequency, given the stopping condition). For Θ bg < 1, this 366 distribution is sharply peaked with a singularity at p = 1 − f w , whereas it drops to 367 zero for large p if Θ bg > 1 (see also the analytical results below).

368
As predicted by the theory, Eq (8) for 0 ≤ p 1 ≤ 1 − f w (see also Appendix A.5). The distribution has a singularity at 405 p 1 = 0 if the corresponding locus mutation rate is smaller than one, Θ 1 < 1. It has 406 a singularity at p 1 = 1 − f w if the corresponding background mutation rate (which 407 is just the mutation rate at the other locus for L = 2) is smaller than one, Θ 2 < 1.

408
The marginal distributions at the major locus, P + fw [p|Θ 1 , Θ 2 ], and the minor locus, 409 P − fw [p|Θ 1 , Θ 2 ], follow from Eq (11) as This equation generalizes the definition of the background mutation rate, Eq (10) As long as Θ bg ≤ 1, we can approximate both the major-locus distribution P  Several partial and complete selective sweeps: Θ bg = 10 Frequency shifts: Θ bg = 100   x-axis: allele frequency Figure 4: Genomic architecture of polygenic adaptation. We distinguish three patterns of architectures with increasing genomic background mutation rate Θ bg : complete sweeps, for Θ bg 0.1, heterogeneous partial sweeps at several loci for 0.1 < Θ bg < 100, and polygenic frequency shifts for Θ bg 100. The plots show the marginal distributions of all loci, ordered according to their allele frequency, i.e. the major locus in red and all following (first, second, third, etc. minors) in blue to green. Lines in respective colors show analytical predictions, Appendix A.4. Simulations were stopped once the populations have adapted to 95% of the maximum mean fitness in each of 10 000 replicates, resulting in an the upper bound for the major locus distribution at, p 1 = 0.95. Simulations for s b = −s d = 0.1. Note the different scaling of the y-axis for different mutation rates.

Relaxing complete redundancy
Several partial and complete selective sweeps: Frequency shifts: Θ bg = 100 p< 3L_4 thB=100  p< 10L_11 thB=100 x-axis: allele frequency relaxed redundancy: 1 st and 2 nd major (= 1 st minor) locus 2 nd , 3 rd ... 9 th minor locus ≥10 th minor 1 st and 2 nd major locus (insets) complete redundancy: major and minor loci (simulations) Figure 5: Relaxed redundancy. Relaxing redundancy such that a single mutant has fitness 1 + 0.5s b/d and only two mutations or more confer the full fitness effect (1 + s b/d ) demonstrates the robustness of our model. As in Fig 4, allele frequency distributions of derived alleles are displayed once the population has reached 95% of maximum attainable mean population fitness. Genomic patterns of adaptation show similar characteristics as with complete redundancy. Due to relaxed redundancy, an additional "major locus" is required to reach the adaptive optimum. As explained in the main text, the distribution of the kth largest locus with complete redundancy therefore corresponds to the distribution of the k + 1st largest locus with relaxed redundancy. Insets in the second column show the same data with the distributions of the two major loci for relaxed redundancy combined (in green).  be derived from diffusion theory and has been known since the early days of 510 population genetics (S. ,   Wollstein and Stephan (2014), which rely on a gene-centered view to study the 549 pattern at a focal locus, irrespective of its role in trait adaptation. In contrast, 550 we use a trait-centered view, which is better suited to describe and distinguish 551 adaptive scenarios. For example, "adaptation by sweeps" refers to a scenario 552 where sweeps happen at some loci, rather than at a specific locus. This point is 553 further discussed in Appendix A.5, where we also display marginal distributions 554 of Eq (8) for fixed loci.

555
The role of the background mutation rate 556 Our results show that the qualitative pattern of polygenic adaptation is predicted 557 by a single compound parameter: the background mutation rate Θ bg (see Eqs (10), (13), (15)), 558 i.e., the population mutation rate for the background of a focal locus within the 559 trait basis. For a large basis, Θ bg is closely related to the trait mutation rate. 560 We can understand the key role of this parameter as follows. As detailed in is also independent of interaction effects due to epistasis or linkage. 576 Since the order of loci is not affected by the deterministic phase of the adaptive 577 process, Θ bg maintains its key role for the adaptive architecture. In the joint 578 frequency distribution, Eq (5) and Eq (8) The result also shows that the adaptive scenario does not depend directly on Hospital (2008)) is large. Second, the fitness function must exhibit strong negative 628 epistasis that allows for alternative ways to reach the trait optimum -and thus 629 produces redundancy (Gaussian stabilizing selection in Chevin and Hospital (2008)). 630 Finally, while the adaptive trajectory depends on the shape of the fitness function, 631 Chevin and Hospital note that it does not depend on the strength of selection on 632 the trait, as also found for our model.     Simulation results (colored dots) for the mean allele frequency ratio are plotted in dependence of the locus population mutation rate Θ l and compared with the analytical prediction (black line). Simulations are stopped when fitness has reached 95% of its maximum. Linkage does not change the results for the ratio of allele frequencies, despite significant build up of linkage disequilibrium with low recombination rates (data not shown). Results for 10 000 replicates standard errors < 0.005 (smaller than symbols).
if the base population is admixed, either due to natural processes or due to 1071 human activity (e.g. breeding from hybrids). For these scenarios, our theoretical 1072 formalism to describe the establishment of mutants during the stochastic phase 1073 (Fig 2) does not apply. In this section, we describe how the formalism can be 1074 extended to cover arbitrary starting frequencies of mutants at the onset of positive 1075 r=0 r=0.01 r=0.1 r=0.5          Extended Yule framework 1077 The Yule process that describes the stochastic phase of the adaptive process 1078 accounts for the mutant copies at all loci that are destined for establishment. In 1079 our framework so far (see the Mathematical Appendix M.2), we have started this process with zero copies. SGV due to mutation-selection-drift balance can still 1081 be produced by such a process if it is started at some time in the past (t < 0).

1082
For general starting frequencies, we can alternatively start this process at time 1083 t = 0, but with mutant copies (immortal lineages) already present. Suppose 1084 that the mutant frequency at locus i at time t = 0 is p i , corresponding to N e p i 1085 mutant copies. Of these, only the n i < N e p i "immortal" mutants (destined for 1086 establishment) are included in the Yule process. Assuming an independent establishment 1087 probability p est per copy, n i is binomially distributed with parameters N e p i and p est . 1088 For the limit distribution of a multi-type Yule process that is started with a non-zero To simplify the analysis, they do not model SGV as a distribution (due to mutation, 1103 selection, and drift), but replace this distribution by its expected value (ignoring 1104 drift). We can apply our scheme with fixed starting frequencies to this case and 1105 thus assess the effect of genetic drift in the starting allele frequency distribution. 1106 We assume equal loci and a starting frequency |µ l /s d | for an (initially deleterious) 1107 mutant allele with selection coefficient s d in the mutation-selection balance. Fig S.3   1108 shows the simulated marginal distributions of the loci with the largest contribution 1109 to the adaptive response (compare Fig 4). We see that the type of the adaptive 1110 architecture is again constant across rows with equal background mutation rate.   p< and p> ΘB=0,01; weak start frq:1; strong start frq:1 p< and p> ΘB=0,01; weak start frq:1; strong start frq:1 p< and p> ΘB=0,1; weak start frq:6; strong start frq:1 p< and p> ΘB=0,1; weak start frq:1; strong start frq:1 p< and p> ΘB=1; weak start frq:6; strong start frq:1   The panels show the adaptive architecture when mutant alleles start from their expected value in mutation-selection balance, without drift. We distribute L · |Θ l /2s d | mutant copies as evenly as possible across all loci. We set −s d = s b /100 = 0.001. Black lines for L = 2 loci show analytical predictions described in the main text (only computationally possible for Θ bg ≤ 1), green lines for Θ bg ≥ 1 show the heuristic prediction for Θ eff bg = 51Θ bg . Finally, gray lines show the marginal distributions when adaptation occurs from mutation-selection-drift balance, compare Fig 4.

1123
To extend our model to diploids, we assume that a single locus that is homozygous 1124 for the mutant allele is sufficient to produce the fully functional mutant phenotype, 1125 while a heterozygous locus produces a mutant that is functional with probability 1 − h. We assume that mutants contribute independently. Thus, if k heterozygous 1127 loci exist, but no homozygous mutant locus, the resulting mutant phenotype will 1128 be functional with probability 1−(1−(1−h)) k = 1−h k . For L = 2 loci, in particular, 1129 the (logarithmic) fitness of genotype G becomes In contrast to the haploid case, the marginal fitnesses are in general not equal.   Θ d l = 10   Θ d l = 100 For a k-locus model with equal mutation rate Θ      Θ bg = 10      where the time dependent coefficient s(t) defines the strength of directional selection. 11 We assume that s(t) < 0 for t < 0, but s(t) > 0 for t > 0, such that the optimal 12 trait value shifts from the wildtype state Z = 0 to the mutant state Z = 1 due to 13 some change in the environment at time t = 0. We also assume that selection is 14 stronger than drift, |N s(t)| 1 for almost all t, but is arbitrary otherwise. 15 We assume that Z is polygenic, with L biallelic loci (wildtype a i and mutant 16 allele A i , i = 1, . . . , L) constituting its genetic basis. While genotype a = (a 1 , a 2 , . . . , a L ) 17 produces the ancestral wildtype Z 0 , all mutant genotypes are fully redundant and 18 produce the mutant phenotype Z 1 , independently of the number of mutations. 19 New mutations from a i to A i occur at a rate µ i per generation, with µ i |s(t)| 20 for almost all t. For the purpose of our model, back mutation from A i to a i can 21 be ignored. The linkage map among loci is arbitrary -unless explicitly specified 22 otherwise. Let p i be the frequency of allele A i , and let f a be the frequency of the 23 wildtype genotype a. Then the mean fitness in the population is whereZ is the trait mean. Since W (Z 1 , t) = s(t)Z 1 is the marginal fitness of any 25 mutant allele, the selection dynamics at the ith locus can be expressed as Thus, the ratio of allele frequencies among loci does not change under selection. 33 differences in (relative) allele frequencies are due to mutation and drift. 35 We are interested in the pattern of allele frequency changes across loci during  Hermisson and Pennings (2005) for this concept). However, our analytical results 50 do not require a static equilibrium and, for a general s(t) < 0 for t < 0, the SGV 51 reflects this non-equilibrium dynamics. 52 As described in the main text, we dissect the adaptive process into two phases. 53 During an initial stochastic phase mutation, selection, and drift lead to the build-up 54 of genetic variation, either from SGV or due to new mutation after time t = 0, 55 as long as allele frequencies p i at all loci are still low. We will describe our 56 approach to this phase in detail in the section on Yule processes below. Once 57 allele frequencies are sufficiently large, genetic drift and recurrent new mutation 58 play only a minor role relative to selection until we reach the end of the rapid 59 adaptive phase. We thus enter a deterministic phase where the dynamics is then 60 well approximated by Eq (M.2b).

61
To relax the stringent redundancy condition of our model, it is natural to assume 63 that a single mutation is not sufficient to produce the full mutant phenotype Z 1 = 1, 64 but only a partial phenotype Z q = q with 0 < q < 1. This makes the marginal 65 fitness of mutant alleles dependent on the genetic background. If genotypes with 66 two or more mutations produce Z 1 , we have where f i is the frequency of the haplotype with a single mutation at locus i.
where f a is the frequency of the ancestral haplotype without mutations. 3. a genotype without a locus that is homozygous for the mutant, but with k loci that are heterozygous has fitness W (t) = s 2 (t) + 2 1−k s 1 (t) − s 2 (t) .
Since 2 where the approximation in the last term assumes that p est (t) 1, which is 153 usually the case unless selection is very strong. = k i p split + p mut,i L j=1 (k j p split + p mut,j ) .

(M.9)
Crucially, these transition probabilities are constant in time and independent of the 164 establishment probability p est (t). As a consequence, they are also independent of 165 the mutant fitness, which only affects the speed of the Yule process (via p est ), but 166 not its sequence of events. 167 We start the process with no mutants and stop it whenever the number of 168 mutants at one of the loci (e.g. locus 1) reaches some number k 1 = n. We are 169 interested in the distribution of the number of mutants k i at the other loci at this 170 time, respectively their ratios k i /n (remember that we already know that these 171 ratios stay invariant during the deterministic phase of the adaptation process). 172 We can prove the following 173 immortal mutant lineages across loci converges to the inverted Dirichlet distribution, where the vector Θ = (Θ 1 , . . . , Θ L ) summarizes the mutation rates and B .

(M.11)
Proof We proceed in three steps.

179
Step 1 Assume that we stop the process when the first locus reaches n > 180 0 lineages. We derive the probability that the process at this time is in state 181 (n, k 2 , . . . , k L ) as follows. We need n + k 2 + · · · + k L events ( 193 Step 2  . . .
because the integrand in this expression is just a Dirichlet density with shifted values of Θ i → Θ i + k i and the right hand side is the corresponding normalization factor. Then using L j=2 (Θ j ) (k j ) (Θ 1 + · · · + Θ L ) (n+k 2 +···+k L ) reduces (M.13) to (M.12). 200 The compound distribution Eq (M.13) can be interpreted as follows: If a random 201 experiment can have a finite number of outcomes (here: mutant lineages at one of 202 L loci), the negative multinomial distribution describes the probability to observe 203 each of these events k i times if we repeat the experiment until a focal event 204 (here: new mutant lineage at the first locus) has occurred n times. While the 205 negative multinomial distribution assumes that all outcomes occur with a fixed 206 probability y i , this probability is itself drawn from a Dirichlet distribution in the 207 Dirichlet-negative-multinomial compound distribution. In the present context, the 208 main advantage of (M.13) over (M.12) is that we can easily perform the limit 209 n → ∞ in this form.

210
Step 3 For large n → ∞, the values of k i /n, i ≥ 2, of the negative multinomial 211 distribution can be replaced by their expectations, 212 We can then transform the density (M.10) from variables y i to the x i (representing 213 the relative mutant frequencies). The entries of the Jacobian matrix (for 2 ≤ i, j ≤ 214 L) are Since this is the sum of an identity matrix (times a factor) and a matrix with 216 identical columns we can easily derive the eigenvalues and thus the determinant, n + k − 1 k y k (1 − y) n Γ(Θ 1 + Θ 2 ) Γ(Θ 1 )Γ(Θ 2 ) y Θ 2 −1 (1 − y) Θ 1 −1 dy and the inverted Dirichlet distribution (M.10) simplifies to a so-called β-prime 220 distribution, (M.14) If we measure the ratio x always relative to the locus with the higher frequency, 222 we obtain a conditioned distribution that is truncated at x = 1.  linkage. For f w > 0, in contrast, any constraint on the distribution of the p i due to 295 the stopping condition will necessarily also depend on the linkage disequilibria. 296 For further analytical progress we now assume that recombination is sufficiently 297 strong that linkage disequilibria can be ignored. We then obtain 298 L j=1 (1 − p j ) = f w (M.19) and the joint allele frequency distribution is given by the following Theorem, which 299 is our main analytical result.

300
Theorem 2 If the adaptive process is stopped at a frequency f w of the ancestral 301 phenotype in the population, and assuming linkage equilibrium among loci, the 302 joint distribution of mutant frequencies on the L-dimensional hypercube is where the δ-function restricts the support of P fw [{p i } i≥1 |Θ] to the (L−1)-dimensional 304 submanifold L j=1 (1 − p j ) = f w .

332
For the singular part, we thus have i dp 3 . . . dp L−1 .
Iterating the substitution procedure for variables p 3 to p L−1 , we arrive at demonstrating the singular behavior for p L → 0 and for p L → 1 − f w . Since 335 the labeling of loci is arbitrary, the assertion follows for all loci. Using this transformation on (M.10), the joint distribution of mutant frequencies (