How asymmetric mating patterns affect the rate of neutral genetic substitution

Introduction A population under neutral drift is expected to accumulate genetic substitutions at a fixed “molecular clock” rate over time. If the population is well-mixed, a classic result equates the rate of substitution per generation to the probability of mutation per birth. However, this substitution rate can be altered if individual birth and death rates vary by class or by spatial location. Methods Here we investigate how mating patterns affect the rate of neutral genetic substitution in a diploid, sexually reproducing population. We employ a general mathematical modeling framework that allows for arbitrary mating pattern and spatial structure. Results We demonstrate that if survival rates and mating opportunities vary systematically across individuals, the rate of neutral substitution can be either accelerated or slowed. In particular, this can occur in populations with uneven sex ratio at birth, or with reproductive skew. Discussion Our results suggest that estimates of the rate of neutral substitution, in species with uneven sex ratio and/or reproductive skew, may need to take asymmetries in mating opportunity and survival into account.


. Introduction
In a well-mixed population, for a genetic locus experiencing rare, neutral mutation, the expected rate of genetic substitution per generation, K, is equal to the probability of mutation per birth, u. This result, first demonstrated by Kimura (1968), can be derived simply as follows: the expected number of substitutions per generation is K = nuρ, where n is the number of alleles (equal to the population size N for haploids, or 2N for diploids), and ρ is the probability that a neutral mutation will become fixed. In a well-mixed population, this neutral fixation probability is ρ = 1/n, leading to K = u. This result underpins the idea of the "molecular clock"-that substitutions at neutral loci can be used to estimate the timing of events in evolutionary history (Ayala, 1997;Bromham and Penny, 2003;Kumar, 2005). The K = u result also extends to spatially structured populations with a high degree of symmetry (Maruyama, 1970;Nagylaki, 1982).
However, Allen et al. (2015) showed that asymmetries in the spatial structure can lead to either ρ > 1/n or ρ < 1/n, thereby altering the rate K of substitutions per generation. This occurs because both the fixation probability and the rate of replacement can vary across spatial locations. Sites that are frequently replaced have a higher rate of new mutations. If a frequently-replaced site is favorable to mutant fixation, this can increase the average fixation probability ρ above 1/n. In this way, asymmetric spatial structure can either speed up or slow down the neutral substitution rate K.
The results of Allen et al. (2015) were derived for haploid, spatially structured populations. Here we show that the same phenomenon can arise for diploid populations with asymmetric mating patterns. We present a general method for calculating the neutral fixation probability ρ under arbitrary mating pattern. This fixation probability can be compared to the corresponding value of ρ WM = 1 2N for a well-mixed population of N individuals. We apply this framework to two hypothetical scenarios: a randomly mating population with uneven sex ratio, and a population with reproductive skew in one sex. In each case we show that asymmetries in the mating pattern and replacement rate can lead to ρ > 1 2N or ρ < 1 2N , accelerating or slowing the rate K of neutral genetic substitution. Specifically, we identify two circumstances in which the neutral substitution rate is altered: (1) the sex ratio is unequal at birth and there is differential mortality between the sexes, or (2) there is reproductive skew in one sex, and differential mortality between high-reproducing and low-reproducing individuals. The first scenario may arise in reptiles with temperature-dependent sex determination (Ferguson and Joanen, 1982;Ewert and Nelson, 1991), while the second may apply to a variety of birds (Wiley, 1973;Petrie et al., 1991;Höglund and Alatalo, 1995;Magrath et al., 2004), mammals (Boesch et al., 2006;Raihani and Clutton-Brock, 2010;Sherman et al., 2017;Higham et al., 2021), and insects (Adams and Atkinson, 2008;Shimoji and Dobata, 2022).

. Modeling framework
We begin by summarizing the mathematical framework from which our results are derived. A formal description of this framework is given in the Supplementary material.

. . Mating and replacement
We consider a diploid population of fixed size N, which can be monoecious (one sex) or dioecious (two sexes). The individuals are indexed I = 1, . . . , N. Depending on the model in question, each index I may also designate other relevant information, such as the sex and/or spatial location of individual I.
Each time-step, some subset of individuals are replaced by new offspring. Each new offspring has two parents (which may be the same individual in the case of self-mating). The set of individuals to be replaced, as well as the parents of each new offspring, are sampled from a fixed probability distribution. This probability distribution can be chosen arbitrarily, subject only to the requirement that it be possible for a mutant allele to spread throughout the population. In this way, our framework can incorporate a wide variety of spatial structures and mating patterns, with no assumption of symmetry or uniformity. This level of generality is achieved using an abstract mathematical representation of population structure, which was developed in previous works (Allen and Tarnita, 2014;Allen et al., 2015;Allen and McAvoy, 2019) and is further extended in the Supplementary material to include arbitrary mating patterns.
The key quantities describing any particular model within this framework are the mating and replacement rates E K IJ , where the FIGURE Neutral drift in a diploid population. In each time-step, some individuals are replaced by the o spring of others. These mating and replacement events are sampled from a fixed probability distribution, which captures all e ects of spatial structure and mating pattern.
The key quantities are the marginal probabilities E K IJ that, in a given time-step, individual K is replaced by the o spring of I and J. There are two allele types, mutant and resident. Each new o spring inherits alleles from its parents according to Mendelian inheritance.
indices I, J, and K run from 1 to N, representing individuals in the population (Figure 1). For I = J, E K IJ is equal to the marginal probability that, in a given time-step, individual K is replaced by the offspring of I and J. The subscript indices in E K IJ are understood as an unordered pair, meaning that E K IJ and E K JI both refer to the same marginal probability. In the case of self-mating, E K II is set equal to twice the marginal probability that K is replaced by the self-mated offspring of I (the factor of two accounts for the doubled chance to spread one's alleles under self-mating). These E K IJ quantities take into account the population's spatial structure and dispersal patterns as well as its mating pattern.
For each individual I = 1, . . . , N, we define the death rate D I as the marginal probability that I is replaced, and the birth rate B I as half the expected offspring number of I. In terms of the mating and replacement rates E K IJ , the birth and death rates are given by The factors of 1 2 in B I and D I arise for different reasons: for D I a 1 2 is needed to avoid double-counting in the J and K indices, while for B I the 1 2 represents the chance of transmitting a particular allele under Mendelian inheritance. We also define the overall rate of turnover, .

. Alleles and states
Within the population we consider two allele variants: a resident allele R, and a neutral mutant allele M. Each new offspring inherit alleles from its parents according to Mendelian inheritance. We assume there is no recurring mutation, so that, over time, either the mutant or the resident allele will become fixed in the population.
The state of the population at any given time can be specified by identifying the genotype (MM, MR, or RR) of each individual I = 1, . . . , N. Since the mutant allele is neutral, the mating and replacement rates E K IJ do not depend on the current population state. The overall process of neutral drift is represented as a Markov chain on the set of population states (see the Supplementary material for details).

. Fixation probability
Over time, the neutral drift process will ultimately converge to a state where either only mutant or only resident alleles are present. We are interested in the fixation probability of the mutant type; that is, the probability of reaching the all-mutant state when starting from a state with only a single mutant allele.

. . Deriving fixation probability
Let ρ I denote the probability that the all-mutant state is reached, when starting from a state with a single mutant allele copy in individual I. The neutral fixation probability ρ I can also be understood as the reproductive value of individual I (Taylor, 1990;Lehmann, 2014;Maciejewski, 2014;Allen and McAvoy, 2019), since it quantifies I's contribution to the future gene pool.
In the Supplementary material we derive the recurrence equation for the reproductive values ρ I : The two terms on the right-hand side correspond to survival and transmission, respectively, of a mutant allele. The 1 2 in the second term reflects the chance of transmission to offspring under Mendelian inheritance. These probabilities of survival or transmission are multiplied by the reproductive values, ρ I or ρ K , of this allele and its copies in the next time-step.
We can write Eq. (4) in simpler form as In order to solve for the fixation probabilities, we also require the equation Equation (6) represents the fact that exactly one of the 2N initial alleles will ultimately spread its descendants throughout the population, and the probability of this occurring for an allele in individual I is ρ I . Together, Eqs. (5) and (6) form a system of equations that uniquely determine ρ I for each I = 1, . . . , N. Although this system involves N + 1 equations, one case of Eq. (5) is always redundant and can be eliminated. If the population can be subdivided into classes (by sex, spatial location, etc.) such that individuals in a given class are interchangeable, only one equation per class is needed. This will be the case in the examples we explore later.
To determine the overall fixation probability, we must take into account where mutations arise. Supposing that mutations occur with constant probability per birth, the likelihood of a new mutant allele arising at individual I is proportional to the turnover rate D I . We therefore suppose the initial mutation has probability D I /B of arising within individual I, for each I = 1, . . . , N. This leads to an overall fixation probability of . . Comparing to well-mixed population A well-mixed, monoecious (one sex) population with uniform random mating has the same value of E K IJ for each triple I, J, K. In this case, Eqs. (1)-(7) yield an overall fixation probability of ρ WM = 1 2N , which we take as our baseline value. We are interested in when and how ρ can differ from ρ WM , in populations for which the replacement rates E K IJ take on non-uniform values due to asymmetric mating patterns and/or spatial structure.
We first note that, by combining Eqs.
(3), (6), and (7), we can write Above, Cov I [D, ρ] is the population covariance of D I with ρ I over I = 1, . . . , N, which is equal by definition to the quantity inside parentheses in the second expression. It follows that ρ > ρ WM if D I and ρ I are positively correlated, ρ < ρ WM if D I and ρ I are negatively correlated, and ρ = ρ WM if they are uncorrelated. In particular, we will have ρ = ρ WM if either D I or ρ I are constant over individuals. These two cases correspond to Results 1 and 2 of Allen et al. (2015), which we extend here to diploid populations. First, Eq. (8)  That is, if mortality is constant across individuals-mathematically, if I,J E K IJ is uniform over K-the rate of neutral substitution is unchanged from the baseline rate.
The second case, constant ρ I , turns out to be equivalent to each individual having birth rate equal to death rate: Result 2. Reproductive values are constant over individuals (i.e., ρ I = 1 2N for all I) if and only if B I = D I for all I.
The condition B I = D I can also be written as J,K E K IJ = J,K E I JK for each I. Intuitively, if each individual has birth rate equal to death rate, then any advantage in reproduction is canceled by a disadvantage in survival. Consequently, mutations are equally likely to become fixed no matter where they arise, and the neutral substitution rate is again unchanged.
Proof of Result 2. Suppose first that ρ I = 1 2N for all I. Substituting in Eq. (5) and invoking Eq. (1), we have and hence B I = D I . Conversely, suppose that D I = B I for all I. By Eq.
(1), this means By inspection, ρ I = 1 2N for all I is a solution to Eq. (9) as well as to Eq. (6). Since the solution to this system is unique, we have ρ I = 1 2N for all I.
Finally, Result 3 of Allen et al. (2015) shows that when birth rates are uniform, the fixation probability can never exceed its well-mixed value. Extended to diploid populations, this becomes: Result 3. If B I is constant over all I = 1, . . . , N then ρ ≤ 1 2N , with equality if and only if the D I are also constant over I.
Constant B I is equivalent to J,K E K IJ being uniform over I. Result 3 implies that uniform birth rates impose a "speed limit" on the neutral substitution rate, and this speed limit is achieved only if death rates are also uniform. In contrast to Results 1 and 2, the proof of Result 3 is rather involved. In the Supplementary material, we reduce this result to the haploid case, which was proven by Allen et al. (2015).

. Scenario : Uneven sex ratio
We now proceed to apply the above framework to example scenarios, showing how ρ can deviate from ρ WM under asymmetric mating and replacement patterns. We first consider a scenario of uneven sex ratio. Uneven sex ratios may arise at birth, as in crocodilians (Ferguson and Joanen, 1982;Lang and Andrews, 1994;González et al., 2019) and turtles (Ewert and Nelson, 1991;Mrosovsky, 1994), or due to differential survival, as is the case in many bird species (Benito and González-Solís, 2007;Donald, 2007).
To model a population with uneven sex ratio, we suppose the set of individuals is partitioned into two subsets, labeled M for males and F for females ( Figure 2A). The numbers of males and females are denoted m and f , respectively. The total population size is N = m + f .
Each time-step, one male and one female are chosen, uniformly at random within each sex, to mate and produce a single offspring. Each offspring replaces a single individual, with each male having probability D M to be replaced, and each female having probability D F . Since exactly one individual is replaced, D M and D F are constrained by the relationship We treat D M and D F as tunable parameters in this model, subject to the above constraint. We observe that D M and D F are equal to the marginal probability that an individual (male or female, respectively) is replaced in one time-step, in accordance with the definition of D I in Section 2.1. The mating and replacement probabilities E K IJ are then obtained by multiplying the probability that a specific male-female pair is chosen ( 1 mf ) by the probability that the offspring replaces a specific male or female individual (D M or D F , respectively). This leads to two distinct values, depending on the sex of the replaced individual: The primes in M ′ and F ′ distinguish the replaced individual from the parent of the same sex. Since a single individual is replaced each time-step, the overall replacement rate is B = 1. For the birth rate of an individual from each sex, we have reflecting the fact that exactly one individual from each sex is chosen to mate. Since each new offspring replaces a deceased individual, the offspring born in each time-step is male with probability mD M and female with probability fD F . The sex ratio at birth is therefore (mD M )/(fD F ). This may differ from the adult sex ratio of m/f due to differential mortality.

. . Deriving fixation probability
Since members of each sex are interchangeable, there are only two cases of the recurrence relation for reproductive values, Eq. (5): Upon substituting from Eq. (11), both equations in Eq. (13) simplify to . /fevo. . This means that the reproductive values of the sexes are in inverse ratio to their frequencies at birth: Solving Eqs. (14) and (15) yields: Substituting into Eq. (7) gives the overall fixation probability: This solution is illustrated in Figure 3.

. . Comparing to well-mixed population
Since the fixation probability in a well-mixed population is Applying Eq. (10), this condition simplifies to We conclude that the fixation probability is unchanged from the baseline value, ρ = ρ WM , if either D M = D F (equal turnover per sex) or mD M = fD F (equal sex ratio at birth). The first case is an instance of Result 1, i.e., ρ = ρ WM if each individual is replaced at the same rate. In the second case, Eqs. (10) and (12) imply This is an instance of Result 2, i.e., ρ = ρ WM if each individual has birth rate equal to death rate. In this case, since the sex ratio is equal at birth (mD M = fD F ), any bias in the adult sex ratio is solely due to differential survival. The minority sex gains an advantage in reproduction, but this is exactly canceled by their higher mortality, leading to equal reproductive values, ρ M = ρ F . According to Eq. (19), ρ > ρ WM (meaning the rate of neutral substitution is accelerated) if and only if D M − D F and mD M − fD F have opposite signs, meaning that the minority sex at birth has higher mortality. Conversely, if the minority sex at birth has lower mortality, then ρ < ρ WM and the rate of neutral substitution is slowed. This can be seen in Figure 3. We observe that if m = f , the two solutions to ρ = ρ WM coincide, and ρ ≤ ρ WM with equality if and only if D M = D F . This is an instance of Result 3, since B M = B F when m = f .

. . Maximizing fixation probability
To determine the extent to which the fixation probability ρ can deviate from its well-mixed value, we apply Lagrange multipliers to Eq. (16), with D M and D F as variables and Eq. (10) as a constraint. This leads to a critical point of Frontiers in Ecology and Evolution frontiersin.org . /fevo. . Substituting into Eq. (17), we obtain the maximal fixation probability: This maximal value corresponds to the peaks of the curves in Figure 3. The ratio of ρ MAX to ρ WM can be written in terms of the sex ratio m/f : As the sex ratio approaches either zero or infinity (the two possible extremes), the ratio ρ MAX /ρ WM converges to 2. Thus, under extreme sex ratios, the rate of neutral substitution may be accelerated (relative to the well-mixed population) by a factor of up to two.
In presenting the model we will refer to reproductive skew in males, but the labels "male" and "female" are arbitrary, and the results apply equally to skew in either sex.
As in the earlier scenario, we consider a population consisting of a set M of males and a set F of females. In this case, however, the males are partitioned into two subsets, M 1 and M 2 , where M 1 males have a greater chance of reproducing. Specifically, the reproductive rate of M 1 males exceeds that of M 2 males by a factor r > 1. The numbers of M 1 males, M 2 males, and females are denoted m 1 , m 2 , and f respectively. The total population size is N = m 1 + m 2 + f . Each time-step, a single offspring is produced by a single malefemale pair. Each particular M 1 male has probability r/(rm 1 + m 2 ) to reproduce, while each particular M 2 male has probability 1/(rm 1 + m 2 ). Each female has probability 1/f to reproduce. This new offspring replaces a deceased adult; each M 1 male, M 2 male, and female has probability D M 1 , D M 2 , and D F , respectively, to be replaced in a given time-step. Since exactly one individual is replaced, these death probabilities are constrained by Since the new offspring replaces the deceased adult, the new offspring is a M 1 male, M 2 male, or female with respective probabilities m 1 D M 1 , m 2 D M 2 , or fD F . In particular, the sex ratio at birth is (m 1 D M 1 + m 2 D M 2 )/(fD F ). From this model description, we derive the following probabilities for a specific male-female pair to replace a specific individual: Frontiers in Ecology and Evolution frontiersin.org . /fevo. .
Each probability E K IJ is the product of three probabilities: that of choosing a male for reproduction, choosing a female for reproduction, and choosing an individual for replacement.

. . Deriving fixation probability
For this model, the recurrence equation for fixation probabilities, Eq. (5), becomes Substituting the marginal probabilities from Eq. (25), Eqs. (26a) and (26b) simplify to Observing that Eqs. (27a) and (27b) differ only in a factor of r on the right hand side, we conclude that This reduces Eq. (27a) to Meanwhile, Eq. (6) gives Solving the system of equations formed by Eqs. (28)-(30), we obtain the solution The overall fixation probability is therefore .

. Comparing to well-mixed population
Comparing our solution for ρ in Eq. (32) to the fixation probability in a well-mixed population of the same size, ρ WM = 1/(2N) = 1/(2m 1 + 2m 2 + 2f ), we find that ρ > ρ WM is equivalent to To simplify, let us suppose that the sex ratio is even both for adults, f = m 1 + m 2 , and at birth, m 1 D M 1 + m 2 D M 2 = fD F = 1 2 . Then Condition (33) reduces to Simplifying using the relation m 1 D M 1 + m 2 D M 2 = 1 2 , we obtain that ρ > ρ WM if and only if There are two cases where ρ = ρ WM , corresponding to the two factors of the left-hand side of Condition (35). Setting the first factor equal to zero leads to This is an example of Result 1 from Section 3.2. Setting the second factor of the left-hand side of Condition (35) equal to zero leads to Frontiers in Ecology and Evolution frontiersin.org . /fevo. . This is an example of Result 2, in that each individual has birth rate equal to its replacement rate. According to Condition (35), ρ > ρ WM (meaning the neutral substitution rate is accelerated) when D M 1 lies between D M 2 and rD M 2 . This arises when M 1 males have higher mortality than M 2 males, but not so high as to negate their reproductive advantage. If D M 1 is smaller than D M 2 or larger than rD M 2 , then ρ < ρ WM and the neutral substitution rate is slowed. This can be seen in Figure 4, where the solution in Eq. (36) is the common intersection point for all curves, while the solution in Eq. (37) varies depending on the parameters.

. . Maximizing fixation probability
To determine how widely the fixation probability can vary for fixed m 1 , m 2 , f , and r, we seek a maximum of the function subject to the constraint of Eq. (24). Applying Lagrange multipliers gives the equations Dividing (39a) by (39b) and simplifying gives while dividing (39a) by (39c) yields Invoking (24) and solving, we obtain the maximizing values Substituting these expressions into (38) gives the maximal fixation probability of Comparing to the well-mixed population, we find the ratio of fixation probabilities is given by To determine the most extreme values of fixation probability, we take r → ∞ in Eq. (43), yielding This is equivalent to Eq. (22) with m replaced by m 1 , since for r → ∞ only M 1 males have the opportunity to reproduce. In this case, an arbitrarily large ratio ρ MAX /ρ WM can be achieved by considering a large number of M 2 males: .

. Maximization under even sex ratio
The preceding section shows that, in our reproductive skew model, the neutral substitution rate may be scaled by any positive factor relative to the well-mixed rate. However, this is achieved only with an extreme male-biased sex ratio, in that m 2 is taken to infinity with m 1 and f held constant. It is also of interest to determine how widely the neutral substitution rate can vary when the overall sex ratio is held even, both at birth and among adults. This is the case depicted in Figure 4.
To this end, we seek to maximize ρ under the constraints that the sex ratio is even among adults, m 1 + m 2 = f , and at birth, m 1 D M 1 + m 2 D M 2 = fD F = 1 2 . Together, these constraints imply D F = 1/(2(m 1 + m 2 )). Substituting this into Eq. (32) and simplifying, we obtain To maximize ρ, we apply Lagrange multipliers in the variables D M 1 and D M 2 with the constraint m 1 D M 1 + m 2 D M 2 = 1/2. We obtain that ρ is maximized when D M 1 = √ rD M 2 ; specifically, Substituting in Eq. (47) gives the maximal fixation probability under these constraints: Eqs. (48)-(49) correspond to the maxima of the curves depicted in Figure 4. Comparing to well-mixed value of ρ WM = 1/(4(m 1 + m 2 )), we have Frontiers in Ecology and Evolution frontiersin.org . /fevo. . In the extreme scenario r → ∞, in which M 1 males completely dominate reproduction, we obtain If, in addition, we take the ratio m 1 /m 2 to zero (meaning M 1 males are a very small minority), then ρ MAX /ρ WM converges to 2, which is an upper bound for ρ MAX /ρ WM in the case of even sex ratios at birth and in adults.

. Discussion
For the fixation probability of a neutral mutation to differ from the well-mixed of ρ WM = 1 2N , three conditions must be met. First, the rates of replacement D I must vary across individuals (or classes of individuals), or else ρ = ρ WM by Result 1. Second, the reproductive values ρ I must vary across individuals; by Result 2 this occurs if and only if at least one individual (or class) has birth rate not equal to its death rate. Third, the variation in D I and ρ I must be correlated, so that 1 N I D I ρ I differs from 1 N I D I 1 N I ρ I . If these three conditions are met, Kimura's (1968) result K = uwhich underlies the molecular clock hypothesis (Kimura, 1983;Ayala, 1997;Bromham and Penny, 2003;Kumar, 2005)-requires a correction factor of ρ/ρ WM . For a diploid population of size N, the corrected rate of neutral substitution is K = 2Nρu.
A key assumption underlying our results is that mutation occurs with constant probability per birth. This fact, combined with unequal replacement rates, means that mutations are more likely to arise in certain subgroups rather than others; if these subgroups also have larger-than-average reproductive value, then the rate of neutral substitution will be accelerated.
It has previously been shown that the rate of genetic substitution can be affected by selection, varying population size (Balloux and Lehmann, 2012), mutation rates that vary by age or sex (Pollak, 1982;Charlesworth, 1994;Lehmann, 2014), or asymmetric spatial structure (Allen et al., 2015). To this list we now add uneven sex ratio, reproductive skew, and other asymmetric mating patterns, in combination with unequal replacement rates. For uneven sex ratio, the correction factor ρ/ρ WM can take any value strictly between 0 to 2 ( Figure 3). Deviations from the well-mixed rate occur when there is both differential survival, D M = D F , and unequal sex ratio at birth, mD M = fD F . If the minority sex at birth has higher (resp., lower) mortality, the neutral substitution rate is accelerated (resp., slowed). According to Eq. (21), the maximal substitution rate occurs when the ratio of replacement rates is inversely proportional to the square root of the sex ratio, D M /D F = f /m-or equivalently, when the adult sex ratio is the square of the sex ratio at birth, m/f = (mD M ) 2 /(fD F ) 2 .
In order to affect the neutral substitution rate in our model, the uneven sex ratio must arise at birth. If, instead, the uneven sex ratio is solely due to differential survival, then the reproductive advantage of the minority sex is canceled by their higher mortality. This limits the range of taxa for which such effects may arise, since unequal sex ratios at birth violate random Mendelian segregation as well as Fisher's (Fisher, 1930) principle that sex ratios should evolve toward unity. Although most bird species have even sex ratio at birth (Donald, 2007), exceptions have been observed, for example, in some tit species (Cichoń et al., 2005) and in collared flycatchers (Ficedula albicollis; Rosivall et al., 2004;Cichoń et al., 2005). A richer source of possible examples are reptiles with temperature-dependent sex determination (Ferguson and Joanen, 1982;Ewert and Nelson, 1991;Lang and Andrews, 1994;Mrosovsky, 1994;González et al., 2019), which often have uneven sex ratio under baseline environmental conditions (Ferguson and Joanen, 1982;Ewert and Nelson, 1991).
For reproductive skew, the correction factor ρ/ρ WM can take any positive value. If we additionally specify that the sex ratio is even at birth and among adults, then ρ/ρ WM varies strictly between 0 and 2. Acceleration of the neutral substitution rate occurs when D M 1 lies between D M 2 and rD M 2 , meaning that the dominant reproducers have higher mortality than others of their sex, but not so high as to negate their reproductive advantage.
Here we have considered effects due to sex and space, but not age. The age structure of a population has clear implications for its neutral substitution rate (Pollak, 1982;Charlesworth, 1994;Lehmann, 2014). Investigating the combined effects of age, sex, and spatial structure is a promising avenue for future work. Incorporating age structure would also allow for greater realism in applying this framework to real-world populations (Charlesworth, 1994;Sample et al., 2018). Age structure could be incorporated into the framework described in the Supplementary material by designating a fixed subset of individuals to be juveniles. The mating and replacement probabilities E K IJ would then be nonzero only if I and J are adults and K is a juvenile. It would also be necessary to specify the probability for each juvenile to survive to adulthood. This would represent a significant-and important-expansion of the mathematical modeling framework developed by Allen and McAvoy (2019).
Our analysis has focused on the rate of neutral substitution at a single genetic locus. It would also be of interest to extend to multiple loci. This would require amending the framework presented in the Supplementary material to keep account of the two alleles at each locus in each (diploid) individual. The probability of recombination in each new offspring would then be a key parameter. Such an expanded framework could bring new methods to bear on the well-studied question of how mating patterns affect linkage disequilibrium (Weir and Cockerham, 1969;Golding and Strobeck, 1980;Nordborg, 2000).
We have also implicitly assumed that the rate of neutral substitution is primarily limited by the appearance rate and fixation probability of neutral mutations. This assumption is reasonable so long as the expected fixation time, T, is significantly less than the waiting time, 1/(2Nuρ), for successful mutations to arise. Outside of this regime, fixation time becomes relevant for the neutral substitution rate (Frean et al., 2013). Fixation times have been studied extensively for spatially structured populations (Frean et al., 2013;Hathcock and Strogatz, 2019;Möller et al., 2019;Tkadlec et al., 2019), and to a lesser extent for populations with asymmetric mating pattern (Paley et al., 2010). Examining the combined effect of fixation probability and time on the rate of neutral substitution, in populations with asymmetric mating pattern, is an interesting avenue for further study.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.