Total consensus under high reproductive-variance conditions

Star-shaped branching patterns of genealogies are common in marine species. High-fecundity marine populations are characterized by low ratios of effective to actual population size, which reflect high variance in reproductive success among parents in mass spawns. When extreme reproduction events occur, offspring from very few parents dominate the population (whereby multiple mergers, or subsets of lineages with star-like trees, arise) and thus, the loss of genetic diversity is significant. Under high reproductive-variance conditions (assuming that reproduction occurs by sampling from the Pareto distribution), this paper explores the distribution of heterozygosity across generations. The result shows that zero heterozygosity is not achieved, implying that the populations may decline without evident loss of genetic variation. It is also found that there are singularities in the heterozygosity distributions. However, in the case of high reproductive variance, the locations of the singular points subtly deviate from those of the case where reproduction occurs by Wright-Fisher sampling.


Introduction
Genetic drift is changes in allele frequencies due to stochasticity (or variance) in the individuals' reproductive success within a population (Wright 1931). Conventionally, the effect on genetic drift is accounted for by the effective population size. Heterozygosity H, defined as the probability that any two randomly selected genes will be different, is the most appropriate measure of genetic variability of a population (Nei 1987). In a haploid population of effective size N e , it follows that at equilibrium, H = 2N e µ/(2N e µ + 1) with neutral mutation rate µ (Kimura & Crow 1964), where there will be N e µ new mutants introduced per generation. If 2N e µ > 1, two or more alleles will usually be maintained. If N e is much less than µ −1 (so H is low), all the genes in the population will usually be the descendants of a single mutant.
The classical model for genetic drift is the Wright-Fisher diffusion, which is appropriate when the reproductive variance is low (Kimura 1955), i.e. the family (or sibship) sizes are all very small compared to the population size N . Genetic drift can be studied using the complementary or dual approach, Kingman's (1982) coalescent, where the probability of more than two lineages merging at each coalescent event is negligible, as N → ∞.
Star-shaped genealogies have been reported in marine populations and explained by a post-glacial expansion following a Pleistocene population bottleneck (e.g. Crandall et al. 2012). Relying on the Wright-Fisher model, these studies suggest that genetic drift has been weak, so allele frequencies have not changed rapidly in the post-expansion population.
Meanwhile, typical effective sizes of marine populations are in the hundreds or low thousands and the N e /N ratios may range from 10 −5 to 10 −2 , which are lower in larger populations (Hauser & Carvalho 2008). Motivated by considering mass spawning species such as marine fishes, there has been interest in situations where, occasionally, a single family is of appreciable size compared with the population (i.e. sweepstakes reproductive success; Hedgecock & Pudovkin 2011). Low N e /N ratios (i.e. high reproductive variances) can alter allele frequencies significantly even for very large N .
Under Wright-Fisher sampling, extreme reproduction events (whereby few individuals contribute most of the offspring to the next generation) will not occur. Although the Kingman coalescent has proven robust to violations of most of its assumptions (Kingman 2000), it drastically fails to approximate the genealogies of species with high reproductive skew (Neher & Hallatschek 2013). Recent progress has been made on describing the genealogy of populations that occasionally have very large families (see Tellier & Lemaire 2014, Grant et al. 2016, Montano 2016, where classical genetic drift is not adequate for describing random changes in allele frequencies. When extreme reproduction events occur, the temporal change in frequencies of alleles can be described by a Wright-Fisher diffusion with jumps or a generalized Fleming-Viot process, which is dual to the backward model, the Λ-or Ξ-coalescent (Donnelly & Kurtz 1999, Birkner et al. 2009, Der et al. 2011. Jumps (or rapid alternations) in the allele frequency by genetic drift forward in time are equivalent to multiple mergers in its backward-time coalescent.
This paper examines the distribution of heterozygosity across generations. The focus is to understand variation of reproductive success (i.e. variation of family sizes) within populations of species with high reproductive skew. Recognizing the importance of maintaining low-frequency alleles (Franklin 1980), it is necessary to investigate whether reproductive skew will result in alleles drifting to fixation (total consensus) within the populations, which is still poorly understood under high reproductive-variance conditions (i.e. low N e /N ratio conditions). I explore the probability to have a total consensus with zero heterozygosity.
It has been reported that the observed histograms of heterozygosity have peaks or singularities (e.g. Fuerst et al. 1977), and comparisons between the measured distributions of H and the predictions of the neutral theory were reviewed by Kimura (1983) as regards the validity of the neutral hypothesis of molecular evolution. Here I show that, although there are singularities in the H-distributions, when reproduction is highly skewed, the locations of the singular points subtly deviate from the expected norm (Higgs 1995) for the infinitely-many-neutral-alleles Wright-Fisher model.

Reproductive variance
Assume a haploid population with a large number of reproducing individuals fixed at N across generations. Now consider Pareto sampling (Huillet 2014, Huillet & Möhle 2021, Niwa 2022. Let Z 1 , . . . , Z N be independent random variables identically distributed according to the probability density with α > 0 and z ≥ 1. Upon normalizing the Z i 's by their sum one defines the weight W i of the term Z i in the sum: where R N /N 1/α , for large N , has a stable distribution called Lévy distribution (see Bouchaud & Georges 1990). Each weight W i gives the probability of reproductive success of individual i, so the i-th family recruits a fraction W i of the offspring generation. The weight is the normalized size of the family. To put it another way, given the population at some generation, for each individual at the following generation, one chooses at random with probability W i one parent i ∈ {1, . . . , N }. Such a power-law offspring-number distribution arises from type-III (exponential) survivorship with family-correlated survival (Reed & Hughes 2002, Niwa et al. 2017. Note that, when the W i 's are identical for all i, the sampling procedure is equivalent to Wright-Fisher sampling. Write ρ(w) for the distribution of weight of families in the population, such that ρ(w) dw is defined as the expected number of families formed of size between w and w + dw: which is the average of the empirical distribution for each generation (i.e. each sample or realization), where δ(·) is the Dirac delta function. The distribution ρ(w) describes among-individual variation of reproductive success that exists in the population. The probability of an individual (randomly sampled from the population) coming from a family of weight w is given by wρ(w). Since the reciprocal of the average weight of families gives the effective number of families (or reproducing lineages) in the population (Wright 1931), one has which gives the proper time-scale for the asymptotic analysis of the ancestral process (Sagitov 1999). The ratio N/N e represents the second moment (i.e. the variation) of reproductive success, which is also the effective family size in the population. If the Pareto sampling is applied to the population, then for large N , Equation (2) reduces to with (Mézard et al. 1984, Derrida & Flyvbjerg 1987, Derrida 1997, Niwa 2022, where Beta(a, b) = Γ(a)Γ(b)/Γ(a + b) is the beta function. When 1 < α < 2, although the distribution ρ(w) as in Equation (3) diverges like w −α−1 for small w, i.e. there are a large number of very small families, the effective family size is of order N 2−α . The population is then dominated by anomalously large families with size of order N 1/α (Bouchaud & Georges 1990).
In the case α < 1, Equation (3) implies that, although families are mostly concentrated around w = 0, these small-w families do not contribute to the total weight, as So any particular one of them has an extremely small weight. The probability wρ(w) is also peaked around w = 1, so there can be one large family in the interval 1/2 < w ≤ 1.

Total consensus
This section examines how heterozygosity fluctuates in a single population (or equivalently, how it varies among loci). Let n i be the number of copies of the i-th allele (i = 1, . . . , K) in the population of size N , so that the frequency of the i-th allele is x i = n i /N . Let φ(x) be the frequency spectrum, where φ(x) dx is defined as the equilibrium number of alleles expected in the frequency class (x, x + dx). When taking a random sample from the population, the probability that the allele obtained has frequency between x and x + dx is g(x) = xφ(x). Let Y be the homozygosity, so that the heterozygosity is H = 1 − Y . The probability distribution Π Y of homozygosity is defined as the average of the empirical distribution The average is taken over all possible partitions of the unit interval to K pieces Under a diffusion approximation (where µ 1 and N 1, but the scaled mutation rate θ = 2N e µ is of order 1), Kimura & Crow (1964) obtained for the infinitely-many-neutral-alleles Wright-Fisher model (in the limit K → ∞ and N → ∞ with constant θ). If Y is close to 1, the sum Y is dominated by the very high-frequency alleles, and many of the low-frequency alleles in the population do not contribute much to Y . Then, the homozygosity distribution is approximated as (4), implying that the distribution of the heterozygosity behaves like H θ−1 close to H = 0. Thus, in the Wright-Fisher model, there is a divergence as H → 0 if θ < 1. In other words, if the scaled mutation rate θ is less than the critical value of 1, no genetic variation (total consensus) in the population can be observed (see also Donnelly & Kurtz 1996, Durrett et al. 1999. To say that the probability is concentrated near H = 0 is to say that most of the time (high-frequency) alleles are near fixation. Under high reproductive-variance conditions, the value of N e is expected to be very low (expressed by the parameter θ < 1). For example, Niwa et al. (2017) obtained α = 1.24 and θ = 0.85 for Japanese sardine (Sardinops melanostictus) mitochondrial DNA cytochrome b sequences. Nevertheless, the Beta(2 − α, α) Λcoalescent predicts a much lower probability to have a total consensus with H = 0 than the Kingman coalescent, both with the same N e given by θ < 1 (Niwa et al. 2017). See also Der et al. (2011) and Montano (2016). However, in yet another class of coalescents arising from the reproductive variation (Equation 3), the Ξcoalescent, the computation of the probability of alleles drifting to fixation is much more challenging. In the following subsections, I generate probability distributions of the heterozygosity by simulation of the process of genetic drift and neutral mutation (which equivalently provides the interlocus variation of heterozygosity). At each generation each individual (gene) experiences a mutation at rate µ = θ/(2N e ), where each new mutation creates an allele that has never before existed in the population.
3.1. Wright-Fisher sampling. I performed simulations of the random sampling and mutation (with two values of θ = 0.5 and 1.0) using a population of N = 10 3 . The reproduction model employed in the simulation is either Wright-Fisher sampling or Pareto sampling with α = 2. Figure 1 shows the histograms of the heterozygosity H obtained after 10 6 generations, which gives the probability distribution of H across many realizations of the process. The numbers of hits on H = 0 for θ = 0.5 and 1.0 are, respectively, 28,693 and 885 in the case of Wright-Fisher sampling, and 1,288 and 1 in the case of Pareto(α = 2) sampling. From Ewens' (1972) sampling formula with θ = 0.5 and 1.0, the probabilities to have a total consensus in a sample of 100 sequences are 8.87 × 10 −2 and 0.01, respectively. Replicate samples (each of 100 sequences) from the Kingman coalescent with θ = 0.5 (resp. 1.0) were generated with the program ms (Hudson 2002), and 89,149 (resp. 10,085) of 10 6 replications showed total consensus. The noisy curves were calculated by explicit simulation of the sampling models, while the gray dashed curves were generated using 10 8 iterations of the recursion (Equation 6 in §4) with the probability distribution g(x) = xφ(x) given in Equation (4) (see also Derrida & Flyvbjerg 1987, Higgs 1995, Niwa 2022 The probability distributions of H have singularities (sharp changes in derivative) at all the values 1 − 1/ with = 2, 3, . . ., which become less pronounced with increasing (Derrida & Flyvbjerg 1987) and also for larger values of θ (Higgs 1995). The pattern of the distribution shown in Figure 1 (in the case of Wright-Fisher sampling or Pareto(α = 2) sampling, both with θ = 0.5) is similar to the observed histogram of heterozygosity over many loci in most species studied by Nei et al. (1976), Fuerst et al. (1977), and Singh & Rhomberg (1987). Due to the nonself-averaging property of homozygosities, i.e. the large interlocus variation in H (Stewart 1976, Higgs 1995, fairly large numbers of loci would have to be examined in order to accurately characterize the heterozygosity of a population. 3.2. Beta coalescent. I simulated 10 4 replications of 100 sequences for the infinitealleles model (θ = 0.5) under the Beta(2 − α, α) coalescent with four values of α = 1.0, 1.25, 1.5, and 1.75 through the time-reversed block counting process (Birkner & Blath 2008). Figure 2 shows that 1.06%, 0.96%, 2.37%, and 4.83% of the samples have, respectively, the total consensus sequences. When 1 < α < 2, genetic variation can be retained for a longer period of time in the populations with lower N e /N ratios (i.e. higher reproductive variances with smaller α). There are singularities in the resulting H-distributions. However, when α = 1.0 and 1.25, the positions shift from the expected value H = 1/2. 3.3. Pareto sampling. I further performed simulations of the random Pareto sampling and mutation (θ = 0.5) using a population of N = 10 3 , and four values of α = 0.9, 1.0, 1.25, and 1.5, where the numbers of hits on H = 0 are, respectively, 0, 0, 1, and 2 during the 10 6 generations of the process. A total consensus will hardly be reached among the population, and low-frequency alleles will usually be maintained.  Figure 3 shows the probability distribution of H across realizations of the process. The heterozygosity distributions are bimodal (M-shaped), and there are no divergences as H → 0. The distribution has a cutoff at the lower end (close to zero), denoted byĤ ( = 0). One then observes a peak in coalescence τ ≈Ĥ/(2µ) generations ago; very few pairs (or groups) of lineages coalesce since then (Neher & Hallatschek 2013). Within the neutral Kingman coalescent framework, a distribution of this kind would be interpreted as a rapid population expansion starting around τ generations ago. However, the size of the population did not change (with very low N e ) in the simulation. The simulation results show that, while the H-distributions have singularities, the singular points for α < 2 shift from the locations (H = 1/2, 2/3, . . .) found for Wright-Fisher sampling (or Pareto sampling with α = 2).

Singularities in the heterozygosity distribution
Now, following Derrida & Flyvbjerg (1987), let us discuss the (shifting of) singularities, assuming that the density function Π Y of homozygosity diverges at Y = 1 −Ĥ with exponent β 1 < 0, i.e. with Π (sing) denoting the dominant part of the singularity in Π Y , Let us begin with K ( 1) different alleles in the population of size N ( 1) and homozygosity Y K . The frequency of the K-th allele is x K , so that the number of copies of the K-th allele in the population is N x K . The N (1 − x K ) remaining individuals of the other K − 1 alleles in the population should have the same statistical properties. Therefore, one can write the recursion relation where Y K−1 corresponds to the homozygosity in the remaining population. So, Y K and Y K−1 have the same limiting distribution Π Y for K → ∞. The distribution Π Y is nothing but the histogram of the Y K , where the sequence Y K is constructed by the random process given by Equation (6). That is, one has a new Y value given byỸ = x 2 + (1 − x) 2 Y , where x is randomly chosen according to a given probability distribution g(x). Iterating this formula requires one random variable x for each new Y value created (cf. also Figure 1). From Equation (6) one sees the Π Y obeys the integral equation Now consider the function For Y close to 1/ −Ĥ/ 2 for = 2, 3, . . ., i.e. Y = 1/ −Ĥ/ 2 + ε with small ε > 0, one has (double-sign corresponds). Consequently the integral in Equation (7) picks up a contribution from Π Y Y 1 −1 −Ĥ ( −1) 2 . Thus, for the singularity where g(1/ ) < ∞ is assumed. Since Y = 1 −Ĥ for Y = 1/2 −Ĥ/4, the singularity of Π Y at Y = 1 −Ĥ gives rise to a singularity at Y = 1/2 −Ĥ/4. This singularity may be viewed as 'inherited' from the singularity at Y = 1 −Ĥ. In general, Y = 1/ −Ĥ/ 2 for Y = 1/( + 1) −Ĥ/( + 1) 2 with = 1, 2, . . ., yielding with β 1 as in Equation (5). An integer exponent in Equation (8) (for β 1 = −1/2) signals a logarithmic singularity. I did not prove the existence of the singularity of Equation (5) at Y = 1 −Ĥ. The simulation results imply that the distribution for allele frequencies has a cutoff at some high frequencyx < 1, i.e. for x close tox, with β 1 being the same exponent as in Equation (5). Now, consider the casex 1, thenx 2 1 −Ĥ. For Y * = 1 −Ĥ − ε (with small ε > 0), from Equation (7) one has Because Π Y is normalized, one sees that from the singularity of g(x) atx given in Equation (9). In summary, the values H = 1 − 1/ +Ĥ/ 2 (with = 1, 2, . . .) are singular points of the probability distribution of heterozygosity. The locations of the singular points depend on α (< 2) as in Equation (1), through the α dependence ofĤ (> 0). Note thatĤ → 0 as α → 2, and in the case α ≥ 2, the locations of the singular points are identical to those found for Wright-Fisher sampling.

Conclusions
The skewed distribution of reproductive success is widely observed among marine species. When reproduction is highly skewed, the value of N e is expected to be very low. The potential consequences highlight the importance of maintaining evolutionary potential. Maintaining low-frequency alleles is important to the longterm maintenance of populations under changing environmental conditions. They are not presently beneficial, although, in the future, they may be of use. Understanding genetic drift in species of conservation concern is essential for designing better conservation strategies. The application of classical population genetics theory would, however, be inappropriate for high-fecundity marine species.
This paper is concerned with the effect of reproductive skew on genetic drift. Reproduction occurs by Pareto(α) sampling. If α ≥ 2, the reproductive variance is low, so the allele frequency follows classical genetic drift, i.e. drift removes allelic variation from the population continuously at a rate inversely proportional to population size (the genetic diversity or the N e of a population is proportional to its size N ). If α < 2, classical genetic drift is not even approximately real in the broad distribution of family sizes. Jumps in the allele frequencies occur intermittently with time-intervals of order N e generations.
N e is a measure of how many individuals contribute to the next generation. In populations with low N e /N ratios (i.e. high reproductive variances), a tiny minority of adults will wind up being the parents of the vast majority of the next generation. We have seen that, even for small N e (less than the critical size expressed by θ = 1), some genetic variation may be observed under the Beta(2 − α, α) Λ-coalescent process. Further, I have provided the first analysis of this topic for Ξ-coalescent processes. The probability distribution of heterozygosity generated by Pareto sampling for reproduction (in the case α < 1 and θ < 1) has a rather sharp cutoff at the lower end, so that a total consensus with H = 0 will never be reached among the population.
Think of the scenario where few lineages substantially contribute to the next generation. This scenario expresses single lineage's rapid expansion, from which it follows that no scaled time would pass and no coalescence occur. In classical terms, rapid growth stops genetic drift, resulting in star-shaped branching patterns of (local) genealogies. However, this scenario does not imply an expansion of the population size, which can remain constant. Rather, an actual population may decline substantially without evident loss of genetic variation. Some studies have documented such a diversity paradox: overexploited and depleted marine fish populations with temporal stability in genetic diversity (see Niwa et al. 2017,and references therein). For example, Ruzzante et al. (2001) observed genetic stability for Newfoundland populations of Atlantic cod (Gadus morhua) even though the population today is only 1% of its former size. The maintenance of genetic diversity in a collapsed stock is not predicted by classical population genetic theory. Under Λ-and Ξ-coalescents, rare alleles can persist for a longer time than under Kingman coalescent, implying a reduced probability of loss (resp. fixation) for very low-frequency (resp. high-frequency) alleles. By contrast, when reproductive variance and N e are low (with θ < 1), alleles at low frequencies are more likely to be lost by drift.
In the infinite-alleles model the number K of distinct alleles actually present in the population fluctuates. The probability distribution Π Y of homozygosity has a different functional form in each interval 1/K to 1/(K + 1), so that the Π Y has singularities (Derrida & Flyvbjerg 1987, Higgs 1995. I have shown that, although there are singularities in the heterozygosity distribution, when reproduction is highly skewed (and θ < 1), the locations of the singular points subtly deviate from the expected norm for the infinitely-many-neutral-alleles Wright-Fisher model. The fact that the total consensus with H = 0 is hard to reach explains the shifting of singularities of the distribution.
The genetic diversity of a population is shaped by its recent demographic history. The shallow genealogy of the star-shaped coalescent might be a signature of rapid population growth following a bottleneck (as occurs, for example, during the last glaciation). Detecting multiple mergers in genetic data is important for understanding which forces have shaped the diversity (e.g. an excess of low-frequency alleles) of a population: post-glacial expansion vs. reproductive skew (with weak vs. strong genetic drift). The compatibility of a sample with the standard (neutral) Kingman coalescent is typically assessed using Tajima's (1989) D-statistic. One approach to identifying multiple mergers is to use the site-frequency spectrum as a summary statistic (see Neher & Hallatschek 2013, Tellier & Lemaire 2014. However, it is not trivial to determine which coalescent model is applied to account for observed genetic data (see Niwa et al. 2016). The shifting of singularities of the heterozygosity distribution is thought to be a signature of high reproductive skew in a population. When examining the interlocus variation of heterozygosity in the sample, by looking at the locations of singular points of the heterozygosity distribution, the subtle deviations (compared with the Wright-Fisher model) could help distinguish between the Kingman and multiple-merger coalescents. It would be interesting to know whether these shifting of singularities can be observed in marine species. Also it would be interesting to know what information about the process of genetic drift (reproductive stochasticity) and neutral mutation is contained in the knowledge of the shifting of singularities of the heterozygosity distribution.