The equilibrium allele frequency distribution for a population with reproductive skew

# The equilibrium allele frequency distribution for a population with reproductive skew

Foo University of Podunk, Timbuktoo    Barfootnotemark:    Baz Somewhere Else    Bof11footnotemark: 1    Hmmfootnotemark:    Ricky Der, Joshua B. Plotkin11footnotemark: 1 Department of Biology, University of Pennsylvania, Philadelphia, PA 19104
###### Abstract

We study the population genetics of two neutral alleles under reversible mutation in the processes, a population model that features a skewed offspring distribution. We describe the shape of the equilibrium allele frequency distribution as a function of the model parameters. We show that the mutation rates can be uniquely identified from the equilibrium distribution, but that the form of the offspring distribution itself cannot be uniquely identified. We also introduce an infinite-sites version of the process, and we use it to study how reproductive skew influences standing genetic diversity in a population. We derive asymptotic formulae for the expected number of segregating sizes as a function of sample size. We find that the Wright-Fisher model minimizes the equilibrium genetic diversity, for a given mutation rate and variance effective population size, compared to all other processes.

## 1 Introduction

Many questions in population genetics concern the role of demographic stochasticity in populations, and its interaction with mutation and selection in determining the fates of allelic types. The foundational work of Fisher, Wright, Haldane, Kimura [\citeauthoryearFisherFisher1958, \citeauthoryearWrightWright1931, \citeauthoryearHaldaneHaldane1932, \citeauthoryearKimuraKimura1994] and others has been instrumental in shaping our intuition about the powerful role that genetic drift plays in evolution, and especially its role in maintaining diversity. This classical theory, and the view of genetic drift as a strong force, emanates from the Wright-Fisher model of replication in a population, and its large-population limit, the Kimura diffusion [\citeauthoryearKimuraKimura1955]. The diffusion approximation has been particularly well-studied, not only because it is mathematically tractable, but also because it is robust to variation in many of the underlying model details. Many discrete population-genetic models, including a large number of Karlin-Taylor and Cannings processes [\citeauthoryearKarlin and McGregorKarlin and McGregor1964, \citeauthoryearCanningsCannings1974, \citeauthoryearEwensEwens2004], share the same diffusion limit as the Wright-Fisher model, and they therefore exhibit qualitatively similar behavior.

Nevertheless, Kimura’s classical diffusion is not appropriate in every circumstance. Its central assumption is the absence of skew in the reproduction process — that is, the assumption that no single individual can contribute a sizable proportion to the composition of the population in a single generation. Several recent studies have suggested that this assumption may be violated in several species, especially in marine taxa but also including many types of plants [\citeauthoryearBeckenbachBeckenbach1994, \citeauthoryearHedgecockHedgecock1994], whose mode of reproduction involves a heavy-tailed offspring distribution.

While the number of empirical studies on heavy-tailed offspring distributions is limited, there is a rich mathematical theory to describe the dynamics of populations with heavy reproductive skew. Beginning with Cannings’ 1974 paper on neutral exchangeable reproduction processes, this literature has led to generalized notions of genetic drift, which subsume the traditional Wright-Fisherian concept of drift. The resulting forward-time continuum limits of such processes generalize the Kimura diffusion. One tractable class of models are the so-called -Fleming-Viot processes, parameterized by a drift measure . The corresponding backward-time, or coalescent theory, for such processes leads to the -coalescents, first defined by Pittman and others [\citeauthoryearPitmanPitman1999, \citeauthoryearSagitovSagitov1999]. Two conspicuous features stand out in this more general theory: processes may have discontinuous sample paths, which feature “jumps” in the frequency of an allele, in contrast to the continuous sample paths of Kimura’s diffusion. Likewise, the coalescents of such processes typically exhibit multiple and even simultaneous mergers, instead of the strictly binary mergers of the classical Kingman coalescent.

Although mathematical aspects of such population processes (such as their construction, existence, uniqueness etc.) have already been described, the specific population-genetic consequences of reproductive skew have only recently begun to be worked out. In many cases, the classical picture of population genetics must be considerably enlarged to accommodate new phenomena — see for example [\citeauthoryearMöhleMöhle2006] on generalizations of the Ewens’ sampling formula, [\citeauthoryearEldon and WakeleyEldon and Wakeley2006b] on linkage disequilibrium in processes with skewed offspring distributions, [\citeauthoryearBirkner and BlathBirkner and Blath2008, \citeauthoryearBirkner, Blath, and SteinruckenBirkner et al.2011] on inference and sampling in the coalescent, and [\citeauthoryearDer, Epstein, and PlotkinDer et al.2012] on the fixation probability of an adaptive allele in the process.

The purpose of this paper is to study the stationary allele frequency distribution for populations with reproductive skew, under neutrality. When there are a finite number of allelic types subject to mutation, allele frequencies evolve to a unique stationary distribution, and our principle aim will be to understand how this distribution depends on the form of reproductive skew, , and how it may depart from the Wright-Fisherian picture.

Whereas a closed-form expression exists for the stationary allele frequency distribution in the (continuum) Wright-Fisher model, very few explicit expressions can be obtained in the general case of an arbitrary drift measure. Instead, we study the stationary distribution indirectly, first by deriving a recurrence relation satisfied by the moments in the two-allele scenario. This relation provides significant information about how the model parameters influence the stationary distribution. In particular, we demonstrate that the mutation parameters are identifiable from the stationary distribution, whereas the form of drift, , is not in general identifiable.

We also study how reproductive skew alters the standing genetic diversity in a population at equilibrium. Some numerical experiments of \citeNMohle2, as well as some asymptotic results of \citeNBerestycki1 for the Beta-coalescent, have suggested that the Wright-Fisher model tends to minimize standing diversity, compared to other offspring distributions. To analyze this behavior, we develop a -version of Kimura’s infinite-sites model, and we study the mean number of segregating sites, , in a sample of size . This measure of genetic diversity is robust in the sense that it is immune to many assumptions of the model and it coincides with the mean number of segregating sites in other infinite-sites models, including Watterson’s fully linked infinite-sites model. We demonstrate that the Wright-Fisher model minimizes diversity amongst all -processes of the same variance-effective population size. In other words, reproductive skew always tends to amplify standing genetic diversity, compared to the classical population-genetic model. We also derive a recursion formula for the mean number of segregating sites, and we use this to obtain asymptotic formulae for the number of segregating sites in large samples.

The remainder of the paper is structured as follows. We start by reviewing discrete population models under reproductive skew. We then describe the forward-time continuum limits of such processes, which can be identified as -Fleming-Viot processes. We develop a recursion equation for the moments of the stationary distribution in the two-allele case, and we use this to determine the identifiability of model parameters. To further examine equilibrium diversity we then introduce a two-allele, infinite-sites model with free recombination, and we study the frequency spectrum of samples from this process. This leads to a recursion formula for the mean number of segregating sites, and theorems concerning the minimization and maximization of diversity among all measures. We conclude by providing a simple intuition for our results, and by placing them in the context of the large literature on reproductive skew.

## 2 Discrete population models with reproductive skew

The -models are generalizations of the classical Wright-Fisher and Moran processes, which incorporate the possibility of large family sizes in the offspring distribution. The characteristics of these processes are most easily understood by studying their continuum limits, described below. Nonetheless, we shall first describe these models and review their properties in a discrete setting, along the lines of the treatment in \citeNEldonWakeley.

We consider a population containing a fixed number of individuals, each of two types. At every time step, a single individual is chosen uniformly from the population and produces a random number offspring, drawn from a distribution of offspring numbers, . The subsequent generation is then comprised of the offspring from the chosen individual supplemented by other individuals, randomly selected without replacement from the remainder of the population. Only a single individual contributes offspring in each reproduction event — the remaining individuals who neither contribute offspring nor die simply persist to the next time step.

When the offspring distribution is concentrated at two individuals, i.e. , this model coincides with the Moran process. More generally, we consider any discrete offspring number distribution supported on the set .

To incorporate mutation, an additional stage is appended after reproduction wherein each individual may mutate to the opposing type, independently and identically with a probability that depends upon the individual’s type, , with . This composite process is graphically depicted in Figure 1. We shall term this discrete process a “generalized” Eldon-Wakeley model.

### 2.1 The transition matrix.

We keep track of the number of individuals in generation of type 1, denoted by . Since is a Markov chain on the states , it possesses an associated transition matrix . The transition matrix can be written as a product of two matrices,

 P=QM (1)

corresponding to the reproduction and mutation stages described above. The rows of the mutation matrix are sums of independent binomial distributions, representing the mutational flux from each class. When all mutation rates are zero , the identity matrix. The matrix describes the neutral genetic drift due to reproduction alone in the -model; its form is more complex, with rows that are mixtures of hypergeometric distributions whose means depend on the offspring distribution . An important quantity is the first row of , called the “offspring distribution” of the process: for . The variance of this distribution, called the “offspring variance” , determines the time-scaling of the continuum limit (see below).

## 3 Continuum approximations of population models with reproductive skew

Analysis of the Moran or Wright-Fisher model is facilitated by a continuum limit, which becomes accurate in the limit of large population size [\citeauthoryearKimuraKimura1955, \citeauthoryearEwensEwens2004]. As described in [\citeauthoryearDer, Epstein, and PlotkinDer et al.2011], it is possible to derive a continuum limit for a significantly larger class of discrete population processes, including the models, without restrictions on the offspring distribution . Similar work in the Cannings case has been developed by \citeNMohle01. While the limiting continuum processes are not, in general, diffusions with continuous sample paths, they are still characterized by an operator , the infinitesimal generator of the continuum process, which reduces to the second-order differential equation of Kimura in the classical Wright-Fisher case.

Continuum approximations involve choosing how to scale time and space, as . Such scalings replace the number of individuals with the frequency , and the generation number by the time , for some choices of sequences . The continuum limit is then the process

 ~Xt=limN→∞1bNX[t/cN] (2)

In the classical Moran model we use the scalings and . In fact, it can be shown that the relationship between the space-scaling and time-scaling is fixed, in the sense that no other relationship leads to non-trivial limiting processes. We wish study allele frequencies, and hence impose the natural scaling . The general theory [\citeauthoryearMöhleMöhle2001, \citeauthoryearDer, Epstein, and PlotkinDer et al.2011] then indicates that the time-scaling must be proportional to

 cN=σ2NN (3)

where is the offspring variance of the process.

Once a time-scaling is fixed, then so is the appropriate scaling regime for the mutation rates, in order to produce a non-trivial balance of mutation and drift. This scaling must satisfy:

 μi=O(σ2N) (4)

In the classical Moran model , which produces the traditional scaling of mutation rate . In other models, such as the models of \citeNEldonWakeley, where and , mutation rates must scale faster to compensate for the increased rate of evolution from the drift process.

### 3.1 The limiting process for a generalized Eldon-Wakeley model.

By applying the techniques of [\citeauthoryearDerDer2010, \citeauthoryearMöhleMöhle2001], one may derive the continuum limit for the generalized Eldon-Wakeley process. These limits are characterized by an operator , and an associated Kolmogorov backward equation, analogous to the diffusion equation of Kimura. We assume that we have a sequence of Eldon-Wakeley models, one for each population size , and each with offspring distribution . We assume the time-scaling and mutational constraints of (3) and (4) so that

 θi=limN→∞2μiσ2N (5)

defines the effective population-wide mutation rate. Under an appropriate condition on the sequence of offspring distributions , there exists a limiting measure which may be derived from as:

 Λ =limN→∞ΛN (6) ΛN(i/N) =(iN)2P(N)U(i),i=0,…,N (7)

and which characterizes the continuum limit. Letting denote the time and state re-scaled process, one can show then that converges to a limiting process that satisfies the backward equation:

 ∂u(x,t)∂t=Gu(x,t),u(x,0)=f(x) (8)

where

 Gu(x)=12(−θ1x+θ2(1−x))∂u∂x+∫10xu(x+(1−x)λ)−u(x)+(1−x)u(x−λx)λ2dΛ(λ) (9)

and where .

The Markov process whose generator is given by (9) is called the forward-time, two-type -Fleming-Viot process.

### 3.2 Intuitive remarks on the generator.

As with the matrix decomposition of (1), the generator of (9) splits into two terms: a portion that describes mutation, independent of the reproduction measure , and an integral portion that describes genetic drift. The term describing mutation coincides with the standard first-order advection term in Kimura’s diffusion equation. The integral term however, generally differs from the Kimura term, and it depends on the drift measure .

Throughout the remainder of this paper we distinguish several important families of -processes. We define the pure processes to be those models for which , the Dirac measure concentrated at a single point , with . Since (9) expresses the generator as an integral decomposition over such Dirac measures, we can think a -process as being a random mixture of these pure processes. Of particular interest are the extreme cases , and — which correspond to the Wright-Fisher process and the so-called “star” processes, respectively. As we will show, these two processes constrain the range of dynamics in models. Another well-studied family in the coalescent literature are the Beta-processes, for which has a Beta distribution.

One can interpret as a “jump” measure controlling the frequency of large family sizes. If is concentrated near zero, then jump sizes are small. In this regime, the integrand behaves like the standard Kimura drift term . For the pure processes, where is concentrated at the point , allele frequencies remain constant for an exponential amount of time, until a bottleneck event in which a fraction of the population is replaced by a single individual. Such events cause the allele frequency to increase instantaneously by the amount , or decrease by . In the most general case of an arbitrary measure , these behaviors are mixed, and the jump events occur at exponential times with a random mean, and are associated with jumps of random size .

If places large mass near zero, the process becomes diffusion-like, with sample paths exhibiting frequent, small jumps. On the other hand, if is mostly concentrated away from zero then allele dynamics are of the “jump and hold” type, with fewer, but more sizable, jumps. Such behavior is most extreme in the star model, whose sample paths are constant until a single jump to absorption.

## 4 The Stationary Distribution of Λ-Processes

In the absence of mutation, , allele frequencies must eventually fix at 0 or 1, and thus any discrete generalized Eldon-Wakeley model possesses a trivial stationary distribution whose concentration at the absorbing states depends on the initial condition. When mutation rates are strictly positive, however, each generalized Eldon-Wakeley process in a population size possesses a unique, non-trivial stationary distribution, , to which the process converges, regardless of the initial condition.

In Figure 2, we plot the stationary distributions for a few -processes, in the case of symmetric mutation . Generally, these distributions have the same qualitative dependence on the mutation rate as the classical Wright-Fisher stationary distribution: they continuously progress from Dirac singularities at the boundaries to distributions concentrated more in the center of the interval, as the mutation rate increases. It is interesting to observe, however, that the non-Wright-Fisherian processes tend to have more mass at intermediate allele frequencies, and less relative mass near the boundaries, than the Wright-Fisherian model. We shall study this phenomenon more precisely, below.

The continuum limit of a sequence of Eldon-Wakeley processes also possesses a unique stationary distribution, . In the Appendix, we demonstrate that

 πN→π,as N→∞ (10)

In other words, the sequence of discrete equilibrium measures converges to the continuum equilibrium distribution. As a result, we can use the continuum equilibrium as a good approximation in large populations.

### 4.1 Moments of the Stationary Distribution.

In the case of two alleles, the stationary allele frequency distribution describes the likelihood of finding the mutant allele at any given frequency, at some time long in the future. We will study the moments of the stationary distributions for processes using a version of the Fokker-Planck equation, analogous to the equation used by Kimura to study the stationary distribution of the Wright-Fisher process. In general, a stationary distribution of a Markov process with generator is the solution to its so-called adjoint Fokker-Planck equation, so that

 ∫10Gu(x)dπ=0, (11)

for every smooth function on . We take as the generator for the process with mutation, given by (9). Although it is difficult to solve for in general, this equation can nonetheless be used to obtain detailed information about the stationary distribution.

To begin, we develop formulae for the moments of the stationary distribution, which will allow us to characterize aspects of standing genetic diversity. Let denote the -th moment of , . Setting into (11) yields an equation for the mean value of the equilibrium, so that

 m1=θ2θ1+θ2 (12)

Next, setting into (11) yields a relation between and ,

 m2=(1+θ2)θ2(θ1+θ2)(1+θ1+θ2) (13)

This recursive process can be continued, because the generator of (9) maps polynomials of degree to polynomials of degree . Thus, we can derive a system of equations that define the moments of . In the appendix, we show that this recursion has the form:

 mk=(k2θ2+ak−1,k)mk−1+∑k−2j=1ajkmjk2(θ1+θ2)+akk (14)

where the coefficients , , , are functions of , and are given by

 akk =∫101−(1−λ)k−kλ(1−λ)k−1λ2dΛ(λ) (15) aj,k =(kj−1)∫10λk−j−1(1−λ)j−1dΛ(λ),j=1,…,k−1 (16)

Initializing this system by (12), and observing that , we see that (14) uniquely determines the moments of , and indeed this equation can be used to solve for any specific moment of the stationary distribution. While it does not appear that the moments can be solved explicitly to produce simple, closed-form expressions as functions of the measure, it is clear that the coefficients are all linear combinations of moments of . Moreover, each moment is always a ratio of polynomials in and .

### 4.2 Identifiability of parameters from equilibrium.

One of the most important questions about the stationary allele frequency distribution is what population-genetic parameters can be identified from it — that is, which parameters of the population can be uniquely determined from data sampled in equilibrium? In the case of processes, the parameters we might wish to infer are the mutation rates, and , as well as the (high-dimensional) drift measure, , which describes the offspring distribution.

The first two moments of the stationary distribution are given by (12) and (13), and they are independent of the drift measure, . Thus, the first- and second-order moments of the stationary distribution for all -processes, and any function of these moments, such as the second-order heterozygosity , must coincide with those of the classic Wright-Fisher model. In the case of symmetric mutation , the third moment is also, remarkably, constant across the -processes, and has the value

 m3=2+θ4+8θ. (17)

The constancy of the first two moments with respect to , and the fact that the mapping from the first two moments to the two mutation parameters given by (12) and (13) is one-to-one, allows us to conclude that, regardless of the underlying reproductive process, the mutation rates are always identifiable from the equilibrium distribution. This is a tremendously productive result — because it means that we can always infer mutation rates from sampled data, even when the offspring distribution of a species is unknown to us.

Conversely, we may ask whether we can identify the form of reproductive process without knowledge of mutation rates — that is, is uniquely identifiable from the stationary distribution alone? It turns out that the answer is negative, as we demonstrate with the following simple example.

Consider the star process with mutation parameters . The generator for this process is

 Gu(x)=12θ(1−2x)u′(x)+(1−x)u(0)−u(x)+xu(1) (18)

The associated stationary distribution is easily derived (see Appendix), and it has a density given by

 dπ1dx=1θ|1−2x|1−θθ (19)

For comparison, the Kimura diffusion has a Dirichlet-type stationary distribution :

 dπ0dx=Γ(2θ)(Γ(θ))2xθ−1(1−x)θ−1 (20)

Note that both distributions (19) and (20) coincide when , despite the enormous difference between the drift measures of these processes. Hence, the map from a given process to its stationary distribution is not one-to-one, and consequently the drift measure cannot generally be identified from the stationary distribution. This is a pessimistic result, because it implies that that the offspring distribution and form of genetic drift cannot in general be inferred from data collected in equilibrium, even when the mutation rate is known.

## 5 An infinite-sites model for the Λ-processes

In order to understand how reproductive skew influences standing genetic diversity, we now develop an infinite-sites version of the process and study its equilibrium behavior. This model generalizes the infinite-sites approach of [\citeauthoryearRoyChoudhury and WakeleyRoyChoudhury and Wakeley2010, \citeauthoryearDesai and PlotkinDesai and Plotkin2008], for the Wright-Fisher model. We will study the sampled site frequency spectrum of our model, under two-way mutation. Our analysis will allow us to quantify our previous observation that the Wright-Fisher model minimizes the amount of standing genetic diversity, amongst all processes. The site frequency spectrum that we will describe in this section, for independent sites, differs from the Watterson-type spectrum for fully linked sites; but our approach nonetheless yields information in that case as well.

We consider an evolving population of large size , following the reproduction dynamics of a neutral forward-time -process, for a fixed measure. We keep track of sites along the genome, each with two possible allelic types under symmetric two-way mutation at rates . The allele dynamics at each site are described by a two-type -process; and the site processes are assumed independent of one another (that is, we assume free recombination).

Let denote the two-allele stationary distribution for the model, given by (11), where the subscript denotes the explicit dependence on the mutation rate. We imagine sampling individuals from the population at equilibrium, assuming . We let , represent the (random) number of sampled individuals at site with a particular allelic type, so that their joint distribution has the form

 P(Y1=y1,…,YL=yL)=L∏i=1∫10(nyi)xyi(1−x)n−yidπθ(x) (21)

The sampled site frequency spectrum [\citeauthoryearSawyer and HartlSawyer and Hartl1992, \citeauthoryearBustamante, Wakeley, Sawyer, and HartlBustamante et al.2001] is defined as the vector

 Zk=L∑i=11Yi=k,k=0,…,n (22)

The variables record the number of sites with precisely (out of ) sampled individuals of a given allelic type. In this sense, the sampled site frequency spectrum represents a discretized version of the stationary distribution . The variables are distributed multinomially on the simplex . The sites are called the segregating sites, representing locations where there is diversity observed in the the sample. Conversely, the sum represents the number of monomorphic sites in the sample.

### 5.1 The infinite-site limit and its Poisson representation.

To study the sampled site frequency spectrum we take the limit of an infinite number of sites, , and we apply a Poisson approximation. We define the genome-wide mutation rate as , and we assume that this mutation rate approaches a constant in the limit of many sites: . In the Appendix, we show that the segregating site variables then converge, as , to a sequence of independent Poisson random variables with means , given by

 cj(n)=limθ→01θ∫10(nj)xj(1−x)n−jdπθ(x),j=1,…,n−1 (23)

The numbers may be interpreted as an infinite-sites sample frequency spectrum. From (23), it is apparent that the means depend on the heterozygotic moments of , and thus, also on the moments of . This representation is thus a generalization of a result of \citeNArindam for the two-allele Wright-Fisher independent-sites model, where , and where the spectrum has the form for .

For the Wright-Fisher model, \citeNArindam has shown that the number of segregating sites in the sample of size , , is a sufficient statistic for , under the independent-sites assumption. This is an important result because the number of segregating sites vastly compresses the information in the frequency spectrum, yet nonetheless contains no loss of information for the purposes of inferring the mutation rate. The Poisson representation of the sample frequency spectrum we have derived shows that remains Poisson distributed even in the general infinite-sites case — under the assumption of site independence. Thus, the sufficiency of for remains true, and consequently possesses desirable qualities for robust estimation of .

### 5.2 Diversity amplification and the number of segregating sites.

The number of segregating sites in a sample is a classic and powerful method to quantify genetic diversity in a population. Here we study how depends on the form of reproduction — that is, on the form of the drift measures . In particular, we will show that the Wright-Fisher model minimizes the expected number of segregating sites in a sample, compared to all other processes. Thus, large family sizes in the offspring distribution will tend to amplify the amount of diversity in a population.

Under the infinite-sites Poisson approximation, the number of segregating sites in a sample of size is Poisson-distributed, and its expected value is

 ESn=En−1∑j=1Zj=Θn−1∑j=1cj(n) (24)

where are the coefficients in (23). The binomial theorem applied to (23) then shows that:

 ESn=limθ↓01θ∫10[1−xn−(1−x)n]dπθ(x) (25)

and so we may interpret the expected number of segregating sites as a type of higher-order heterozygosity statistic of the stationary distribution. According to (23), is a linear combination of moments of , of order at most . It follows that the average number of segregating sites may be evaluated by the recursion (14) and it can be expressed as rational functions of moments of . The first several such expressions are listed below:

 ES2 =Θ (26) ES3 =32⋅Θ (27) ES4 =∫10(5λ2−14λ+11)dΛ(λ)∫10(6−8λ+3λ2)dΛ(λ)⋅Θ (28) ES5 =52∫10(5−6λ+2λ2)dΛ(λ)∫10(6−8λ+3λ2)dΛ(λ)⋅Θ (29) ES6 =12∫10(2608λ2−1558λ+411−2428λ3+1312λ4−388λ5+49λ6)dΛ(λ)∫10(6−8λ+3λ2)(15−40λ+45λ2−24λ3+5λ4)dΛ(λ)⋅Θ (30) ⋮

These expressions for the expected number of segregating sites become extremely complex for larger sample sizes . Nevertheless, we can use asymptotic methods to study how diversity is expected to behave in large sample sizes. We will address two primary questions. First, how does does the expected number of segregating sites, , grow as a function of the sample size, for a given drift-measure ? And second, which reproduction processes maximize and minimize , for fixed ?

In the appendix, we use the moment recursion (14) to derive the following recursion for the sequence

 ESn=nΘ2ann+n−1∑j=1ajnannESj (31)

with the initial value , and where are given by (15) and (16). We can use this relation to obtain detailed information about both as a function of the sample size , and as a function of the underlying measure.

Consider first the pure processes , in which a single individual may replace a given fixed fraction of the population. Then we can prove from (31) that (see Appendix):

 ESn=λΘ2n+O(logn),n→∞ (32)

Two features are of interest in this asymptotic expression. First, the average number of segregating sites grows linearly with sample size in a pure process, as opposed to logarithmically as in the Wright-Fisher case. Second, the rate of linear growth depends on the jump fraction, , so that asymptotically, diversity is maximized for large replacement fractions , and correspondingly minimized when this fraction is small.

Equation (32) can be generalized to a larger class of measures. If is any probability measure whose support excludes a neighborhood of zero, then we have the asymptotic formula:

 ESn=C(Λ)Θ⋅n+O(logn),n→∞ (33)

where . This equation shows that linear growth of the expected number of segregating sites is characteristic of any process whose drift measure is bounded away from zero — that is, any process that does not contain a component of the Wright-Fisher process. This result allows us to determine which reproduction processes maximize and minimize the average number of segregating sites , for a given value of . In the appendix, we prove the following optimization principle: for each sample size , the diversity maximizing and minimizing processes within the class of all processes must in fact be pure processes, i.e. where is concentrated at a single point. It follows then from (32) that, asymptotically, the Wright-Fisher model () minimizes, and the star-model () maximizes, respectively, the mean number of segregating sites amongst all -processes.

Although these results apply in the limit of large sample sizes, we conjecture that the Wright-Fisher and star models are also the extremal diversity processes for any sample size, . From the optimization principle stated above, it suffices to check this statement within the restricted class of pure processes. In Figure 3, we show as a function of the jump-size parameter for the pure models, for a few values of . These results confirm that the Wright-Fisher model minimizes , whereas the star model maximizes , over all -models. We have conducted numerical studies which support this proposition more generally, even for very small sample sizes. In this sense, the Wright-Fisher model and star models are extremal processes, and, for a given effective variance population size, respectively minimize and maximize the expected genetic diversity in any sample.

## 6 Discussion

We have studied the stationary distribution of a very general class of population models, under mutation. We have focused on understanding the interaction between the form of the offspring distribution and resulting form of genetic drift it engenders, as well as the shape of the stationary distribution. We have demonstrated that the mutation rate can always be uniquely identified from the stationary distribution, even when the drift measure is unknown (as it always will be, in practice). However, the form of the drift measure cannot always be uniquely identified from equilibrium properties of the process — and so may require dynamic data to determine its specific form.

The stationary allele frequency distribution of the Wright-Fisher process is extremal, in a sense, within the class of -processes. Specifically, the Wright-Fisher model exhibits greater probability mass near very high and low allele frequencies. This observation was formalized by analyzing a infinite-sites model, in which we found that the mean number of segregating sites in a sample is indeed minimized by the Wright-Fisher process.

Our results can be placed in the context of a nascent literature that views the Wright-Fisher process as an extremal model within the large space of possible population processes. Aside from the diversity-minimization property we have demonstrated here, it has previously been observed, for instance, that among -processes, the Wright-Fisher model minimizes the fixation probability of an adaptive allele [\citeauthoryearDer, Epstein, and PlotkinDer et al.2012], minimizes the time to absorption for new mutants [\citeauthoryearDer, Epstein, and PlotkinDer et al.2011], and, among generalized coalescents, possesses the fastest rate of “coming down from infinity” [\citeauthoryearBerestycki, Berestycki, and LimicBerestycki et al.2010]. The basic intuition behind all these results revolves around the type of sample paths possessed by different processes. A typical sample path in the Kimura diffusion undergoes a high frequency of small jumps (in fact, is continuous), and thus new mutants persist for only generations before being eliminated by genetic drift. By contrast, in a general model with the same variance effective population size, large jumps in the sample path may occur, but with lower frequency, thereby lengthening the absorption time — for example, up to order generations in the pure processes. Since the mean number of segregating sites in the entire population is the product of the genomic mutation rate and the expected absorption time for a new mutant, standing genetic diversity must increase when reproductive skew is present.

Although we have presented results only within the class of -processes, many of our formulae — for example (14) – can be generalized to the set of all Cannings models. We expect the diversity-minimization property of the Wright-Fisher model will hold even within this larger family.

The infinite-sites model of the -process we have developed here differs from the Watterson infinite-sites model typically encountered in coalescent theory, in two respects. First, we have assumed free recombination and hence independent sites, whereas in Watterson’s model sites are tightly linked. Second, we assume two-way mutation between alternative alleles at each site, whereas Watterson’s model features one-way mutation at each site away from the existing type. Nonetheless, some of the results derived for our site-independent, infinite-sites model extend to the Watterson, linked infinite-sites processes as well.

In general, the (random) number of segregating sites in a sample is a function of the dependency structure among sites. For example, in the simple Wright-Fisherian case, independence of sites gives rise to a Poisson distribution for , compared to a sum of geometric random variables in the case of no recombination [\citeauthoryearEwensEwens2004]. However, as \citeNWatterson75 has already remarked, the mean value of is generally robust to the recombination structure of an infinite-sites model. If denote the allelic distributions at sites, then (22) shows that the expected number of segregating sites is a function only of the marginal distributions of , instead of their joint distribution. Thus the expected number of segregating sites in a sample is unaffected by linkage. Likewise, the distinction between one-way and two-way mutation (and folded and unfolded spectra) does not alter the mean number of segregating sites other than by a possible overall scaling.

Because is such a common measure of genetic diversity, our results have some connections to the literature on -coalescents. Recently, \citeNBerestycki12arxiv showed that, for those measures whose coalescent comes down from infinity, the (random) number of segregating sites in a sample of size for the Watterson model has the asymptotic law:

 Sn∫n0qψ−1(q)dq→Θ (34)

where is the Laplace exponent of the measure, defined as:

 ψ(q)=∫10exp(−qλ)−1+qλλ2dΛ (35)

The authors conjectured that (34) holds more generally, even when does not come down from infinity. In this respect, our asymptotic result (33) for — derived for measures bounded away from zero (and thus always fail to come down from infinity) is evidence in favor of their more general conjecture, in the case not covered by the hypotheses of their theorem. For under such assumptions, the Laplace exponent has the expression

 ψ(q)∼q∫10λ−1dΛ (36)

which implies from (34) that

 ∫n0qψ−1(q)dq∼n⋅(∫10λ−1dΛ)−1 (37)

which is proportional with our formula (33) for . Finally, returning to the case of independent sites, developed in this paper, it is also true that the distributional convergence of (34) holds, a fact which follows from the Poisson representation for .

In our analysis of the expected number of segregating sites, we have concentrated on the two extreme cases— the Wright-Fisher case, for which is known to grow logarithmically in the sample size , and the case of pure processes (and more generally those processes whose drift measure support excludes zero), for which we have demonstrated linear growth of . Nevertheless, the recurrence relation (31) can be used to analyze intermediary cases as well, for example the Beta processes, in which the density of behaves like a power-law in the vicinity of zero. For such reproduction measures, growth in diversity with sample size will lie somewhere between the logarithmic and the linear cases.

## 7 Appendix

### 7.1 The Stationary Distribution of Processes with Reproductive Skew.

Let be a sequence of discrete generalized Eldon-Wakeley processes, one for each population size , converging to a continuum process , under the state and time re-normalization of (3), (4). If we suppose that each discrete process operates under strictly positive mutation rates , then it is easily verified that the associated forward-time transition matrices possess strictly positive entries, and thus, from the Perron-Frobenius theorem, there exists a unique stationary distribution for each process . A standard argument, using the fact that the sequence is tight, shows that there is a subsequence converging to a probability measure which is a stationary distribution for (see \citeNEthier+Kurtz86, for example). This argument indeed demonstrates that any weak limit point of is a stationary distribution ; below, through the moment recursion, this distribution is uniquely characterized, and hence every weakly convergent subsequence of converges to , thus .

### 7.2 Derivation of a Recursion for the Moments of the Stationary Distribution.

Let be the stationary distribution for the two-type process, which satisfies (11). In this section we obtain a recursion formula for the moments of .

Define the operator . Setting , we have

 Lxk =((1−λ)k−1)xk−(1−λ)kxk+1+x(x+(1−x)λ)k =((1−λ)k−1)xk−(1−λ)kxk+1+xk+1(1−λ)k+kxk(1−λ)k−1λ+k−2∑j=0(kj)xj+1(1−λ)jλk−j =((1−λ)k−1+kλ(1−λ)k−1)xk+k−1∑j=1(kj−1)xj(1−λ)j−1λk−j+1 (38)

Separating the generator (9) into the mutation and pure-drift portions, we define the latter to be the operator

 GDu(x)=∫101λ2Lu(x)dΛ(λ). (39)

If we write

 GDxk=k∑j=1bjkxj (40)

Then substituting (38) into (39), and then comparing the coefficients to (40) , we have

 bkk =∫10(1−λ)k+kλ(1−λ)k−1−1λ2dΛ(λ) (41) bjk =(kj−1)∫10λk−j−1(1−λ)j−1dΛ(λ),j

Let be any stationary distribution of the process. Then according to (11),

 ∫10[12(−θ1x+θ2(1−x))kxk−1+GDxk]dπ(x)=0 (43)

Using the expansion for above, we derive:

 (bkk−k2(θ1+θ2))mk+(k2θ2+bk−1,k)mk−1+k−2∑j=1bjkmj=0 (44)

where is the -th moment of the . This is equivalent to formulae (14), (15) and (16), where , and for .

### 7.3 Derivation of the star-process stationary distribution.

Consider the probability measure on , with density

 dμdx=1θ|1−2x|1−θθ (45)

The generator for the star-process undergoing symmetric mutation is , and the space of twice continuously differentiable functions is a core for . Noting that the density satisfies the equation and integrating by parts, it is readily verified that for every . Thus is a stationary distribution for the process, and is further unique as established by the moment recursion (44).

### 7.4 Poisson representation of the infinite-sites model.

In this section we show that the segregating site variables in the independent sites model converge to a sequence of Poisson random variables with means , given by (23). We make use of the structure of the moments of the stationary distribution as found in (14). First we require a preliminary lemma.

###### Lemma 1.

Let be the stationary distribution of a process undergoing symmetric mutation . Then there exist constants , for ,

 pj(θ)=∫10(nj)xj(1−x)n−jdπθ(x)=cjθ+o(θ),θ↓0 (46)
###### Proof.

The lemma is equivalent to saying and has a derivative at . First observe from (14) that under , all the moments of the stationary distribution are differentiable in for all , and hence is differentiable everywhere. Also

 |pj(θ)|≤(nj)∫10x(1−x)dπθ(x)=(nj)2θ1+2θ (47)

hence as . ∎

Now to establish the Poisson representation, it is enough to apply the well-known Poisson approximation to the multinomial distribution. We use:

###### Theorem 1.

[\citeauthoryearMcDonaldMcDonald1980]. If is multinomial with parameters , and are independent Poissons with means , then

 ∥(Z1,…,Zn−1)−(V1,…,Vn−1)∥≤2L(n−1∑j=1pj)2 (48)

where is the total variation norm of measures.

The Poisson representation is now obvious, since by Lemma 1 and therefore the right-hand side of (48) goes to zero as . Since are converging to a sequence of independent Poisson distributions with finite means , where are as in Lemma 1, so must .

### 7.5 The average number of segregating sites.

In this section, we study the number of segregating sites in the infinite-sites model, deriving the recursion (31) for the average diversity measure , and use it to obtain asymptotic expressions for diversity.

Let be the stationary distribution of the process under two way symmetric mutation . Define as the heterozygosity measure

 Hn(θ)=n−1∑j=1∫10(nj)xj(1−x)n−jdπθ(x) (49)

Applying the binomial theorem,

 Hn(θ)=∫10(1−xn−(1−x)n)dπθ(x)=∫10(1−2xn)dπθ(x) (50)

where the second equality follows from symmetry of the stationary distribution. Now, define the diversity measure . We have from (25):

 sn=limθ↓01θHn(θ) (51)

Under symmetric mutation , the recursion formulae for moments (14) reads

 mn(θ)=nθ+∑n−1j=1ajnmj(θ)nθ+ann (52)

From symmetry of the stationary distribution and differentiability of , there are numbers such that . Inserting this into the right-hand side of (52), and expanding in a Taylor series, we obtain, by comparing the first-order coefficients, a recursion for :

 vn=−n/4+∑n−1j=1ajnvjann (53)

The equations (50) and (51) imply that . Thus the corresponding recursion for is

 sn=n/2+∑n−1j=1ajnsjann (54)

where we initialize . Observe that by the binomial theorem, one has the relation

 ann=n−1∑j=1ajn (55)

Thus the numbers define a probability measure on the set . By studying this measure and the recurrence relation defining , we may derive the asymptotics for .

Now suppose that the underlying process is associated with a measure with support bounded away from zero. Then from (15), , for some . Therefore, . Using this estimate in the recursion (51), and defining ,

 sn=An/2+An−1∑j=1ajnsj+O(γn) (56)

Now, we shall prove

###### Theorem 2.
 sn=Cn+O(logn),n→∞ (57)

where

 C=A2−2A∫10λ−2(1−λ)dΛ (58)
###### Proof.

Let . We find, using the explicit expressions (16),

 Cn+gn=An/2+An−1∑j=1∫10λ−2(nj−1)λn−(j−1)(1−λ)j−1dΛ⋅(Cj+gj)+O(γn) (59)

The right-hand side of the above is

 An/2+An+1∑j=1∫10Cjλ−2(nj−1)λn−(j−1)(1−λ)j−1dΛ+An−1∑j=1∫10gjλ−2(nj−1)λn−(j−1)(1−λ)j−1dΛ+O(γn) (60)

where we have subtracted two binomial terms and exponentially bounded them. Using the formula for the mean of a binomial random variable with parameters in the first series, this is equivalent to

 An/2+An∫10C(1−λ)λ−2dΛ+An−1∑j=1∫10gjλ−2(nj−1)λn−(j−1)(1−λ)j−1dΛ+O(1) (61)

Inserting this expression into (59), dividing by , and using the definition of , we find that satisfies the relation

 gn=An−1∑j=1∫10gjλ−2(nj−1)λn−(j−1)(1−λ)j−1dΛ+O(1/n) (62)

Completing the series with two binomial terms once more,

 gn=A∫10λ−2(Eλ,ng)dΛ+O(1/n) (63)

where the expectation occurs with respect to a binomial random variable with parameters . We may further interpret the integral term as an expectation of under a mixture of binomials with mixing weights .

We can now prove that that . To be explicit in (63), let be a constant such that

 (64)

for all . We shall show, by induction, that , for large . By enlargening , the base case is satisfied. Assume that the statement it is true up to . Then using Jensen’s inequality on the mixture in (63),

 gn+1≤Blog(A∫10λ−2(1−λ)ndΛ)+B/n (65)

If is a number smaller than , then we can, keeping in mind the definition of , further bound this expression by:

 gn+1 ≤Blog((1−ϵ)n)+B/n (66) =Blog((1−ϵ)n)+Blog(n+1)−Blog(n)+O(1/n2) (67) ≤Blog(n+1) (68)

Thus Theorem 2 is proved. ∎

The constant in Theorem 2 is algebraically equivalent to in (33).

### 7.6 An Optimization Principle for the average number of segregating sites.

In this section we prove the following theorem:

###### Theorem 3.

For every sample size , and for fixed mutation rate , the minimum and maximum values of over the class of processes is achieved within the class of pure Eldon-Wakeley processes; that is, where , for .

###### Proof.

It is evident from the recursion for in (54) that the average number of segregating sites must have the form:

 ESn/Θ=∫10fn(λ)dΛ∫10gn(λ)dΛ (69)

for some functions (in fact, polynomials) . Because is positive for every pure Eldon-Wakeley process (where ), we can without loss of generality assume that . For such functions, we have the lemma:

###### Lemma 2.

For positive functions defined on , we have the inequalities:

 minλf(λ)g(λ)≤∫10f(λ)dΛ∫10g(λ)dΛ≤maxλf(λ)g(λ) (70)
###### Proof.

We prove the lower bound; the upper bound is established in the same way. Say is a two-point measure: . Then elementary manipulations show

 p1f(λ1)+p2f(λ2)p1g(λ1)+p2g(λ2)≥min {f(λ1)/g(λ1),f(λ2)/g(λ2)} (71)

By induction one easily generalize to measures concentrated at any finite number of points. Finally, the full case is obtained by taking weak limits of measures concentrated at a finite number of points. ∎

Returning to the proof of the main theorem, one sees that Lemma 2 immediately implies the result, since is precisely for the case . ∎

The optimization principle can be refined. Let be any probability measure whose support excludes a neighborhood of zero, and suppose that and are the smallest and largest values, respectively, on which is supported, i.e. , and . Then

 λmin2≤C(Λ)≤λmax2 (72)

In conjunction with (33), this shows that the asymptotic growth rate in for a -process whose drift measure is supported on the interval can be lower and upper-bounded by the rates of growth in of two pure processes, with parameters