# A new and almost perfectly accurate approximation of the eigenvalue effective population size of a dioecious population: comparisons with other estimates and detailed proofs

Thierry de Meeûs, C. Noûs

► To cite this version:

Thierry de Meeûs, C. Noûs. A new and almost perfectly accurate approximation of the eigenvalue effective population size of a dioecious population: comparisons with other estimates and detailed proofs. Peer Community Journal, 2023, 3, e51 [50 p.]. 10.24072/pcjournal.280. hal-04148819

HAL Id: hal-04148819

https://hal.science/hal-04148819v1

Submitted on 11 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Published
2023-05-31

**Cite as**
Thierry De Meeûs and Camille Noûs (2023) *A new and almost perfectly accurate approximation of the eigenvalue effective population size of a dioecious population: comparisons with other estimates and detailed proofs*, Peer Community Journal, 3: e51.

**Correspondence**
thierry.demeeus@ird.fr

**Peer-review**
Peer reviewed and
recommended by
PCI Evolutionary Biology,
https://doi.org/10.24072/pci.evolbiol.100651

(cc) BY

This article is licensed
under the Creative Commons
Attribution 4.0 License.

# A new and almost perfectly accurate approximation of the eigenvalue effective population size of a dioecious population: comparisons with other estimates and detailed proofs

Thierry De Meeûs<fnref n="1" /> and Camille Noûs²

Volume 3 (2023), article e51

https://doi.org/10.24072/pcjournal.280

## Abstract

The effective population size is an important concept in population genetics. It corresponds to a measure of the speed at which genetic drift affects a given population. Moreover, this is most of the time the only kind of population size that empirical population genetics can give access to. Dioecious populations are expected to display excesses of heterozygosity as compared to monoecious panmictic populations, as measured by Wright's *F*<sub>IS</sub>. It can be shown that these excesses are negatively correlated with the population size. This is why *F*<sub>IS</sub> can be used to estimate the eigenvalue effective population size of dioecious populations. In this paper, we propose a new approximation that provides a very accurate estimate of the eigenvalue effective population size of a dioecious population as a function of the real population size. We then explore the accuracy of different *F*<sub>IS</sub>-based methods using the leading eigenvalue of transition matrices or coalescence. It appears that the eigenvalue-based method provides more accurate results in very small populations, probably due to approximations made by the coalescence approach that are less valid in such situations. We also discuss the applicability of this method in the field.

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Univ Montpellier, Cirad, IRD, Intertryp, Montpellier, France, ²Cogitamus laboratory (https://www.cogitamus.fr/), France</footnote>
## Introduction

A convenient way to measure the speed at which a given population loses its genetic diversity by genetic drift is to compute its effective population size $N_e$ (Vitalis & Couvet, 2001b). Several formal definitions exist. They all refer to an ideal population that follows all Castle-Weinberg assumptions (Castle, 1903; Weinberg, 1908) (see De Meeûs et al. (2021) for an explanation why this labelling is fairer than the more popular Hardy-Weinberg), except for the size of the population that is limited to $N_e$. It means a self-compatible diploid monoecious and panmictic population of size $N_e$, with no selection, no migration, no mutation and discrete generations, and where alleles for the next generation are binomially sampled from the $2N_e$ available ones, which causes genetic drift at a specific speed, proportional to $1/(2N_e)$. Such a population is also known as following the Wright-Fisher (WF) model (Crow & Kimura, 1970), as opposed to the Castle-Weinberg model, where the population is of infinite size, and thus without genetic drift. Some approaches focus on the rate of inbreeding increase, the rate of heterozygosity loss, the variation of allele frequencies from one generation to the other (Vitalis & Couvet, 2001b), or the coalescence time (Balloux & Lehmann, 2003; Balloux et al., 2003; Balloux, 2004; Nomura, 2008). This led authors to define the inbreeding effective population size, which refers to the speed at which inbreeding evolves, the eigenvalue effective population size (see appendices 1-3 to see the detailed analytical tools and Appendix 4 to see why it was named as such), the variance (of allele frequencies from one generation to the next) effective population size (Crow & Kimura, 1970; Vitalis & Couvet, 2001b; Ewens, 2004) and the coalescence (or coancestry) effective population size (Balloux & Lehmann, 2003; Balloux et al., 2003; Balloux, 2004; Nomura, 2008) (see below), respectively. In all cases, the effective population size is computed for a given population of census size $N$, which deviates from an ideal population (following WF) at one or several of the properties defined above. Because of these deviations, genetic drift operates at a faster rate, or sometimes at a slower rate, than the same population if it fulfilled the ideal conditions. The effective population size of such a non-ideal population is the ideal population of size $N_e$ that would drift at the same speed as the non-ideal one, also known as the size of a population following WF and drifting at the same speed as the focal population (Vitalis & Couvet, 2001b).

Many species have separate sexes. Several authors have investigated the impact that dioecy and sex ratio have on effective population size. In this note, we review some of these results and we then derive a new and apparently more accurate approximation for the eigenvalue effective population size of a dioecious population. We also propose another estimator of $N_e$ from Wright's $F_{IS}$. We compare the relative performances of the different methods. Several appendices present the proofs of all equations used in the main text. These appendices are extremely detailed so that almost anybody willing to understand precisely how a given result was obtained, here and in the cited literature, can easily access this knowledge. Nevertheless, more skilled readers will probably not need to read any of those.

## Classic results from the literature

The effective population size of a dioecious population has been defined in different ways. In Wright's book (Wright, 1969) (page 197), the approximate (eigenvalue) effective size of a population of size $N$, with $N_f$ females and $N_m$ males is:

$$ (1) \qquad N_e = \frac{4N_f N_m}{N} $$

Nevertheless, in the same book (page 197 again), a better approximation is suggested, and a quick proof can be found in Felsenstein's book (pages 266-267) (Felsenstein, 2019) (see also below, equation 15), for the eigenvalue effective population size $N_e$:

$$ (2) \qquad N_e = \frac{4N_f N_m}{N} + \frac{1}{2} $$

More recently, Balloux (Balloux, 2004), computed the coalescence effective population size as:
$$ (3) \qquad N_e = \frac{4N_f N_m}{N} + \frac{1}{2} + \frac{1}{2N} $$

From equations (2) and (3), and for sex-ratios (SR) that are not too female biased (e.g. $N_f - N_m < \sqrt{N/2}$, for Equation 2), one can see that dioecy tends to slightly increase $N_e$. This is obviously a consequence of the supplementary delay required for two alleles to become identical by descent in the same individual, since selfing cannot occur. Another consequence is that dioecious populations are expected to display heterozygote excesses (Robertson, 1965). This led to formalizing the expected deviation of heterozygote frequency in dioecious populations, as measured by Wright's $F_{IS}$ (Wright, 1965), which may provide a simple tool to estimate the effective population size, assuming even sex ratios. Using simple algebra on observed and expected heterozygosity, Pudovkin et al. (1996) proposed the following equation (see Appendix 5 for a detailed proof) for the eigenvalue effective population size:

$$ (4) \qquad N_e = -\frac{1}{2F_{IS}} $$

They also proposed a supposedly more accurate approximation with their equation 4 (but see also Appendix 5):

$$ (5) \qquad N_e = -\frac{1}{2F_{IS}} + \frac{1}{2(1-F_{IS})} $$

Balloux (2004) proposed another solution, based on the coalescent effective population size and requiring quite cumbersome analytic treatments, which are detailed in Appendix 6:

$$ (6) \qquad N_e = -\frac{1}{2F_{IS}} - \frac{F_{IS}}{1+F_{IS}} $$

## The general model of a dioecious pangamic population

For now, and unless specified otherwise, we assume a dioecious diploid population of constant size and sex-ratio, with discrete generations, no mutation and no migration. At each generation, alleles that will be present in an individual of generation $t$ were randomly drawn with replacement in the pool of gametes of their two parents (or from infinite pools of gametes), and females and males are polygamous and mate randomly. For a dioecious population with $N_f$ females and $N_m$ males, probabilities of identity within individuals $Q_{I(t)}$ and between individuals within the subpopulation $Q_{S(t)}$ at generation $t$ respectively are (see for instance Balloux (2004), equation 14 or Felsenstein (2019), pages 266-267):

$$ (7) \qquad Q_{I(t)} = Q_{S(t-1)} $$

and

$$ (8) \qquad Q_{S(t)} = \frac{1}{4} \left[ \frac{1}{N_f} \left( \frac{1}{2} + \frac{1}{2} Q_{I(t-1)} \right) + \left( 1 - \frac{1}{N_f} \right) Q_{S(t-1)} \right] + \frac{1}{4} \left[ \frac{1}{N_m} \left( \frac{1}{2} + \frac{1}{2} Q_{I(t-1)} \right) + \left( 1 - \frac{1}{N_m} \right) Q_{S(t-1)} \right] \\ + \frac{1}{2} Q_{S(t-1)} $$

Equation 7 is straightforward: inbreeding of individuals at any generation comes from the genetic similarity between their parents. Equation 8 is less intuitive. We want to compute the probabilities of identity by descent between two alleles of two different individuals taken at random at generation $t$, and assuming random mating of parents and a great number of matings $n_M$ (i.e. $n_M \to \infty$). This way, the probability of mating between two individuals remains independent of previous copulas these may have been involved in. The probability that two alleles of generation $t$ come from two females of generation $t-1$ is $1/4$, from two males is also $1/4$, and from one male and one female is $1/2$. If both come from two females, the probability they came from the same mother is $1/N_f$ (i.e. found in full or half sibs) or two different females is $1-1/N_f$. In the same diploid individual, the probability to sample twice the same allele is $1/2$. Otherwise,
two different alleles are sampled from the same individual with probability 1/2, but in that case they are identical by descent with probability $Q_{I(t-1)}$. If the two alleles came from two different females, the probability that these two alleles are identical by descent is $Q_{S(t-1)}$. The same reasoning applies to the two-males case. For the one-female-one-male case, the probability that the two alleles are identical is also $Q_{S(t-1)}$.

Combining equations 7 and 8, we get:

$$ 
\begin{aligned}
Q_{S(t)} &= \frac{1}{4} \left[ \frac{1}{N_f} \left( \frac{1}{2} + \frac{1}{2} Q_{S(t-2)} \right) + \left( 1 - \frac{1}{N_f} \right) Q_{S(t-1)} \right] \\
&\quad + \frac{1}{4} \left[ \frac{1}{N_m} \left( \frac{1}{2} + \frac{1}{2} Q_{S(t-2)} \right) + \left( 1 - \frac{1}{N_m} \right) Q_{S(t-1)} \right] + \frac{1}{2} Q_{S(t-1)} \\
Q_{S(t)} &= Q_{S(t-1)} \frac{1}{4} \left( 1 - \frac{1}{N_f} + 1 - \frac{1}{N_m} + 2 \right) + Q_{S(t-2)} \frac{1}{4} \left( \frac{1}{2N_f} + \frac{1}{2N_m} \right) + \frac{1}{4} \left( \frac{1}{2N_f} + \frac{1}{2N_m} \right)
\end{aligned}
$$

<figure><img src="image_3.png" /></figure>

$$ Q_{S(t)} = Q_{S(t-1)} \frac{1}{4} \left( \frac{4N_f N_m - N_m - N_f}{N_f N_m} \right) + Q_{S(t-2)} \frac{1}{8} \left( \frac{N_f + N_m}{N_f N_m} \right) + \frac{1}{8} \left( \frac{N_f + N_m}{N_f N_m} \right) $$

<figure><img src="image_4.png" /></figure>

$$ Q_{S(t)} = Q_{S(t-1)} \frac{1}{4} \left( \frac{4N_f N_m - N}{N_f N_m} \right) + Q_{S(t-2)} \frac{1}{8} \left( \frac{N}{N_f N_m} \right) + \frac{1}{8} \left( \frac{N}{N_f N_m} \right) $$

Using the genetic diversity of the subpopulation $H_s=1-Q_s$ yields simpler expressions:

$$ 1 - H_{S(t)} = (1 - H_{S(t-1)}) \frac{1}{4} \left( \frac{4N_f N_m - N}{N_f N_m} \right) + (1 - H_{S(t-2)}) \frac{1}{8} \left( \frac{N}{N_f N_m} \right) + \frac{1}{8} \left( \frac{N}{N_f N_m} \right) $$

$$ H_{S(t)} = 1 - (1 - H_{S(t-1)}) \frac{1}{4} \left( \frac{4N_f N_m - N}{N_f N_m} \right) - (1 - H_{S(t-2)}) \frac{1}{8} \left( \frac{N}{N_f N_m} \right) - \frac{1}{8} \left( \frac{N}{N_f N_m} \right) $$

<figure><img src="image_5.png" /></figure>

$$ H_{S(t)} = 1 - \frac{1}{4} \left( \frac{4N_f N_m - N}{N_f N_m} \right) + \frac{1}{4} H_{S(t-1)} \left( \frac{4N_f N_m - N}{N_f N_m} \right) - \frac{1}{8} \left( \frac{N}{N_f N_m} \right) + \frac{1}{8} H_{S(t-2)} \left( \frac{N}{N_f N_m} \right) - \frac{1}{8} \left( \frac{N}{N_f N_m} \right) $$

<figure><img src="image_6.png" /></figure>

$$ H_{S(t)} = H_{S(t-1)} \frac{1}{4} \left( \frac{4N_f N_m - N}{N_f N_m} \right) + H_{S(t-2)} \frac{1}{8} \left( \frac{N}{N_f N_m} \right) + 1 - \frac{1}{4} \left( \frac{4N_f N_m - N + N}{N_f N_m} \right) $$

<figure><img src="image_7.png" /></figure>

$$ (9) \qquad H_{S(t)} = H_{S(t-1)} \frac{4N_f N_m - N}{4N_f N_m} + H_{S(t-2)} \frac{N}{8N_f N_m} $$

Let $\lambda$ be:

$$ (10) \qquad \lambda = \frac{H_{S(t)}}{H_{S(t-1)}} $$

Assuming $\lambda$ to be constant from one generation to the next, and dividing all terms by $H_{S(t-1)}$, Equation 10 writes:

$$ \lambda = \frac{4N_f N_m - N}{4N_f N_m} + \frac{1}{\lambda} \frac{N}{8N_f N_m} $$

<figure><img src="image_8.png" /></figure>

$$ \lambda^2 - \frac{4N_f N_m - N}{4N_f N_m} \lambda = \frac{N}{8N_f N_m} $$
<figure><img src="image_9.png" /></figure>

$$ \lambda^2 - 2 \frac{4N_f N_m - N}{8N_f N_m} \lambda + \left( \frac{4N_f N_m - N}{8N_f N_m} \right)^2 = \frac{N}{8N_f N_m} + \left( \frac{4N_f N_m - N}{8N_f N_m} \right)^2 $$

<figure><img src="image_10.png" /></figure>

$$ \left( \lambda - \frac{4N_f N_m - N}{8N_f N_m} \right)^2 = \frac{N}{8N_f N_m} + \left( \frac{4N_f N_m - N}{8N_f N_m} \right)^2 $$

<figure><img src="image_11.png" /></figure>

$$ \lambda = \frac{4N_f N_m - N}{8N_f N_m} \pm \sqrt{\frac{N}{8N_f N_m} + \left( \frac{4N_f N_m - N}{8N_f N_m} \right)^2} $$

As $\lambda$ is positive, the single positive (leading) root of this equation is:

<figure><img src="image_12.png" /></figure>

$$ \lambda = \frac{4N_f N_m - N}{8N_f N_m} + \sqrt{\frac{N}{8N_f N_m} + \left( \frac{4N_f N_m - N}{8N_f N_m} \right)^2} $$

<figure><img src="image_13.png" /></figure>

$$ \lambda = \frac{1}{2} - \frac{N}{8N_f N_m} + \sqrt{\frac{N}{8N_f N_m} + \left( \frac{1}{2} - \frac{N}{8N_f N_m} \right)^2} $$

<figure><img src="image_14.png" /></figure>

$$ \lambda = \frac{1}{2} - \frac{N}{8N_f N_m} + \frac{1}{2} \sqrt{\frac{N}{2N_f N_m} + \left(1 - \frac{N}{4N_f N_m}\right)^2} $$

<figure><img src="image_15.png" /></figure>

$$ \lambda = \frac{1 - \frac{N}{4N_f N_m} + \sqrt{\frac{N}{2N_f N_m} + 1 + \left(\frac{N}{4N_f N_m}\right)^2 - 2\frac{N}{4N_f N_m}}}{2} $$

(11)

$$ \lambda = \frac{1 - \frac{N}{4N_f N_m} + \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2}}{2} $$

Note that the same results can be obtained with the leading eigenvalue of the transition matrix describing the evolution of genetic identities (Appendix 7).

For a monogenic panmictic population:

$$ H_{S(t)} = H_{S(t-1)} \left(1 - \frac{1}{2N_e}\right) $$

and in that case:

(12)

$$ \lambda = 1 - \frac{1}{2N_e} $$

We now need to combine Equations 11 and 12 to get:

$$ 1 - \frac{1}{2N_e} = \frac{1 - \frac{N}{4N_f N_m} + \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2}}{2} $$

<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /></figure>

$$ \frac{1}{2N_e} = 1 - \frac{1 - \frac{N}{4N_f N_m} + \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2}}{2} $$

<figure><img src="image_18.png" /></figure>

$$ \frac{1}{2N_e} = \frac{2 - 1 + \frac{N}{4N_f N_m} - \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2}}{2} $$

<figure><img src="image_19.png" /></figure>

$$ \frac{1}{N_e} = 1 + \frac{N}{4N_f N_m} - \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2} $$

<figure><img src="image_20.png" /></figure>

$$ N_e = \frac{4N_f N_m}{4N_f N_m + N - 4N_f N_m \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2}} $$

(13)

$$ N_e = \frac{\frac{4N_f N_m}{4N_f N_m \left[ 1 - \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2} \right] + N}} $$

## A new approximation

According to Taylor-MacLaurin's expansion series, $\sqrt{1 + X} \approx 1 + \frac{1}{2}X - \frac{1}{8}X^2$ (see Appendix 8 for a detailed proof). We can thus approximate the square root in equation 13 with:

$$ \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2} \approx 1 + \frac{1}{2}\left(\frac{N}{4N_f N_m}\right)^2 - \frac{1}{8}\left(\frac{N}{4N_f N_m}\right)^4 + \dots $$

The quantity $4N_f N_m$ is the lowest for the most uneven sex-ratios, like $N_m=1$ and $N_f=N-1$. In that case:

$$ -\frac{1}{8}\left(\frac{N}{4N_f N_m}\right)^4 + \dots = -\frac{1}{8}\left(\frac{N}{4(N-1)}\right)^4 + \dots = -\frac{1}{32}\left(\frac{N}{N-1}\right)^4 + \dots \ll 1 + \frac{1}{2}\left(\frac{N}{4N_f N_m}\right)^2 $$

We can then consider that:

$$ \sqrt{1 + \left(\frac{N}{4N_f N_m}\right)^2} \approx 1 + \frac{1}{2}\left(\frac{N}{4N_f N_m}\right)^2 $$

Eq 13 can thus write:

$$ N_e \approx \frac{4N_f N_m}{4N_f N_m \left[ 1 - 1 - \frac{1}{2} \left( \frac{N}{4N_f N_m} \right)^2 \right] + N} $$

<figure><img src="image_21.png" /></figure>
$$N_e \approx \frac{4N_f N_m}{-2N_f N_m \left(\frac{N}{4N_f N_m}\right)^2 + N}$$

<figure><img src="image_22.png" /></figure>

$$N_e \approx \frac{4N_f N_m}{N - \frac{N^2}{8N_f N_m}}$$

<figure><img src="image_23.png" /></figure>

$$(14) \qquad N_e \approx \frac{4N_f N_m}{N} \times \frac{1}{1 - \frac{N}{8N_f N_m}}$$

Using Taylor-MacLaurin again we can see that: $1/(1-X)=1+X+X^2+X^3+...$ (Appendix 8).

We can thus rewrite the second term of the denominator of equation 14:

$$\frac{1}{1 - \frac{N}{8N_f N_m}} \approx 1 + \frac{N}{8N_f N_m} + \left(\frac{N}{8N_f N_m}\right)^2 + \left(\frac{N}{8N_f N_m}\right)^3 + \dots$$

Using this in equation 14 yields:

$$N_e \approx \frac{4N_f N_m}{N} \times \frac{1}{1 - \frac{N}{8N_f N_m}} = \frac{4N_f N_m}{N} \left( 1 + \frac{N}{8N_f N_m} + \left( \frac{N}{8N_f N_m} \right)^2 + \left( \frac{N}{8N_f N_m} \right)^3 + \dots \right)$$

<figure><img src="image_24.png" /></figure>

$$N_e = \frac{4N_f N_m}{N} + \frac{1}{2} + \frac{1}{2} \frac{N}{8N_f N_m} + \frac{1}{2} \left( \frac{N}{8N_f N_m} \right)^2 + \frac{1}{2} \left( \frac{N}{8N_f N_m} \right)^3 + \dots$$

For the same reasons as those given above, in this equation, terms that are squared, cubed etc... can be neglected, and we then get:

$$(15) \qquad N_e \approx \frac{4N_f N_m}{N} + \frac{1}{2} + \frac{1}{4} \frac{N}{4N_f N_m}$$

Note that if we neglect $N/(16N_f N_m)$, we obtain equation 2.

For an even sex ratio, we get:

$$(16) \qquad N_e \approx N + \frac{1}{2} + \frac{1}{4N}$$

Equations 15 and 16 are a little different from Balloux's equations 18 and 10 (Balloux, 2004), respectively:

$$\begin{cases} N_e = \frac{4N_f N_m}{N} + \frac{1}{2} + \frac{1}{2N} \\ N_e \approx N + \frac{1}{2} + \frac{1}{2N} \end{cases}$$
The reasons for this discrepancy between these two sets of equations are unclear due to the lack of details in Balloux' paper. For his equation 10, Balloux simply cites Wright's book (Wright, 1969) without mentioning the page or the equation number. A glance at Wright's book only provided a stronger approximation (page 212, equation 8.4): $N_e=N+1/2$, without giving much details (but see Felsenstein's book page 266-267 (Felsenstein, 2019)). Appendix 6 provides a detailed proof for Balloux's equation 10 in the (simpler) case of even sex-ratios.

In Figure 1, it can be seen that the first approximation found in Wright's book (Equation 1), as in all population genetics textbooks, strongly underestimates $N_e$, except for very big populations (as expected), as compared to other approximations. Wright's second equation and Balloux's one seem to display an equivalently small bias, though in varying directions for Balloux's equation, depending on the sex-ratio. This can be seen with a study of the sign of $\Delta N_e$ (see Appendix 9), where we notice that with the unique valid root of $\Delta N_e$ ($SR_2$), Balloux's equation will provide an over-estimate when $SR>SR_2$, an under-estimate when $SR<SR_2$ and will be accurate when $SR=SR_2=3-2\sqrt{2}\approx0.1716$ (e.g. $SR\approx1/6$; $SR\approx4/23$). From there, it is easily deduced that, in Figure 1, positive values of $\Delta N_e$ correspond to $SR>SR_2$, negative ones to $SR<SR_2$, and close or very close to accuracy around this threshold (for instance 35/204 is very close to $SR_2$). This bias is very small when $N_e>10$ (Figure 1). Finally, the new simplified estimate proposed in Equation 15 perfectly fits to Equation 13, except for very small $N_e<4$ where a very small overestimate can be noticed (Figure 1).

<figure><img src="image_25.png" /><figcaption>Figure 1 - Comparisons of the performances of different approximations of effective population size ($N_e$) in dioecious populations with uneven (left) and even (right) sex-ratios, as compared to equation 13 (Eq 13) ($N_e$). Performance was measured as $\Delta e=(N_e-a-N_e)/N_e$, with Wright's equations 1 and 2 (Eq1 and Eq2 respectively), Balloux (Eq3), and the new simplified estimate of the present paper (Eq15).</figcaption></figure>

## Estimating the effective population size from Wright's $F_{IS}$

As seen above with equations 4, 5 and 6, heterozygote excesses expected in pangamic dioecious populations as computed by Wright's $F_{IS}$, can give access to an estimate of $N_e$ from genotypic data. In the following, we propose other $F_{IS}$ based estimates.

Let $H_{exp}$ and $H_{obs}$ be the expected (under Castle-Weinberg expectations) and the observed proportion of heterozygotes in the population, respectively. We can express these proportions as the probabilities of drawing at random two different alleles in the population $H_{S(t)}$ and in an individual $H_{I(t)}$ respectively at any generation t: $H_{exp}=H_{S(t)}$ and $H_{obs}=H_{I(t)}$. Finally, from equation 7, we can see that $H_{I(t)}=H_{S(t-1)}$. If we take Nei's parametric definition of $F_{IS}$ (Nei & Chesser, 1983):

$$F_{IS} = \frac{H_{exp} - H_{obs}}{H_{exp}}$$

<figure><img src="image_26.png" /></figure>
$$F_{IS} = 1 - \frac{H_{obs}}{H_{exp}}$$

<figure><img src="image_27.png" /></figure>

$$ (17) \qquad F_{IS} = 1 - \frac{H_{S(t-1)}}{H_{S(t)}} $$

We can combine equations 10 and 17 to obtain:

$$ (18) \qquad F_{IS} = 1 - \frac{1}{\lambda} $$

Now, combining equation 18 with equation 12 yields:

$$ F_{IS} = 1 - \frac{1}{1 - \frac{1}{2N_e}} $$

<figure><img src="image_28.png" /></figure>

$$ F_{IS} = \frac{1 - \frac{1}{2N_e} - 1}{1 - \frac{1}{2N_e}} $$

<figure><img src="image_29.png" /></figure>

$$ F_{IS} = \frac{\frac{1}{2N_e}}{1 - \frac{1}{2N_e}} $$

<figure><img src="image_30.png" /></figure>

$$ F_{IS} = -\frac{1}{2N_e - 1} $$

<figure><img src="image_31.png" /></figure>

$$ 2N_e F_{IS} - F_{IS} = -1 $$

<figure><img src="image_32.png" /></figure>

$$ N_e = \frac{F_{IS} - 1}{2F_{IS}} $$

<figure><img src="image_33.png" /></figure>

$$ (19) \qquad N_e = \frac{1}{2} - \frac{1}{2F_{IS}} $$

This result is the same as Equation 4 plus one half of an individual.

We can also express N as a function of FIS in a dioecious population with an even sex-ratio (Appendix 10):

$$ (20) \qquad N = -\frac{1+F_{IS}}{2F_{IS}} $$

Here, N can be called the effective number of breeders, which must not be confused with the effective population size. Notice that equation 20 differs from equation 4 by the subtraction of half an individual.

If we use Equation 20 in Equation 16 and rearrange the formula we get (Appendix 10):

$$ (21) \qquad N_e \approx -\frac{1}{2F_{IS}} - \frac{F_{IS}}{2(1+F_{IS})} $$

We can here note that, if we use equation 20 in equation 2 (Wright (1969), page 197), with an even sex-ratio, we obtain equation 4 (first FIS based estimate of Pudovkin et al. (1996)).

We have estimated FIS based Ns under various population sizes and sex-ratios, using the expected FIS computed in Appendix 11:
$$ (22) \qquad F_{IS} \approx \frac{N_m+N_f}{N_m+N_f+8N_mN_f} $$

In the Figure 2, it can be seen that, as compared to expected values (Equation 13) Equations 5 (Pudovkin et al., second equation) and 19 (simple equation of the present paper) tend to over-estimate $N_e$ as compared to other estimates, unless population size becomes big enough. Equations 4 (Pudovkin et al., first equation) and 6 (Balloux) present an equivalent bias but in the opposite direction (under-estimate and over-estimate, respectively). This bias is only visible for small $N_e$'s (i.e. $N_e<10$). Interestingly, the average of these two is exactly equation 21 of the present study, which fits perfectly with the expected one, with an extremely small over-estimate for the smallest $N_e<3$.

<figure><img src="image_34.png" /><figcaption>**Figure 2** - Comparisons of the performances of different $F_{IS}$ based (Equation 22) estimates of effective population size ($N_{e,a}$) in dioecious populations, as compared to equation 13 (Eq 13) ($N_e$), measured as $\Delta_e=(N_{e,a}-N_e)/N_e$, with Pudovkin et al. equations 1 and 2 (Eq 4 and Eq 5 respectively), Balloux (Eq 6), and the new simplified estimates of the present paper: the most simple (Eq 19), and the final one (Eq 21).</figcaption></figure>

## Discussion

Genetic drift can be influenced by several factors other than dioecy, population size and sex-ratio: reproductive variance, non-overlapping generations, changes in population size and/or sex ratio, subdivision and selection. Such complications may lead to very cumbersome algebra if one wanted to present a more general expression for the effective population size.

The main goals of the present study were, in decreasing order of importance. 1) to present an improved approximation of the classic case where only dioecy and/or uneven sex-ratio have an effect and compare it to previous approximations; 2) to present detailed derivations of both the old and new results that would be accessible to most readers; and 3) to provide an improved $F_{IS}$ estimate deriving from this new approximation, to be used in empirical population genetics studies instead of the two most often used methods: Pudovkin et al. (1996) and Balloux (2004) (11 and 8 citations per year respectively according to Google scholar's based research in the computer program Publish and Perish 8.8.4275.8412 (Harzing, 2007)).

We do not expect any natural population to closely follow the model we explored in the present work. Nevertheless, excluding reproductive system variation and selection, any combination of the factors mentioned above will tend to reduce $N_e$ and $F_{IS}$ accordingly, and most probably their variance, with not much harm to the average expectations. Partial selfing (e.g. deriving from some kinds of automictic
parthenogenesis of females (De Meeûs et al., 2007)), or partial sib-mating, should be easily detected, produce generalized heterozygote deficits and thus exclude our method. Clonal propagation (e.g. through a special kind of endomitotic automixis of females (De Meeûs et al., 2007)) should also be easy to detect (De Meeûs et al., 2006) and again the *F*<sub>IS</sub> based method would be dismissed. Selection is locus specific and should only affect one or two of the loci used, which consequently should be easy to detect and exclude. Subdivision can have two effects: if there is a Wahlund effect, this should be easy to detect (Manangwa et al., 2019); and if not, highly subdivided populations should exhibit effective sub-population sizes that are very close to the one that these would exhibit if totally isolated (because rare immigrants are not expected to display much influence), and if not, subsamples should all converge toward the total effective population size, which should be easily detected too (Ravel et al. (2023), see also Waples & Do (2010), and Waples & England (2011)). Additionally, in most situation, empiricists have no idea of the effective sex-ratio, of the scenarios regarding how generations overlap, or how population size fluctuates across generations. Consequently, complication of estimates will neither allow an easy understanding of the mathematical developments, which was an important goal of the present work, nor take into account with certainty the real scenarios that occurred in the population under investigation. This is why, even if the new estimate will hardly give significantly different values as compared to the previous ones, we think it is still better using the theoretical one that is closer to the exact expected value (equation 21) and interpretable on more sound biological means (see below).

The new approximation proposed in Equation 15 is equivalent to what is expected in large dioecious populations (Equation 1), plus half an individual, plus half of a coalesced individual in a large dioecious population. As far as we know, such added quantities were never discussed before (but see Felsenstein (2019), page 266). We can here try to provide some biological interpretation of such quantities. One half of heterozygosity is lost when an individual reproduces by selfing. In a panmictic population (i.e. monoecious) of size *N*, a proportion 1/*N* of gametes are produced by selfing (Rousset, 1996), in which case half the genes coalesce in the progeny (*P*<sub>CM</sub>=1/2*N*. This can happen in the *N* individuals. Hence (1/*N*)×(1/2)×*N*=1/2 individuals are produced with coalesced genes per generation through selfing in a WF population. This means that one half of such coalescent event does not happen when random selfing is forbidden, as it is necessarily the case in dioecious populations. Additionally, in very big dioecious populations, *N*<sub>e</sub>≈4*N<sub>f</sub>*N<sub>m</sub>/*N* (Equation 1). For any diploid population, the instantaneous probability of coalescence is *P*<sub>C</sub>=1/(2*N*<sub>e</sub>) (see Laporte & Charlesworth (2002), Equation 7; or Nomura (2008), Equation 3). Consequently, for a very big dioecious population, this probability becomes (Equation 1) *P*<sub>CBD</sub>=(1/2)×*N*/(4*N<sub>f</sub>*N<sub>m</sub>). The number of individuals concerned are those that inherited twice the same allele from their grand-parents, which is (1/*N*<sub>f</sub>)×(1/4)×*N*<sub>f</sub> for females and (1/*N*<sub>m</sub>)×(1/4)×*N*<sub>m</sub> for males, hence ½ individuals. This yields *P*<sub>CBD</sub>×(1/2) individuals. In small dioecious populations, such coalescent events hardly happen, because as long as polymorphism is kept, male and female parents that mate randomly can only rarely have sampled twice the same alleles from the same grand-parent. These two differences with 1) panmictic populations, and 2) big dioecious populations, may explain Equation 15. In other words, if we call *N*<sub>e_BD</sub> the effective population size of a very big dioecious population, *N*<sub>NCNS</sub> the number of individuals that cannot be coalescent due to the absence of selfing and *N*<sub>NCSD</sub> the number of individuals that cannot be coalescent in a dioecious population due to the limited number of possible mates, then, in small dioecious populations, the effective population size is *N*<sub>e_SD</sub> = *N*<sub>e_BD</sub> + *N*<sub>NCNS</sub> + *N*<sub>NCSD</sub>.

Interestingly, the highly sophisticated, and fairly complicated to compute, Balloux's equation (Balloux, 2004), Equation 3 in the present paper, did not perform better than Wright's second equation (Equation 2), and worse than our Equation 15. As shown in Appendix 6, the coalescence effective population size was obtained after neglecting different terms at several successive steps of the analytical process. To be as accurate as Equation 21, Equation 2 indeed requires *N*<sub>e</sub> > 10 and/or a sufficient number of generations (e.g. 10). As seen from the Figure 1, this seems to indeed apply as long as *N*<sub>e</sub>>12. No explanations were provided for the abstract notion of the coalescence effective population size and the way used to weight approximated coalescent times computed at different hierarchical levels (e.g. individuals, subpopulations, etc...). What we were interested in, in this paper, was to compute the local effective population size, i.e. the one that affects the speed of polymorphism loss within subpopulations. In that case, the eigenvalue effective population size may be the most accurate.

Regarding *F*<sub>IS</sub>-based estimates, the fact that Pudovkin et al., 2<sup>nd</sup> equation (Pudovkin et al. (1996), equation 5) did not perform well, probably comes from the confusion between effective population size
and the number of individuals (or of breeders), at different steps of the analytical procedure. Equation 19 provided similar results as equation 5, though with a slightly stronger bias and is thus too biased also. Pudovkin et al., second equation (equation 5 of the present manuscript), is quite popular in empirical population genetics studies, and is the one implemented in NeEstimator (Do et al., 2014). It presented underestimates, even when $N_e>20$. Balloux's equation (equation 6), also popular, suffered from an overestimation of $N_e$, in a symmetric position as compared to underestimations of Equation 4 (Pudovkin et al., first equation). For both, the biases are small, particularly so when $N_e>4$. Interestingly, following the steps described in Balloux's paper (Balloux, 2004), but replacing the coalescence approach by the leading eigenvalue approach, provided the most accurate $F_{IS}$-based estimate of the effective population size in dioecious populations (Equation 21). It appeared to exactly correspond to the average between Equations 4 (Pudovkin et al., first equation) and 6 (Balloux).

It is worth recalling that the $F_{IS}$-based estimates given in Equation 21 assumes an even sex ratio. Nevertheless, strongly biased sex-ratios will affect $F_{IS}$ accordingly and should not have strong consequences on the estimate of effective population sizes, as suggested by the Figure 2.

We may also bear in mind that although random mating was assumed, we did not specify any reproductive strategy (mono or polygamy). Indeed, Equation 8 assumes polygamy, but monogamy is known to lead to the same result as polygamy, as demonstrated pages 267-268 in Felsenstein's book (Felsenstein, 2019), and in Appendix 12. The only difference is that, in monogamous populations, the sex-ratio of individuals that reproduce is necessarily even. Consequently, monogamy can prevent a possible high variance in male mating success, which would reduce $N_e$. But monogamy cannot produce an increase of $N_e$ as compared to pangamic polygamy. In that sense, and everything else being equal, gibbons (which are monogamous) should preserve better their genetic diversity than gorillas, but not better than bonobos (assuming bonobos are pangamous).

It is worth mentioning that these computations were based on accurate (exact) measures of $F_{IS}$. Unbiased estimators of F-statistics (Weir & Cockerham, 1984) suffer from large variances (Robertson & Hill, 1984). It is thus likely that deviation from the real value will have a large impact on $F_{IS}$-based estimates of $N_e$, especially for small expected ("real") values. It can be seen that $N_{e-Eq4} < N_{e-Eq21} < N_{e-Eq6} < N_{e-Eq5} < N_{e-Eq19}$. From there, it can be expected that with small underestimations of $F_{IS}$, $N_{e-Eq6}$ will be closer to the real value; then $N_{e-Eq5}$ for bigger underestimations, and finally $N_{e-Eq19}$ for the strongest underestimations. On the contrary, overestimations of $F_{IS}$ will necessarily lead $N_{e-Eq4}$ to stay closer to the real $N_e$. However, the differences are expected to be quite small, particularly so as compared to Pudovkin 1 (Equation 4) and Balloux (Equation 6), especially for the highest values (e.g. $N_e>6$): Nevertheless, not knowing what the real $F_{IS}$ should be, it is probably wiser using the less biased estimate, i.e. $N_{e-Eq21}$.

It is also worth mentioning that $F_{IS}$ should be estimated from adults, as the genetic structure in immature individuals may considerably differ from the one they would display in the pool of adults that survived.

The fact that our equation 21 outperformed other equations for $N_e<4-6$ may suggest strong limitations in the practical applicability of this performance since such systems may be expected to quickly undergo extinction. In addition to the fact that it is generally preferable to work with the most accurate equation, these results are likely to be especially pertinent for certain types of biological systems that are able to persist for extended periods despite having very small effective population sizes. For instance, the populations of the parasitoid wasp *Nasonia vitripensis*, depending on the distribution of its host (parasitic flies), often display systematic mating of females with their brothers (Werren, 1980). For the mite *Varroa destructor*, a female enters a brood cell, which she caps, where she feeds on the bee larva and then gives birth to a haploid male, which later mates with its diploid sisters, laid by the mother from fertilized eggs from a previous mating that occurred before the colonization of the brood cell, or after mating with her son (which may happen for 30% of females) (Beaurepaire et al., 2017; Häußermann et al., 2020). In both cases, males are produced by arrhenotokous parthenogenesis (unfertilized haploid eggs), meaning that many populations of these organisms probably display very small $N_e$, and even smaller than 1 in some instances. We did not undertake an extensive review of similar cases, as it is not in the scope of the present paper, but such kind of situations may not be rare in dioecious parasitic organisms like parasitoid hymenoptera, mites or nematodes.

According to recent papers based on simulations (i.e. perfect data), $F_{IS}$-based single sample (or subsample) estimates of $N_e$ are not the most accurate (Wang, 2009; Do et al., 2014; Wang, 2016). According
to Do et al. (2014), the linkage disequilibrium (LD) based estimate (Waples, 2006), appeared to perform better than the co-ancestry method (CoA) (Nomura, 2008) and the *F*<sub>IS</sub>-based method (Equation 5) (Pudovkin et al., 1996). According to more recent simulation studies (Wang, 2016), the sibship frequency based estimate (SF) (Wang, 2009) seemed to provide more accurate results than the previous ones. No comparison was ever made with an alternative method based on one and two locus identity measures (1&2LI) (Vitalis & Couvet, 2001c, 2001b), implemented by the software Estim 1.2.2 (Vitalis, 2002), updated from Estim 1 (Vitalis & Couvet, 2001a). Based on simulations, the 1&2LI method provided accurate (though slightly underestimated) results, especially when more than four loci were used (Vitalis & Couvet, 2001b). Again, no simulation study exhaustively compared all available one-sample estimates. This would require replicated simulations of different scenarios of population structure (Island or stepping stone models with varying subpopulation number, sub-population sizes and immigration rates), different kinds of loci (microsatellite like or SNP like loci) with varying number of loci, number of alleles and mutation rates, and with or without amplification problems (null alleles, stuttering, short allele dominance or allelic dropouts), and varying sampling strategies. A comparison with temporal methods (Nei & Tajima, 1981; Pollak, 1983; Wang & Whitlock, 2003; Jorde & Ryman, 2007) might also prove interesting, though the number of generations between two samples of the same site will add another relevant parameter to explore (Waples & Do, 2010). This will obviously require much more work to undertake, which is beyond the scope of the present paper.

We nevertheless undertook a quick simulation study with Easypop (Balloux, 2001). We simulated single isolated and randomly mating dioecious populations, with varying sex-ratio, at 100 independent loci with a KAM model of mutation with *K*=100 possible allelic states and a mutation rate of *u*=0.00001, and 100 generations. All simulations started with maximum diversity. We then computed effective population sizes. We computed *F*<sub>IS</sub> with Fstat (Goudet, 1995). For these simulations, most of the averaged *F*<sub>IS</sub> across loci were positive and therefore could not be used to estimate *N*<sub>e</sub>. We then preferred computing the average across loci displaying a negative *F*<sub>IS</sub>. For NeEstimator analyses (LD and Coancestries), we assumed polygamy and kept estimates excluding alleles less frequent than 5% (LD method). For Estim (1&2LI), we assumed panmixia. For Colony (SF), we generated data using Create (Coombs et al., 2008) and assumed polygamy and some inbreeding, as this may occur at unknown level in real data. Figure 3 illustrates what kind of variations could be observed from one parameter set to the other and from one method to the other. It suggests that some kind of average across methods may allow grasping the range of actual effective population sizes of sub-populations from genotyped sub-samples.

<figure><img src="image_35.png" /><figcaption>**Figure 3** - Effective population size estimates (*N*<sub>e</sub>) with five different methods as compared to the expected value (Equation 22), for different simulated populations with varying numbers of females (*N*<sub>f</sub>) and males (*N*<sub>m</sub>).</figcaption></figure>
Regarding real data, quoting Nomura, "combined estimate of several independent estimates is expected to improve the precision of separate estimates" (Nomura, 2008). For each method, one could compute the average $N_e$ across subsamples of the same population, ignoring undefined values (negative or infinite), note the maximum and minimum values obtained and keep the number of usable values as a weight. Finally, the grand average (across methods) and average minimum and maximum, all weighted by the number of usable values obtained in each method, could be computed. For more clarity, a template of this method can be found in the file "TemplateRhipicephalusFstatResNeFiveMethods.xlsx", coming from the analysis of cattle tick populations from New-Caledonia (De Meeûs et al., 2010). With all data, average $N_e \approx 120$ in minimax$\approx[80, 200]$. When excluding the two most extreme values $N_e \approx 50$ in minimax$\approx[10, 110]$. Using the harmonic mean, as suggested by Nomura (Nomura, 2008), $N_e \approx 20$ in minimax$\approx[10, 30]$. Simulation studies could be used to identify an estimator that more accurately approximates the eigenvalue effective population size of genotyped populations.

Temporal data are rarely available (but see Palstra & Ruzzante (2008)), but when these are, they give access to different estimates, which may be usefully included in the computations of averages and magnitude of variations.

Undefined $N_e$ may correspond to very big values. Thus, ignoring these may lead to underestimates. They may also correspond to the variance of estimate of the parameter used, like $F_{IS}$, as mentioned above. This possible flaw may be attenuated by the use of repeated subsamples and independent loci. Waples and Do proposed to include negative $N_e$ as such in the computation of an harmonic mean, with weights proportional to reciprocals of variances (Waples & Do, 2010). Nevertheless, on the tick data set, this strategy ended with globally negative (and then unsound) values for these populations (not shown), which are expected to display important population sizes (i.e. $120 \le N_e \le 1200$ (Koffi et al., 2006)). As already discussed, this will require a more thorough exploration through simulations of various kinds of populations.

To conclude, even if the differences with some other equations are not very big, the new approximation proposed here appeared almost perfect and biologically relatively sound, and, wherever it is used for, we suggest to use it instead of the previous and more biased estimates found in the literature.

## Acknowledgements

Thierry de Meeûs would like to thank François Rousset, François Ziegler and members of the Maxima-discuss list for recurrently answering to his queries and for giving him tips that really unblocked the advances of the present work. Special thanks to Joseph Felsenstein who gave the idea of using Taylor-MacLaurin series, which triggered the transformation of a small paragraph of a book into a full article. We also are grateful to Robin Waples, Jay Taylor and two anonymous referees, whose comments considerably helped improving this manuscript, as Myriam Heuertz for her significant editorial inputs. This work was entirely financed by the Institut de Recherche pour le Développement (IRD) (recurring subsidies of UMR Intertryp). Preprint version 6 of this article has been peer-reviewed and recommended by Peer Community In Evolutionary Biology (https://doi.org/10.24072/pci.evolbiol.100651, Baer (2023)).

## Data availability

Maxima scripts are provided in the section Scripts at the end of the manuscript. Simulations used for the graphic of Figure 3 are available in the file "Simulations-1st-pop.zip". This file, the template "TemplateRhipicephalusFstatResNeFiveMethods.xlsx", with the corresponding dataset "BoophilusAdultsDataCattle.txt", and the file VarDif.pdf (computation of the variance of a difference) are also available online at https://doi.org/10.5281/zenodo.7945369 (De Meeûs, 2023).

## Conflict of interest disclosure

The authors declare that they have no financial conflict of interest with the content of this article. Thierry de Meeûs is one of the PCI Infections administrators.
References

Baer C (2023) All you ever wanted to know about $N_e$ in one handy place. *Peer Community in Evolutionary Biology*, **100651**. https://doi.org/10.24072/pci.evolbiol.100651

Balloux F (2001) EASYPOP (version 1.7): A computer program for population genetics simulations. *Journal of Heredity*, **92**, 301-302. https://doi.org/10.1093/jhered/92.3.301

Balloux F (2004) Heterozygote excess in small populations and the heterozygote-excess effective population size. *Evolution*, **58**, 1891-1900. https://doi.org/10.1554/03-692

Balloux F, Lehmann L (2003) Random mating with a finite number of matings. *Genetics*, **165**, 2313-2315. https://doi.org/10.1093/genetics/165.4.2313

Balloux F, Lehmann L, De Meeûs T (2003) The population genetics of clonal and partially clonal diploids. *Genetics*, **164**, 1635-1644. https://doi.org/10.1093/genetics/164.4.1635

Beaurepaire AL, Krieger KJ, Moritz RFA (2017) Seasonal cycle of inbreeding and recombination of the parasitic mite *Varroa destructor* in honeybee colonies and its implications for the selection of acaricide resistance. *Infection Genetics and Evolution*, **50**, 49-54. https://doi.org/10.1016/j.meegid.2017.02.011

Castle WE (1903) The laws of heredity of Galton and Mendel, and some laws governing race improvement by selection. *Proceedings of the American Academy of Arts and Sciences*, **39**, 223-242. https://doi.org/10.2307/20021870

Coombs JA, Letcher BH, Nislow KH (2008) CREATE: a software to create input files from diploid genotypic data for 52 genetic software programs. *Molecular Ecology Resources*, **8**, 578-580. https://doi.org/10.1111/j.1471-8286.2007.02036.x

Crow JF, Kimura M (1970) *An Introduction to Population Genetics Theory*. The Blackburn Press, Caldwell, New-Jersey.

De Meeûs T (2023) Dataset of the article "A new and almost perfectly accurate approximation of the eigenvalue effective population size of a dioecious population: comparisons with other estimates and detailed proofs" (Version 2) [Data set]. *Zenodo*. https://doi.org/10.5281/zenodo.7945369

De Meeûs T, Chan CT, Ludwig JM, Tsao JI, Patel J, Bhagatwala J, Beati L (2021) Deceptive combined effects of short allele dominance and stuttering: an example with *Ixodes scapularis*, the main vector of Lyme disease in the U.S.A. *Peer Community Journal*, **1**, e40. https://doi.org/10.24072/pcjournal.34

De Meeûs T, Koffi BB, Barré N, de Garine-Wichatitsky M, Chevillon C (2010) Swift sympatric adaptation of a species of cattle tick to a new deer host in New-Caledonia. *Infection Genetics and Evolution*, **10**, 976-983. https://doi.org/10.1016/j.meegid.2010.06.005

De Meeûs T, Lehmann L, Balloux F (2006) Molecular epidemiology of clonal diploids: A quick overview and a short DIY (do it yourself) notice. *Infection Genetics and Evolution*, **6**, 163-170. https://doi.org/10.1016/j.meegid.2005.02.004

De Meeûs T, Prugnolle F, Agnew P (2007) Asexual reproduction: Genetics and evolutionary aspects. *Cellular and Molecular Life Sciences*, **64**, 1355-1372. https://doi.org/10.1007/s00018-007-6515-2

Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR (2014) NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size ($N_e$) from genetic data. *Molecular Ecology Resources*, **14**, 209-214. https://doi.org/10.1111/1755-0998.12157

Ewens WJ (2004) *Mathematical Population Genetics: I. Theoretical Introduction, 2nd Edition*. Interdisciplinary Applied Mathematics, Volume 27 (Vol. 27). Springer, New York.

Felsenstein J (2019) *Theoretical Evolutionary Genetics*. Department of Genome Sciences and Department of Biology, University of Washington, Seattle, Washington.

Goudet J (1995) FSTAT (Version 1.2): A computer program to calculate F-statistics. *Journal of Heredity*, **86**, 485-486. https://doi.org/10.1093/oxfordjournals.jhered.a111627

Harzing AW (2007) Publish or Perish, available from https://harzing.com/resources/publish-or-perish.

Häußermann CK, Giacobino A, Munz R, Ziegelmann B, Palacio MA, Rosenkranz P (2020) Reproductive parameters of female *Varroa destructor* and the impact of mating in worker brood of *Apis mellifera*. *Apidologie*, **51**, 342-355. http://doi.org/10.1007/s13592-019-00713-9

Horn RA, Johnson CR (2013) *Matrix Analysis, Second Edition*. Cambridge University Press, Cambridge, UK.

Jorde PE, Ryman N (2007) Unbiased estimator for genetic drift and effective population size. *Genetics*, **177**, 927-935. https://doi.org/10.1534/genetics.107.075481
Koffi BB, De Meeûs T, Barré N, Durand P, Arnathau C, Chevillon C (2006) Founder effects, inbreeding and effective sizes in the Southern cattle tick: the effect of transmission dynamics and implications for pest management. *Molecular Ecology*, **15**, 4603-4611. https://doi.org/10.1111/j.1365-294x.2006.03098.x

Laporte V, Charlesworth B (2002) Effective population size and population subdivision in demographically structured populations. *Genetics*, **162**, 501-519. https://doi.org/10.1093/genetics/162.1.501

Manangwa O, De Meeûs T, Grébaut P, Segard A, Byamungu M, Ravel S (2019) Detecting Wahlund effects together with amplification problems : cryptic species, null alleles and short allele dominance in Glossina pallidipes populations from Tanzania. *Molecular Ecology Resources*, **19**, 757-772. https://doi.org/10.1111/1755-0998.12989

Nei M, Chesser RK (1983) Estimation of fixation indices and gene diversities. *Annals of Human Genetics*, **47**, 253-259. https://doi.org/10.1111/j.1469-1809.1983.tb00993.x

Nei M, Tajima F (1981) Genetic drift and estimation of effective population size. *Genetics*, **98**, 625-640. https://doi.org/10.1093/genetics/98.3.625

Nomura T (2008) Estimation of effective number of breeders from molecular coancestry of single cohort sample. *Evolutionary Applications*, **1**, 462-474. https://doi.org/10.1111/j.1752-4571.2008.00015.x

Palstra FP, Ruzzante DE (2008) Genetic estimates of contemporary effective population size: what can they tell us about the importance of genetic stochasticity for wild population persistence? *Molecular Ecology*, **17**, 3428-3447. https://doi.org/10.1111/j.1365-294X.2008.03842.x

Pollak E (1983) A new method for estimating the effective population size from allele frequency changes. *Genetics*, **104**, 531-548. https://doi.org/10.1093/genetics/104.3.531

Pudovkin AI, Zaykin DV, Hedgecock D (1996) On the potential for estimating the effective number of breeders from heterozygote excess in progeny. *Genetics*, **144**, 383-387. https://doi.org/10.1093/genetics/144.1.383

Ravel S, Mahamat MH, Ségard A, Argiles-Herrero R, Bouyer J, Rayaisse J-B, Solano P, Mollo BG, Pèka M, Darnas J, Belem AMG, Yoni W, Noûs C, De Meeûs T (2023) Population genetics of *Glossina fuscipes fuscipes* from southern Chad. *Peer Community Journal*, **3**, e31. https://doi.org/10.24072/pcjournal.257

Robertson A (1965) The interpretation of genotypic ratios in domestic animal populations. *Animal Production*, **7**, 319-324. https://doi.org/10.1017/S0003356100025770

Robertson A, Hill WG (1984) Deviations from Hardy-Weinberg proportions - Sampling variances and use in estimation of inbreeding coefficients. *Genetics*, **107**, 703-718. https://doi.org/10.1093/genetics/107.4.703

Rousset F (1996) Equilibrium values of measures of population subdivision for stepwise mutation processes. *Genetics*, **142**, 1357-1362. https://doi.org/10.1093/genetics/142.4.1357

Rousset F (2004) *Genetic Structure and Selection in Subdivided Populations*. Princeton University Press, Princeton.

Vitalis R (2002) Estim 1.2-2: a computer program to infer population parameters from one- and two-locus gene identity probabilities, updated from Vitalis and Couvet (2001), *Molecular Ecology Notes*, **1**, 354-356.

Vitalis R, Couvet D (2001a) ESTIM 1.0: a computer program to infer population parameters from one- and two-locus gene identity probabilities. *Molecular Ecology Notes*, **1**, 354-356. https://doi.org/10.1046/j.1471-8278.2001.00086.x

Vitalis R, Couvet D (2001b) Estimation of effective population size and migration rate from one- and two-locus identity measures. *Genetics*, **157**, 911-925. https://doi.org/10.1093/genetics/157.2.911

Vitalis R, Couvet D (2001c) Two-locus identity probabilities and identity disequilibrium in a partially selfing subdivided population. *Genetical Research*, **77**, 67-81. https://doi.org/10.1017/s0016672300004833

Vodopivec A (2017) wxMaxima, a document based interface for the computer algebra system Maxima, distributed under the GPL license, donloadable at https://wxmaxima-developers.github.io/wxmaxima/. version **17.10.1**

Wang JL (2009) A new method for estimating effective population sizes from a single sample of multilocus genotypes. *Molecular Ecology*, **18**, 2148-2164. https://doi.org/10.1111/j.1365-294X.2009.04175.x

Wang JL (2016) A comparison of single-sample estimators of effective population sizes from genetic marker data. *Molecular Ecology*, **25**, 4692-4711. https://doi.org/10.1111/mec.13725

Wang JL, Whitlock MC (2003) Estimating effective population size and migration rates from genetic samples over space and time. *Genetics*, **163**, 429-446. https://doi.org/10.1093/genetics/163.1.429
Waples RS (2006) A bias correction for estimates of effective population size based on linkage disequilibrium at unlinked gene loci. *Conservation Genetics, 7*, 167-184. https://doi.org/10.1007/s10592-005-9100-y

Waples RS, Do C (2010) Linkage disequilibrium estimates of contemporary $N_e$ using highly variable genetic markers: a largely untapped resource for applied conservation and evolution. *Evolutionary Applications, 3*, 244-262. https://doi.org/10.1111/j.1752-4571.2009.00104.x

Waples RS, England PR (2011) Estimating contemporary effective population size on the basis of linkage disequilibrium in the face of migration. *Genetics, 189*, 633-644. https://doi.org/10.1534/genetics.111.132233

Weinberg W (1908) Über den Nachweis der Verebung beim Menschen. *Jahresheft des Vereins fur Vaterlundische Naturkunde in Wurttemberg, 64*, 368-382. https://archive.org/details/b30613000/page/370/mode/2up

Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. *Evolution, 38*, 1358-1370. https://doi.org/10.1111/j.1558-5646.1984.tb05657.x

Werren JH (1980) Sex ratio adaptations to local mate competition in a parasitic wasp. *Science, 208*, 1157-1159. https://doi.org/10.1126/science.208.4448.1157

Wright S (1965) The interpretation of population structure by F-statistics with special regard to system of mating. *Evolution, 19*, 395-420. https://doi.org/10.1111/j.1558-5646.1965.tb01731.x

Wright S (1969) *Evolution and the genetics of Populations Volume 2: The Theory of Gene Frequencies*. The University of Chicago Press, Chicago.

Appendices

**Appendix 1: Matrix multiplication, identity matrix, matrix determinant and matrix inversion.**

By convention, matrices and vectors appear in bold, while scalars write in italics, and matrix multiplication is noted by a point.

Let matrix **A** and vector **x** be:

$$ \mathbf{A} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \quad \text{and} \quad \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} $$

Multiplying **A** by **x** yields a new vector:

$$ (A1-1) \qquad \mathbf{A} \cdot \mathbf{x} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} ax_1 + bx_2 \\ cx_1 + dx_2 \end{pmatrix} $$

With $\mathbf{B}=\begin{pmatrix} e & f \\ g & h \end{pmatrix}$, then:

$$ (A1-2) \qquad \mathbf{A} \cdot \mathbf{B} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} e & f \\ g & h \end{pmatrix} = \begin{pmatrix} ae + bg & af + bh \\ ce + dg & cf + dh \end{pmatrix} $$

Please, note that most of the time **A**.**B**≠**B**.**A**, since:

$$ \mathbf{B} \cdot \mathbf{A} = \begin{pmatrix} e & f \\ g & h \end{pmatrix} \begin{pmatrix} a & b \\ c & d \end{pmatrix} = \begin{pmatrix} ae + cf & de + df \\ ag + ch & cg + dh \end{pmatrix} $$

The identity matrix **I** must verify:

$$ (A1-3) \qquad \left\{ \begin{array}{l} \mathbf{A} \cdot \mathbf{I} = \mathbf{I} \cdot \mathbf{A} = \mathbf{A} \\ \mathbf{I} \cdot \mathbf{x} = \mathbf{x} \end{array} \right. $$

$$ \text{Let } \mathbf{I} = \begin{pmatrix} I_{11} & I_{12} \\ I_{21} & I_{22} \end{pmatrix}. $$
Then:

$$ \begin{pmatrix} a & b \\ c & d \end{pmatrix} \cdot \begin{pmatrix} I_{11} & I_{12} \\ I_{21} & I_{22} \end{pmatrix} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} $$

<figure><img src="image_36.png" /></figure>

$$ \begin{pmatrix} aI_{11} + bI_{21} & aI_{12} + bI_{22} \\ cI_{11} + dI_{21} & cI_{12} + dI_{22} \end{pmatrix} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} $$

<figure><img src="image_37.png" /></figure>

$$ \begin{cases} aI_{11} + bI_{21} = a \\ aI_{12} + bI_{22} = b \\ cI_{11} + dI_{21} = c \\ cI_{12} + dI_{22} = d \end{cases} $$

<figure><img src="image_38.png" /></figure>

$$ \begin{cases} I_{11} = \frac{a - bI_{21}}{a} \\ I_{12} = \frac{b(1 - I_{22})}{a} \\ c \frac{a - bI_{21}}{a} + dI_{21} = c \\ c \frac{b(1 - I_{22})}{a} + dI_{22} = d \end{cases} $$

<figure><img src="image_39.png" /></figure>

$$ (A1-4) \qquad \begin{cases} I_{11} = \frac{a-bI_{21}}{a} \\ I_{12} = \frac{b(1-I_{22})}{a} \\ I_{21} \left(d - \frac{bc}{a}\right) = c - \frac{ca}{a} \\ I_{22} \left(d - \frac{bc}{a}\right) = d - \frac{bc}{a} \end{cases} $$

For $a \neq 0$, then equation A1-4 writes:

$$ \begin{cases} I_{11} = \frac{a - bI_{21}}{a} \\ I_{12} = \frac{b(1 - I_{22})}{a} \\ I_{21} \left( \frac{ad - bc}{a} \right) = 0 \\ I_{22} \left( \frac{ad - bc}{a} \right) = \frac{ad - bc}{a} \end{cases} $$

<figure><img src="image_40.png" /></figure>

$$ (A1-5) \qquad \begin{cases} I_{11} = \frac{a-bI_{21}}{a} \\ I_{12} = \frac{b(1-I_{22})}{a} \\ I_{21}(ad-bc) = 0 \\ I_{22} = \frac{ad-bc}{ad-bc} \end{cases} $$

Expression $ad-bc$ is the determinant of matrix $A$, $Det(A)$. Equation A1-5 has a solution only if $Det(A) \neq 0$, in which case A1-5 writes:
$$ \begin{cases} I_{11} = 1 \\ I_{12} = 0 \\ I_{21} = 0 \\ I_{22} = 1 \end{cases} $$

$$ \mathbf{I} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} $$

<figure><img src="image_41.png" /></figure>

From there, it is easy to understand that the identity matrix of any dimension $n$ is a squared matrix with diagonal numbers equal to 1 and others to 0:

$$ \mathbf{I}_n = \begin{pmatrix} 1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 1 \end{pmatrix} $$

$$ \text{For } n=3, \mathbf{I}_3 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

We now have the tools to find the inverse of matrix $\mathbf{A}$, $\mathbf{A}^{-1}$, which must verify $\mathbf{A}^{-1} \cdot \mathbf{A}=\mathbf{I}$. If we set

$$ \mathbf{A}^{-1} = \begin{pmatrix} \alpha & \beta \\ \gamma & \delta \end{pmatrix}, \text{ we can write:} $$

$$ \begin{pmatrix} \alpha & \beta \\ \gamma & \delta \end{pmatrix} \cdot \begin{pmatrix} a & b \\ c & d \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} $$

<figure><img src="image_42.png" /></figure>

$$ \begin{cases} \alpha a + \beta c = 1 \\ \alpha b + \beta d = 0 \\ \gamma a + \delta c = 0 \\ \gamma b + \delta d = 1 \end{cases} $$

<figure><img src="image_43.png" /></figure>

$$ \begin{cases} \alpha = \frac{1 - \beta c}{a} \\ \frac{1 - \beta c}{a} b + \beta d = 0 \\ \gamma = -\frac{\delta c}{a} \\ -\frac{\delta c}{a} b + \delta d = 1 \end{cases} $$

<figure><img src="image_44.png" /></figure>

$$ \begin{cases} \alpha = \frac{1 - \beta c}{a} \\ \beta \left(d - \frac{bc}{a}\right) = -\frac{b}{a} \\ \gamma = -\frac{\delta c}{a} \\ \delta \left(d - \frac{cb}{a}\right) = 1 \end{cases} $$

<figure><img src="image_45.png" /></figure>
<figure><img src="image_46.png" /></figure>

$$ \left\{ \begin{aligned} \alpha &= \frac{1 - \beta c}{a} \\ \beta \left(d - \frac{bc}{a}\right) &= \frac{-\frac{b}{a}}{\frac{ad - bc}{a}} \\ \gamma &= -\frac{\delta c}{a} \\ \delta \left(d - \frac{cb}{a}\right) &= \frac{1}{\frac{ad - bc}{a}} \end{aligned} \right. $$

(A1-6)

$$ \left\{ \begin{aligned} \alpha &= \frac{1 - \beta c}{a} \\ \beta &= \frac{-\frac{b}{a}}{\frac{ad - bc}{a}} \\ \gamma &= -\frac{\delta c}{a} \\ \delta &= \frac{1}{\frac{ad - bc}{a}} \end{aligned} \right. $$

For $a \neq 0$, equation A1-6 becomes:

<figure><img src="image_47.png" /></figure>

$$ \left\{ \begin{aligned} \alpha &= \frac{1 - \beta c}{a} \\ \beta &= \frac{-b}{ad - bc} \\ \gamma &= -\frac{\delta c}{a} \\ \delta &= \frac{a}{ad - bc} \end{aligned} \right. $$

<figure><img src="image_48.png" /></figure>

$$ \left\{ \begin{aligned} \alpha &= \frac{1 + \frac{b}{ad - bc} c}{a} \\ \beta &= \frac{-b}{ad - bc} \\ \gamma &= - \frac{ad - bc c}{a} \\ \delta &= \frac{a}{ad - bc} \end{aligned} \right. $$

<figure><img src="image_49.png" /></figure>

$$ \left\{ \begin{aligned} \alpha &= \frac{ad - bc + bc}{\frac{ad - bc}{a}} \\ \beta &= \frac{-b}{ad - bc} \\ \gamma &= - \frac{ac}{ad - bc} \\ \delta &= \frac{a}{ad - bc} \end{aligned} \right. $$
$$ 
\left\{ 
\begin{aligned}
\alpha &= \frac{d}{ad - bc} \\
\beta &= \frac{-b}{ad - bc} \\
\gamma &= \frac{-c}{ad - bc} \\
\delta &= \frac{a}{ad - bc}
\end{aligned}
\right. 
$$

Consequently, the reverse of matrix A is:

$$ A^{-1} = \frac{1}{ad - bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} $$

We know that ad-bc=Det(A), thus:

$$ (A1-7) \qquad A^{-1} = \frac{1}{\det(A)} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} $$

When Det(A)=0, A is singular. When Det(A)≠0, A is nonsingular. A nonsingular matrix is necessarily squared.

To compute the inverse of a 3×3 matrix, it is easier using a mathematical software as Maxima (Vodopivec, 2017) (command invert(A)).

## Appendix 2: eigenvalues and eigenvectors

For the sake of space saving and simplicity, we will take the example of a 2×2 matrix X and a two lines column vector e:

$$ 
\left\{ 
\begin{aligned}
\mathbf{X} &= \begin{pmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \end{pmatrix} \\
\mathbf{e} &= \begin{pmatrix} e_1 \\ e_2 \end{pmatrix}
\end{aligned}
\right. 
$$

If λᵢ is an eigenvalue and eᵢ the corresponding eigenvector of matrix X, then they must satisfy the equation: X.eᵢ=λᵢ.eᵢ. We can translate this into the system of equations:

$$ 
\begin{cases} 
x_{11} \times e_1 + x_{12} \times e_2 = \lambda e_1 \\ 
x_{21} \times e_1 + x_{22} \times e_2 = \lambda e_2 
\end{cases} 
$$

Excluding the trivial solution where e₁=e₂=0, we can rewrite the preceding equations as:

$$ 
\left\{ 
\begin{aligned}
x_{11} + x_{12} \times \frac{e_2}{e_1} &= \lambda \\
x_{21} + x_{22} \times \frac{e_2}{e_1} &= \lambda \frac{e_2}{e_1}
\end{aligned}
\right. 
$$

$$ 
\left\{ 
\begin{aligned}
x_{11} + x_{12} \times \frac{e_2}{e_1} &= \lambda \\
x_{21} &= \frac{e_2}{e_1}(\lambda - x_{22})
\end{aligned}
\right. 
$$

$$ 
\left\{ 
\begin{aligned}
x_{11} + x_{12} \times \frac{e_2}{e_1} &= \lambda \\
\frac{x_{21}}{\lambda - x_{22}} &= \frac{e_2}{e_1}
\end{aligned}
\right. 
$$

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>
<figure><img src="image_53.png" /></figure>

$$ \left\{ \begin{array}{l} x_{11} + x_{12} \times \frac{x_{21}}{\lambda - x_{22}} = \lambda \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

$$ \left\{ \begin{array}{l} \frac{(\lambda - x_{22})x_{11} + x_{12} \times x_{21}}{\lambda - x_{22}} = \lambda \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_54.png" /></figure>

$$ \left\{ \begin{array}{l} x_{11} \times \lambda - x_{11}x_{22} + x_{12} \times x_{21} = \lambda^2 - x_{22} \times \lambda \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_55.png" /></figure>

$$ \left\{ \begin{array}{l} \lambda^2 - (x_{11} + x_{22}) \times \lambda = -x_{11}x_{22} + x_{12} \times x_{21} \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_56.png" /></figure>

$$ \left\{ \begin{array}{l} \lambda^2 - x_{11} \times \lambda - x_{22} \times \lambda + x_{11}x_{22} - x_{12} \times x_{21} = 0 \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_57.png" /></figure>

$$ \left\{ \begin{array}{l} -\lambda(x_{11} - \lambda) + x_{22}(x_{11} - \lambda) - x_{21} \times x_{12} = 0 \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_58.png" /></figure>

$$ (A2-1) \qquad \left\{ \begin{array}{l} (x_{11} - \lambda)(x_{22} - \lambda) - x_{21} \times x_{12} = 0 \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

Knowing that the determinant of a matrix A Det(A)=Det $\begin{pmatrix} a & b \\ c & d \end{pmatrix}$=ad-cd, Eq A2-1 writes:

<figure><img src="image_59.png" /></figure>

$$ \left\{ \begin{array}{l} \text{Det} \begin{pmatrix} x_{11} - \lambda & x_{12} \\ x_{21} & x_{22} - \lambda \end{pmatrix} = 0 \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

$$ \left\{ \begin{array}{l} \text{Det} \left( \begin{pmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \end{pmatrix} - \begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix} \right) = 0 \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_60.png" /></figure>

$$ (A2-2) \qquad \left\{ \begin{array}{l} \text{Det} \left( \begin{pmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \end{pmatrix} - \lambda \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \right) = 0 \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

The matrix $I=\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$ is called the identity matrix (Appendix 1).
The first line of EqA2-2 writes $\text{Det}(A)-\lambda.I=0$ and corresponds to the so called characteristic equation of matrix A. We can solve EqA2-2:

$$ \left\{ \begin{array}{l} \text{Det}\left( \begin{pmatrix} x_{11} - \lambda & x_{12} \\ x_{21} & x_{22} - \lambda \end{pmatrix} \right) = 0 \\ \\ \frac{x_2}{x_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_61.png" /></figure>

$$ \left\{ \begin{array}{l} (x_{11} - \lambda)(x_{22} - \lambda) - x_{21}x_{12} = 0 \\ \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_62.png" /></figure>

$$ \left\{ \begin{array}{l} -\lambda x_{11} - \lambda x_{22} + \lambda^2 + x_{11}x_{22} - x_{21}x_{12} = 0 \\ \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_63.png" /></figure>

$$ \left\{ \begin{array}{l} \lambda^2 - 2\lambda \frac{1}{2}(x_{11} + x_{22}) + \left[ \frac{1}{2}(x_{11} + x_{22}) \right]^2 = -(x_{11}x_{22} - x_{21}x_{12}) + \left[ \frac{1}{2}(x_{11} + x_{22}) \right]^2 \\ \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_64.png" /></figure>

$$ \left\{ \begin{array}{l} \left[ \lambda - \frac{1}{2}(x_{11} + x_{22}) \right]^2 = -\text{Det}(A) + \left[ \frac{1}{2}(x_{11} + x_{22}) \right]^2 \\ \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_65.png" /></figure>

$$ \left\{ \begin{array}{l} \lambda - \frac{1}{2}(x_{11} + x_{22}) = \mp \sqrt{-\text{Det}(A) + \left[ \frac{1}{2}(x_{11} + x_{22}) \right]^2} \\ \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_66.png" /></figure>

$$ \left\{ \begin{array}{l} \lambda = \frac{1}{2}(x_{11} + x_{22}) \mp \sqrt{-\text{Det}(A) + \left[ \frac{1}{2}(x_{11} + x_{22}) \right]^2} \\ \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_67.png" /></figure>

$$ \left\{ \begin{array}{l} \lambda = \frac{1}{2}(x_{11} + x_{22}) \mp \sqrt{-x_{11}x_{22} + x_{21}x_{12} + \frac{1}{4}(x_{11}^2 + x_{22}^2 + 2x_{11}x_{22})} \\ \\ \frac{e_2}{e_1} = \frac{x_{21}}{\lambda - x_{22}} \end{array} \right. $$

<figure><img src="image_68.png" /></figure>
$$ (A2-3) \qquad \left\{ \begin{aligned} \lambda &= \frac{(x_{11}+x_{22})\mp\sqrt{x_{11}^2+x_{22}^2+4x_{21}x_{12}-2x_{11}x_{22}}}{2} \\ \frac{e_2}{e_1} &= \frac{x_{21}}{\lambda-x_{22}} \end{aligned} \right. $$

We have two eigenvalues:

$$ (A2-4) \qquad \left\{ \begin{aligned} \lambda_1 &= \frac{x_{11}+x_{22}+\sqrt{x_{11}^2+x_{22}^2+4x_{21}x_{12}-2x_{11}x_{22}}}{2} \\ \lambda_2 &= \frac{x_{11}+x_{22}-\sqrt{x_{11}^2+x_{22}^2+4x_{21}x_{12}-2x_{11}x_{22}}}{2} \end{aligned} \right. $$

Note that a solution exists only if $\lambda \neq x_{ii}$ ($i=1$ or $2$), and if $e_1 \neq 0$ or $e_2 \neq 0$. For a 2x2 matrix, if a solution exists for its characteristic equation, it has two eigenvalues, i.e. the same number as the dimension of the matrix: $\lambda_1$ and $\lambda_2$. For each eigenvalue, we can find an infinite collection of eigenvectors that all satisfy:

$$ \begin{cases} x_{11} \times e_1 + x_{12} \times e_2 = \lambda e_1 \\ x_{21} \times e_1 + x_{22} \times e_2 = \lambda e_2 \end{cases} $$

$$ \begin{cases} (x_{11} - \lambda) \times e_1 = -x_{12} \times e_2 \\ (x_{22} - \lambda) \times e_2 = -x_{21} \times e_1 \end{cases} $$

$$ \begin{cases} e_2 = -\frac{(x_{11} - \lambda)}{x_{12}} \times e_1 \\ (x_{22} - \lambda) \times e_2 = -x_{21} \times e_1 \end{cases} $$

Then, for $e_1=1$, a first pair of eigenvectors would be:

$$ (A2-5) \qquad \left\{ \begin{aligned} \mathbf{e}_1 &= \begin{pmatrix} 1 \\ \frac{\lambda_1 - x_{11}}{x_{12}} \end{pmatrix} \\ \mathbf{e}_2 &= \begin{pmatrix} 1 \\ \frac{\lambda_2 - x_{11}}{x_{12}} \end{pmatrix} \end{aligned} \right. $$

We can go back to EqA2-2 to obtain eigenvalues as function of $x_{21}$ (as is the case in certain textbooks), this leads to:

$$ \left\{ \begin{aligned} \mathbf{e}_1 &= \begin{pmatrix} \frac{\lambda_1 - x_{22}}{x_{21}} \\ 1 \end{pmatrix} \\ \mathbf{e}_2 &= \begin{pmatrix} \frac{\lambda_2 - x_{22}}{x_{21}} \\ 1 \end{pmatrix} \end{aligned} \right. $$

## Appendix 3: Matrix power and diagonalization

Computing matrix powers is difficult, except for diagonal matrices. Indeed, using equation A1-2, it is easy to see that:

$$ \begin{array}{l} (a \quad 0)^2 = (a \quad 0)(a \quad 0) = \begin{pmatrix} a^2 + 0 & 0 + 0 \\ 0 + 0 & 0 + d^2 \end{pmatrix} \\ \Leftrightarrow \\ (A3-1) \qquad (a \quad 0)^t = (a^t \quad 0)^t \end{array} $$
For any other squared matrix A, it may thus be useful to diagonalize it, if one wants to compute any power of it. In Horn and Johnson's book, page 59 (Horn & Johnson, 2013), we are invited to solve the equation:

$$ (A3-2) \qquad \mathbf{P}^{-1} \cdot \mathbf{A} \cdot \mathbf{P} = \mathbf{D} $$

where P is an invertible matrix and D a diagonal matrix.
Let A, P and D, v₁ and v₂ be:

$$ \left\{ \begin{array}{l} \mathbf{A} = \begin{pmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \end{pmatrix} \\ \mathbf{P} = \begin{pmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{pmatrix} \\ \mathbf{D} = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \\ \mathbf{v}_1 = \begin{pmatrix} p_{11} \\ p_{21} \end{pmatrix} \\ \mathbf{v}_2 = \begin{pmatrix} p_{12} \\ p_{22} \end{pmatrix} \end{array} \right. $$

We can also write $\mathbf{P}=(\mathbf{v}_1 \ \mathbf{v}_2)$. We can thus rewrite equation A(3.2) as:

$$ \mathbf{P}^{-1} \cdot \mathbf{A} \cdot (\mathbf{v}_1 \ \mathbf{v}_2) = \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} $$

<figure><img src="image_69.png" /></figure>

$$ \mathbf{P} \cdot \mathbf{P}^{-1} \cdot \mathbf{A} \cdot (\mathbf{v}_1 \ \mathbf{v}_2) = \mathbf{P} \cdot \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} $$

<figure><img src="image_70.png" /></figure>

$$ \mathbf{P} \cdot \mathbf{P}^{-1} \cdot \mathbf{A} \cdot (\mathbf{v}_1 \ \mathbf{v}_2) = \begin{pmatrix} p_{11} & p_{12} \\ p_{21} & p_{22} \end{pmatrix} \cdot \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} $$

<figure><img src="image_71.png" /></figure>

$$ (\mathbf{A} \cdot \mathbf{v}_1 \ \mathbf{A} \cdot \mathbf{v}_2) = \begin{pmatrix} \lambda_1 \cdot p_{11} & \lambda_2 \cdot p_{12} \\ \lambda_1 \cdot p_{21} & \lambda_2 \cdot p_{22} \end{pmatrix} $$

<figure><img src="image_72.png" /></figure>

$$ (\mathbf{A} \cdot \mathbf{v}_1 \ \mathbf{A} \cdot \mathbf{v}_2) = (\lambda_1 \cdot \mathbf{v}_1 \ \lambda_2 \cdot \mathbf{v}_2) $$

<figure><img src="image_73.png" /></figure>

$$ \begin{cases} \mathbf{A} \cdot \mathbf{v}_1 = \lambda_1 \cdot \mathbf{v}_1 \\ \mathbf{A} \cdot \mathbf{v}_2 = \lambda_2 \cdot \mathbf{v}_2 \end{cases} $$

We recognize what we saw about eigenvalues and eigenvectors in Appendix 2, meaning that matrix S is a combination of eigenvectors of A, and D is a diagonal matrix with matrix A's eigenvalues on the diagonal from the bigger (top left) to the smallest (bottom right).

From there, computing the power of any matrix A is relatively easy. Indeed, if we have $\mathbf{P}^{-1} \cdot \mathbf{A} \cdot \mathbf{P}=\mathbf{D}$, then this also writes $\mathbf{P} \cdot \mathbf{P}^{-1} \cdot \mathbf{A} \cdot \mathbf{P}^{-1}=\mathbf{P} \cdot \mathbf{D} \cdot \mathbf{P}^{-1} \Leftrightarrow \mathbf{A}=\mathbf{P} \cdot \mathbf{D} \cdot \mathbf{P}^{-1}$. From there, computing $\mathbf{A}^t$ is easy:

$$ \mathbf{A}^t = (\mathbf{P} \cdot \mathbf{D} \cdot \mathbf{P}^{-1}) \cdot (\mathbf{P} \cdot \mathbf{D} \cdot \mathbf{P}^{-1}) \cdot (\mathbf{P} \cdot \mathbf{D} \cdot \mathbf{P}^{-1}) \cdot (\mathbf{P} \cdot \mathbf{D} \cdot \mathbf{P}^{-1}) \cdots (\mathbf{P} \cdot \mathbf{D} \cdot \mathbf{P}^{-1}) $$

<figure><img src="image_74.png" /></figure>

$$ \mathbf{A}^t = \mathbf{P} \cdot \mathbf{D} \cdot (\mathbf{P}^{-1} \cdot \mathbf{P}) \cdot \mathbf{D} \cdot (\mathbf{P}^{-1} \cdot \mathbf{P}) \cdot \mathbf{D} \cdot (\mathbf{P}^{-1} \cdot \mathbf{P}) \cdots (\mathbf{P}^{-1} \cdot \mathbf{P}) \cdot \mathbf{D} \cdot \mathbf{P}^{-1} $$

<figure><img src="image_75.png" /></figure>

$$ \mathbf{A}^t = \mathbf{P} \cdot \mathbf{D} \cdot \mathbf{I} \cdot \mathbf{D} \cdot \mathbf{I} \cdot \mathbf{D} \cdot \mathbf{I} \cdots \mathbf{I} \cdot \mathbf{D} \cdot \mathbf{P}^{-1} $$

where I is the identity matrix. From there, we can compute:

$$ (A3-3) \qquad \mathbf{A}^t = \mathbf{P} \cdot \mathbf{D}^t \cdot \mathbf{P}^{-1} $$

Consequently, we can use equation (A3-3) to calculate the power of any diagonalizable square matrix.
We can now derive some other properties of eigenvalue-eigenvector pairs (eigenpairs).

For the 2x2 matrix A, with the eigenpairs $\lambda_i$ and $\mathbf{e}_i$, where $i$ stands for 1 or 2, $A.e_i=\lambda_i.e_i$. Then $A^2.e_i=A.(A.e_i)=A.(\lambda_i.e_i)=\lambda_i.A.e_i=\lambda_i^2.e_i$. It follows that:

$$ (A3-4) \qquad A^t.e_i = \lambda_i^t.e_i $$

Let $\mathbf{v}$ be a vector composed of a combination of eigenvectors of matrix $A$ so that $\mathbf{v}=\sum_i x_i.e_i$, where the $x_i$'s are scalars that can be computed. We can then write:

$$ \begin{align*} A.\mathbf{v} &= A.\sum_i x_i.e_i \\ A.\mathbf{v} &= \sum_i x_i.A.e_i = \sum_i x_i.\lambda_i.e_i \\ A.A.\mathbf{v} &= A.\sum_i x_i.\lambda_i.e_i \\ A^2.\mathbf{v} &= \sum_i x_i.\lambda_i^2.e_i \end{align*} $$
$$ (A3-5) \qquad A^t.\mathbf{v} = \sum_i x_i.\lambda_i^t.e_i $$

This property can be used to any power function of matrices. In particular, for the matrix $S=(I-\gamma A)$, which should be invertible, it is easy to see that the eigenpairs of $S$ are (in decreasing order of the hierarchy) $\lambda_1'=1-\gamma\lambda_2$, $\mathbf{e}_1'=\mathbf{e}_2$ and $\lambda_2'=1-\gamma\lambda_1$, $\mathbf{e}_2'=\mathbf{e}_1$. Indeed, if we take the eigenvector $\mathbf{e}_i$, then: $S.e_i=\mathbf{e}_i-\gamma A.e_i \Leftrightarrow S.e_i=\mathbf{e}_i-\gamma\lambda_i.e_i \Leftrightarrow S.e_i = (I-\gamma\lambda_i).\mathbf{e}_i$ (QED).

Using equation A3-5, we can write:

$$ \begin{align*} S^t.\mathbf{v} &= \sum_i x_i (1 - \gamma \lambda_i)^t \cdot \mathbf{e}_i \\ \text{This logically yields:} \\ (I - \gamma A)^{-1} \cdot \mathbf{v} &= \sum_i x_i (1 - \gamma \lambda_i)^{-1} \cdot \mathbf{e}_i \end{align*} $$
$$ (A3-6) \qquad (I - \gamma A)^{-1} \cdot \mathbf{v} = \sum_i \frac{1}{1 - \gamma \lambda_i} x_i \cdot \mathbf{e}_i $$

Equation A3-6 corresponds to equation A.10 of Rousset's book (page 219), which was given without any detailed proof.

## Appendix 4: eigenvalue effective population size

This notion refers to the evolution of heterozygosity, or more exactly genetic diversity, defined as the probability, at generation $t$, to draw randomly two alleles that are not identical by descent, from one generation to the other, and labelled $H_{t-1}$ and $H_t$. If the evolution of a population by genetic drift has reached a steady state, the ratio of $H_t/H_{t-1}$ remains constant generation after generation and has been shown to correspond to the leading eigenvalue of the transition matrix for the evolution of vectors of genetic identity probabilities (see below).

Let $Q_i$ be the probability of identity within diploid individuals, and $Q_S$, the probability of identity between two alleles from two individuals of the same population. In an ideal population of size $N$, and thus under panmixia (thus here $Q_i=Q_S$), we can set the system of equations:
<figure><img src="image_76.png" /></figure>

$$ 
\left\{ 
\begin{aligned}
Q_{I_t} &= \frac{1}{N} \left( \frac{1}{2} + \frac{1}{2} Q_{I_{t-1}} \right) + \left( 1 - \frac{1}{N} \right) Q_{S_{t-1}} \\
Q_{S_t} &= \frac{1}{N} \left( \frac{1}{2} + \frac{1}{2} Q_{I_{t-1}} \right) + \left( 1 - \frac{1}{N} \right) Q_{S_{t-1}}
\end{aligned}
\right. 
$$

<figure><img src="image_77.png" /></figure>

$$ 
\left\{ 
\begin{aligned}
Q_{I_t} &= \frac{1}{2N} Q_{S_{t-1}} + \left( 1 - \frac{1}{N} \right) Q_{S_{t-1}} + \frac{1}{2N} \\
Q_{S_t} &= \frac{1}{2N} Q_{S_{t-1}} + \left( 1 - \frac{1}{N} \right) Q_{S_{t-1}} + \frac{1}{2N}
\end{aligned}
\right. 
$$

$$ 
\left\{ 
\begin{aligned}
Q_{I_t} &= \left( 1 - \frac{1}{2N} \right) Q_{S_{t-1}} + \frac{1}{2N} \\
Q_{S_t} &= \left( 1 - \frac{1}{2N} \right) Q_{S_{t-1}} + \frac{1}{2N}
\end{aligned}
\right. 
$$

If we replace identities with their corresponding values in terms of genetic non-identity (hence diversity), we obtain:

<figure><img src="image_78.png" /></figure>

$$ 
\left\{ 
\begin{aligned}
1 - H_{I_t} &= \left( 1 - \frac{1}{2N} \right) \left( 1 - H_{S_{t-1}} \right) + \frac{1}{2N} \\
1 - H_{S_t} &= \left( 1 - \frac{1}{2N} \right) \left( 1 - H_{S_{t-1}} \right) + \frac{1}{2N}
\end{aligned}
\right. 
$$

(A4-1)

$$ 
\left\{ 
\begin{aligned}
H_{I_t} &= \left( 1 - \frac{1}{2N} \right) H_{S_{t-1}} \\
H_{S_t} &= \left( 1 - \frac{1}{2N} \right) H_{S_{t-1}}
\end{aligned}
\right. 
$$

Assuming a steady state, so that $H_{S_t}/H_{S_{t-1}}=H_{S_{t-1}}/H_{S_{t-2}}=\lambda$, we can set:

$$ \lambda = 1 - \frac{1}{2N} $$

Let us now define the vector $\mathbf{H}_t$ and transition matrix $\mathbf{A}$, as:

$$ 
\mathbf{H}_t = \begin{pmatrix} H_{I_t} \\ H_{S_t} \end{pmatrix} \\
\mathbf{A} = \begin{pmatrix} 0 & 1 - \frac{1}{2N} \\ 0 & 1 - \frac{1}{2N} \end{pmatrix}
$$

Using these, equation A1-1 also writes: $\mathbf{H}_t=\mathbf{A}\mathbf{H}_{t-1} \Leftrightarrow \mathbf{H}_t=\mathbf{A}^2\mathbf{H}_{t-2} \Leftrightarrow \mathbf{H}_t=\mathbf{A}^t\mathbf{H}_0$. This also writes:

(A4-2)

$$ 
\begin{pmatrix} H_{I_t} \\ H_{S_t} \end{pmatrix} = \begin{pmatrix} 0 & 1 - \frac{1}{2N} \\ 0 & 1 - \frac{1}{2N} \end{pmatrix}^t \cdot \begin{pmatrix} H_1 \\ H_S \end{pmatrix}
$$

were $H_1$ and $H_S$ are genetic diversities at time 0.

We can decompose $\mathbf{H}_0$ as a function of eigenvectors of $\mathbf{A}$ ($\mathbf{e}_1$ and $\mathbf{e}_2$) (see appendix 3) and scalars $c_1$ and $c_2$ such as:

(A4-3)

$$ \begin{pmatrix} H_1 \\ H_S \end{pmatrix} = c_1 \cdot \mathbf{e}_1 + c_2 \cdot \mathbf{e}_2 $$
Combining Equations A4-2 and A4-3 yields:

$$ \binom{H_{1t}}{H_{St}} = \begin{pmatrix} 0 & 1 - \frac{1}{2N} \\ 0 & 1 - \frac{1}{2N} \end{pmatrix}^t \cdot (c_1 \cdot \mathbf{e}_1 + c_2 \cdot \mathbf{e}_2) $$

<figure><img src="image_79.png" /></figure>

$$ (A4-4) \qquad \binom{H_{1t}}{H_{St}} = c_1 \cdot \begin{pmatrix} 0 & 1 - \frac{1}{2N} \\ 0 & 1 - \frac{1}{2N} \end{pmatrix}^t \cdot \mathbf{e}_1 + c_2 \cdot \begin{pmatrix} 0 & 1 - \frac{1}{2N} \\ 0 & 1 - \frac{1}{2N} \end{pmatrix}^t \cdot \mathbf{e}_2 $$

Following what we know from the properties of eigenpairs (see Appendices 2 and 3), using equation A3-5, we can rewrite equation A4-4 as:

$$ (A4-5) \qquad \binom{H_{1t}}{H_{St}} = c_1 \cdot \lambda_1^t \cdot \mathbf{e}_1 + c_2 \cdot \lambda_2^t \cdot \mathbf{e}_2 $$

With eigenpairs $\lambda_i$ and $\mathbf{e}_i$ of matrix $\mathbf{A}=\begin{pmatrix} 0 & 1 - \frac{1}{2N} \\ 0 & 1 - \frac{1}{2N} \end{pmatrix}$, using EqA2-4, we can compute the two eigenvalues of matrix $\mathbf{A}$:

$$ \left\{ \begin{aligned} \lambda_1 &= \frac{0 + 1 - \frac{1}{2N} + \sqrt{0^2 + \left(1 - \frac{1}{2N}\right)^2 + 4 \times 0 \times \left(1 - \frac{1}{2N}\right) - 2 \times 0 \times \left(1 - \frac{1}{2N}\right)}}{2} \\ \lambda_2 &= \frac{0 + 1 - \frac{1}{2N} - \sqrt{0^2 + \left(1 - \frac{1}{2N}\right)^2 + 4 \times 0 \times \left(1 - \frac{1}{2N}\right) - 2 \times 0 \times \left(1 - \frac{1}{2N}\right)}}{2} \end{aligned} \right. $$

<figure><img src="image_80.png" /></figure>

$$ \left\{ \begin{aligned} \lambda_1 &= \frac{1 - \frac{1}{2N} + \sqrt{\left(1 - \frac{1}{2N}\right)^2}}{2} \\ \lambda_2 &= \frac{1 - \frac{1}{2N} - \sqrt{\left(1 - \frac{1}{2N}\right)^2}}{2} \end{aligned} \right. $$

<figure><img src="image_81.png" /></figure>

$$ (A4-6) \qquad \left\{ \begin{aligned} \lambda_1 &= 1 - \frac{1}{2N} \\ \lambda_2 &= 0 \end{aligned} \right. $$

For the eigenvectors of $\mathbf{A}$, using equation A2-5, we have:

$$ \left\{ \begin{aligned} \mathbf{e}_1 &= \begin{pmatrix} \frac{1}{1 - \frac{1}{2N} - 0} \\ \frac{1 - \frac{1}{2N} - 0}{1 - \frac{1}{2N}} \end{pmatrix} \\ \mathbf{e}_2 &= \begin{pmatrix} \frac{1}{0 - 0} \\ \frac{1 - \frac{1}{2N}}{1 - \frac{1}{2N}} \end{pmatrix} \end{aligned} \right. $$
<figure><img src="image_82.png" /></figure>

$$ (A4-7) \qquad \begin{cases} \mathbf{e}_1 = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \\ \mathbf{e}_2 = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \end{cases} $$

Combining A4-7 and A4-6 with A1-5 we can write:

$$ \binom{H_{I_t}}{H_{S_t}} = c_1 \cdot \left(1 - \frac{1}{2N}\right)^t \cdot \binom{1}{1} $$

<figure><img src="image_83.png" /></figure>

$$ (A4-8) \qquad \begin{cases} H_{I_t} = c_1 \cdot \left(1 - \frac{1}{2N}\right)^t \\ H_{S_t} = c_1 \cdot \left(1 - \frac{1}{2N}\right)^t \end{cases} $$

From there, we can also easily see that, for any $c_1 \neq 0$:

$$ (A4-9) \qquad \frac{H_{S_t}}{H_{S_{t-1}}} = \frac{c_1 \cdot \left(1 - \frac{1}{2N}\right)^t}{c_1 \cdot \left(1 - \frac{1}{2N}\right)^{t-1}} = \left(1 - \frac{1}{2N}\right) = \lambda_1 $$

The ratio of genetic diversities of generation *t* and *t-1* is indeed the leading eigenvalue of the transition matrix describing the evolution of genetic diversities (and of genetic identities as well) (QED).

We can also determine $c_1$ and $c_2$, if genetic diversities at time 0 are known. From equations A4-3 and A4-7, we know that:

$$ \begin{cases} H_I = c_1 + c_2 \\ H_S = c_1 \end{cases} $$

<figure><img src="image_84.png" /></figure>

$$ (A4-10) \qquad \begin{cases} c_1 = H_S \\ c_2 = H_I - H_S \end{cases} $$

Now if we combine equations A1-8 and A1-10 we can compute, for $H_S$ (which is here the same as $H_I$):

$$ \begin{cases} H_{I_t} = H_S \cdot \left(1 - \frac{1}{2N}\right)^t \\ H_{S_t} = H_S \cdot \left(1 - \frac{1}{2N}\right)^t \end{cases} $$

This results confirms that, at any generation $H_I=H_S$, and hence $c_2=0$. We can then simply write, for the Wright-Fisher model:

$$ H_{I_t} = H_S \cdot \left(1 - \frac{1}{2N}\right)^t $$

<figure><img src="image_85.png" /></figure>

$$ (A4-11) \qquad H_{S_t} = H_S \cdot \lambda_1^t $$

where $H_S$ is the local genetic diversity at time 0, and $\lambda_1$ is the leading eigenvalue of the transition matrix A.
Now, for any transition matrix $\mathbf{A} = \begin{pmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \end{pmatrix}$, where the $x_{ij}$ are probabilities (e.g. of identity), with eigenvalues $\lambda_1$ and $\lambda_2$ and corresponding eigenvectors $\mathbf{e}_1$ and $\mathbf{e}_2$ (Appendix 2), we can use the same approach and obtain:

$$ (A4-12) \qquad \begin{pmatrix} H_{I_t} \\ H_{S_t} \end{pmatrix} = c_1 \cdot \lambda_1^t \cdot \mathbf{e}_1 + c_2 \cdot \lambda_2^t \cdot \mathbf{e}_2 $$

From equation A2-4, it is very easy to see that $\lambda_1+\lambda_2=x_{11}+x_{22}$, and that $\lambda_1>\lambda_2$. From matrix $\mathbf{A}$ defined for the WF model, we can see that all $x_{ij}$ are probabilities that only sum to 1 for the Castle-Weinberg model and hence $x_{11}+x_{22}\le 1$. In the case of WF, the difference between $\lambda_1$ and $\lambda_2$ is even very big (equation A4-6). Consequently, when $t$ becomes big enough $\lambda_1^t>>\lambda_2^t$ and equation A1-12 can be approximated as:

$$ (A4-13) \qquad \begin{pmatrix} H_{I_t} \\ H_{S_t} \end{pmatrix} \approx c_1 \cdot \lambda_1^t \cdot \mathbf{e}_1 $$

Combining A4-13 with A2-5 yields:

$$ (A4-14) \qquad \left\{ \begin{array}{l} H_{I_t} \approx c_1 \cdot \lambda_1^t \\ H_{S_t} \approx c_1 \cdot \lambda_1^t \frac{\lambda_1 - x_{11}}{x_{12}} \end{array} \right. $$

From there, and for any $H_x$ (X=I or S), it is straightforward that the ratio $H_t/H_{t-1}=\lambda_1$. From equation A4-6, we know that the leading eigenvalue of the transition matrix of a fully panmictic model (i.e. WF) is $\lambda_e=1-1/(2N_e)$ and thus, the effective population size of any non-reference population will follow $\lambda_1=1-1/(2N_e)$, which is equivalent to:

$$ 1 - \lambda_1 \approx \frac{1}{2N_e} $$

Consequently, the eigenvalue effective population size of any population will be:

$$ (A4-15) \qquad N_e \approx \frac{1}{2(1-\lambda_1)} $$

where $\lambda_1$ is the leading eigenvalue of the transition matrix, describing the evolution of genetic diversities (or same wise of genetic identities) from one generation to the other, for that population.

All these detailed explanations leading to equation A4-15 provide the same result as equation 3.105 in Ewen's book (page 120) (Ewens, 2004), which was given with much more elliptic explanations.

It is also worth noting that equation A4-15 is only accurate when $t$ is big enough, or when the population has reached a steady state so that the ratio $H_t/H_{t-1}$ becomes constant and equal to $\lambda_1$.

## Appendix 5: Pudovkin et al.'s methods to compute $N_e$

Let $p_f$ and $p_m$ be allele frequencies of one of two alleles at a given locus in females and males respectively, in a population with an even sex-ratio. Then, in the progeny, the proportion of heterozygotes observed should be $H_{exp-dio}=p_f(1-p_m)+(1-p_f)p_m$. In this population, the frequency of this allele will be $(p_f+p_m)/2$. Consequently, the expected frequency of heterozygotes under the panmictic (monoecious) model in the progeny ($H_{expmon}$) would be:

$$ H_{exp-mon} = 2 \left( \frac{p_f + p_m}{2} \right) \left( 1 - \frac{p_f + p_m}{2} \right) \Leftrightarrow H_{exp-mon} = (p_f + p_m) \left( \frac{2 - p_f - p_m}{2} \right) $$

$$ H_{exp-mon} = (p_f + p_m) \left( \frac{2 - p_f - p_m}{2} \right) \Leftrightarrow H_{exp-mon} = (p_f + p_m) \left( \frac{2 - p_f - p_m}{2} \right) $$
$$H_{\text{exp-mon}} = \frac{2p_f + 2p_m - p_f^2 - p_m^2 - 2p_f p_m}{2}$$

<figure><img src="image_86.png" /></figure>

$$H_{\text{exp-on}} = p_f + p_m - p_f p_m - \frac{p_f^2 + p_m^2}{2}$$

<figure><img src="image_87.png" /></figure>

$$H_{\text{exp-mon}} = p_f(1-p_m) + p_m(1-p_f) - \frac{p_f^2 + p_m^2 - 2p_f p_m}{2}$$

<figure><img src="image_88.png" /></figure>

$$H_{\text{exp-mon}} = H_{\text{exp-dio}} - \frac{(p_f - p_m)^2}{2}$$

<figure><img src="image_89.png" /></figure>

$$(A5-1) \qquad H_{\text{exp-dio}} = H_{\text{exp-mon}} + \frac{(p_f - p_m)^2}{2}$$

Please note that $H_{\text{exp-dio}}$ and $H_{\text{exp-mon}}$ here correspond to $H_{\text{obs}}$ and $H_{\text{exp}}$ respectively in (Pudovkin et al., 1996). There is thus an observed heterozygote excess in the progeny.

The quantity $p_f p_m$ can be considered as a random variable with average 0 over all possible parental groups. If we consider that the frequency of the first allele was $p$ in the parental population, then the average of $(p_f p_m)^2$ is the variance of a difference in allele frequencies between two binomial samples of size $N$ for each gender ($N$ alleles in females and $N$ in males=2$N$ alleles in total). The variance of frequency of a given allele randomly taken in a population of size $N$, is $p(1-p)/N$, where $p$ is the frequency of the allele in the parents. The variance of a difference between two uncorrelated (e.g. independent) random variables is the sum of individual variances (see the file VarDif.pdf), here $p(1-p)/(N)+p(1-p)/(N)=2p(1-p)/N$. Equation (A5-1), for the entire space of possible outcomes can thus write:

$$(A5-2) \qquad H_{\text{exp-dio}} = H_{\text{exp-mon}} + \frac{p(1-p)}{N}$$

If we replace $p(1-p)$ as $H_{\text{exp-mon}}/2$, $N$ by $N_e$ and rearrange equation (A5-2), one obtains:

$$H_{\text{exp-dio}} = H_{\text{exp-mon}} + \frac{\frac{H_{\text{exp-mon}}}{2}}{N_e} - \frac{H_{\text{exp-mon}}}{2N_e} = H_{\text{exp-mon}} - H_{\text{exp-dio}}$$

<figure><img src="image_90.png" /></figure>

$$(A5-3) \qquad N_e = -\frac{\frac{H_{\text{exp-mon}}}{2}}{2H_{\text{exp-mon}}-H_{\text{obsexp-dio}}}$$

The parametric formula of Wright's $F_{IS}$ can be written as (Nei & Chesser, 1983):

$$(A5-4) \qquad F_{IS} = \frac{H_{\text{exp}} - H_{\text{obs}}}{H_{\text{exp}}}$$

If we combine equations (A5-3) and (A5-4), replacing $H_{\text{exp}}$ with $H_{\text{exp-mon}}$ and $H_{\text{obs}}$ with $H_{\text{exp-dio}}$, we obtain the same result (with $F_{IS}$) as equation (3) in Pudovkin et al.'s paper (if we replace $F_{IS}$ by $-D$):

$$(A5-5) \qquad N_e = -\frac{1}{2F_{IS}}$$

In their appendix, Pudovkin et al. then used a sleight of hand. They set $N=N_e$ again, $2p(1-p)=H_{t-1}$ and $H_{\text{exp-mon}}=H_t$, and used the equation $\lambda=H_t/H_{t-1}$, citing Kimura and Crow's book (Crow & Kimura, 1970). Then, with $p(1-p)=H_{t-1}/2$, $H_{t-1}=H_t/\lambda$ and $H_{\text{exp}}=H_t$, we can rewrite A5-2:
$$H_{\text{exp-dio}} = H_t + \frac{H_{t-1}}{2N_e}$$

<figure><img src="image_91.png" /></figure>

$$H_{\text{exp-dio}} = H_t + H_t \frac{1}{2\lambda N_e}$$

<figure><img src="image_92.png" /></figure>

$$(A5-6) \qquad H_{\text{exp-dio}} = H_t \left( 1 + \frac{1}{2\lambda N_e} \right)$$

Pudovkin et al. used another sleight of hand from equation 3.11.8 (page 104) from Crow and Kimura's book (Crow & Kimura, 1970), replacing subpopulation sizes N by N<sub>e</sub> (again) and obtained:

$$(A5-7) \qquad \lambda = \frac{N_e - 1 + \sqrt{N_e^2 + 1}}{2N_e}$$

The way Pudovkin et al used this equation may be inaccurate because Crow and Kimura's equation refers to the number of individuals, not the effective population size. Nevertheless, if we combine equations A5-6 and A5-7 we obtain:

$$H_{\text{exp-dio}} = H_t \left( 1 + \frac{1}{2 \frac{N_e - 1 + \sqrt{N_e^2 + 1}}{2N_e} N_e} \right)$$

<figure><img src="image_93.png" /></figure>

$$H_{\text{exp-dio}} = H_t \left( 1 + \frac{1}{N_e - 1 + \sqrt{N_e^2 + 1}} \right)$$

<figure><img src="image_94.png" /></figure>

$$H_{\text{exp-dio}} - H_t = H_t \left( 1 + \frac{1}{N_e - 1 + \sqrt{N_e^2 + 1}} - 1 \right)$$

<figure><img src="image_95.png" /></figure>

$$(A5-8) \qquad \frac{H_{\text{exp-dio}} - H_t}{H_t} = \frac{1}{2N_e - 1 + \sqrt{N_e^2 + 1}}$$

From equation A5-4, and setting that H<sub>t</sub> is the expected heterozygote frequency in the progeny, hence H<sub>exp-mon</sub>, we can rewrite equation A5-8 as:

$$-F_{\text{IS}} = \frac{1}{N_e - 1 + \sqrt{N_e^2 + 1}}$$

<figure><img src="image_96.png" /></figure>

$$-F_{\text{IS}}N_e + F_{\text{IS}} - F_{\text{IS}}\sqrt{N_e^2 + 1} = 1$$

<figure><img src="image_97.png" /></figure>

$$-F_{\text{IS}}N_e + F_{\text{IS}} - 1 = F_{\text{IS}}\sqrt{N_e^2 + 1}$$

<figure><img src="image_98.png" /></figure>

$$(-F_{\text{IS}}N_e + F_{\text{IS}} - 1)(-F_{\text{IS}}N_e + F_{\text{IS}} - 1) = F_{\text{IS}}^2(N_e^2 + 1)$$
$$ \Leftrightarrow F_{IS}^2 N_e^2 - F_{IS}^2 N_e + F_{IS} N_e - F_{IS} N_e + F_{IS}^2 - F_{IS} + F_{IS} N_e - F_{IS} + 1 - F_{IS}^2 N_e^2 - F_{IS}^2 = 0 $$

$$ \Leftrightarrow -2F_{IS}^2 N_e + 2F_{IS} N_e - 2F_{IS} + 1 = 0 $$

$$ \Leftrightarrow 2F_{IS}N_e(1 - F_{IS}) = 2F_{IS} - 1 $$

$$ \Leftrightarrow N_e = -\frac{1 - 2F_{IS}}{2F_{IS}(1 - F_{IS})} $$

$$ \Leftrightarrow N_e = -\frac{1 - F_{IS} - F_{IS}}{2F_{IS}(1 - F_{IS})} $$

$$ \Leftrightarrow N_e = -\frac{1 - F_{IS}}{2F_{IS}(1 - F_{IS})} + \frac{F_{IS}}{2F_{IS}(1 - F_{IS})} $$

$$ \Leftrightarrow \quad (A5-9) \qquad N_e = -\frac{1}{2F_{IS}} + \frac{1}{2(1-F_{IS})} $$

Considering that $F_{IS}=-D$, we get:

$$ \Leftrightarrow \quad (A5-10) \qquad N_e = \frac{1}{2D} + \frac{1}{2(D+1)} $$

Equation A5-10 is the same as Pudovkin et al. (1996) equation 4.

## Appendix 6: Coalescent effective population size in a dioecious pangamic population

Let $Q_I$ and $Q_S$ be the probabilities that the same allele is sampled twice, either in one individual or in two distinct individuals from the same population. Let $u$ be the mutation rate per generation in an infinite allele model where each mutation event produces a new allele that never existed before (no homoplasy). Then, for a dioecious population with an even sex-ratio and random mating, we can set the following recurrences between generation $t$ and $t-1$ (Equations 7 and 8 with an even sex-ratio and mutation rate $u$) (see equations 7 and 8 with $N_f=N_m$).

$$ \Leftrightarrow \quad (A6-1) \qquad \left\{ \begin{array}{l} Q_{I(t)} = (1-u)^2 Q_{S(t-1)} \\ Q_{S(t)} = (1-u)^2 \left[ \frac{1}{N} \left( \frac{1}{2} + \frac{1}{2} Q_{I(t-1)} \right) + \left( 1 - \frac{1}{N} \right) Q_{S(t-1)} \right] \end{array} \right. $$

$$ \Leftrightarrow \quad (A6-1) \qquad \left\{ \begin{array}{l} Q_{I(t)} = (1-u)^2 Q_{S(t-1)} \\ Q_{S(t)} = (1-u)^2 \frac{1}{2N} Q_{I(t-1)} + (1-u)^2 \left(1 - \frac{1}{N}\right) Q_{S(t-1)} + (1-u)^2 \frac{1}{2N} \end{array} \right. $$

Let $Q_t$ the vector of genetic identities at time $t$ and $A$ be the squared matrix of transition for genetic identities, $v$ the corresponding vector of residuals, and $I$ the identity matrix. If $\gamma=(1-u)^2$ is the probability that two alleles taken at random did not mutate, then we can write:

$$ \Leftrightarrow \quad (A6-2) \qquad Q_t = \gamma A Q_{t-1} + \gamma v $$

For the example of a dioecious population with even sex ratio this would yield (see equation A6-1):
$$ (A6-3) \qquad \left\{ \begin{array}{l} \mathbf{Q}_t = \begin{pmatrix} Q_{I(t)} \\ Q_{S(t)} \end{pmatrix} \\[1em] \mathbf{A} = \begin{pmatrix} 0 & 1 \\ \frac{1}{2N} & \left(1 - \frac{1}{N}\right) \end{pmatrix} \\[1em] \mathbf{v} = \begin{pmatrix} 0 \\ \frac{1}{2N} \end{pmatrix} \end{array} \right. $$

Equation A6-2 is equivalent to:

$$ \Leftrightarrow \qquad \mathbf{Q}_t = \gamma \mathbf{A} (\gamma \mathbf{A} \mathbf{Q}_{t-2} + \gamma \mathbf{v}) + \gamma \mathbf{v} $$

$$ \Leftrightarrow \qquad \mathbf{Q}_t = \gamma^2 \mathbf{A}^2 \mathbf{Q}_{t-2} + \gamma^2 \mathbf{A} \mathbf{v} + \gamma \mathbf{v} $$

$$ \Leftrightarrow \qquad \mathbf{Q}_t = \gamma^2 \mathbf{A}^2 (\gamma \mathbf{A} \mathbf{Q}_{t-3} + \gamma \mathbf{v}) + \gamma^2 \mathbf{A} \mathbf{v} + \gamma \mathbf{v} $$

$$ \Leftrightarrow \qquad \mathbf{Q}_t = \gamma^3 \mathbf{A}^3 \mathbf{Q}_{t-3} + \gamma^3 \mathbf{A}^2 \mathbf{v} + \gamma^2 \mathbf{A} \mathbf{v} + \gamma \mathbf{v} $$

$$ \mathbf{Q}_t = \gamma^{(t-1)} \mathbf{A}^{(t-1)} \mathbf{Q}_1 + \gamma^{(t-1)} \mathbf{A}^{(t-2)} \mathbf{v} + \gamma^{(t-2)} \mathbf{A}^{(t-3)} \mathbf{v} + \dots + \gamma^2 \mathbf{A}^1 \mathbf{v} + \gamma^1 \mathbf{A}^0 \mathbf{v} $$

Assuming that equilibrium values has been reached at time $t$ ($t \to \infty$):

$$ (A6-4) \qquad \mathbf{Q} = \gamma^{(t-1)} \mathbf{A}^{(t-1)} \mathbf{Q}_1 + \left( \sum_{i=1}^{\infty} \gamma^t \mathbf{A}^{(t-1)} \mathbf{v}_i \right) $$

(A6-4)

We can see that the second term in these equations will increase with $t$, albeit at a diminishing rate, while the first term will decrease with $t$. Hence, if inbreeding within individuals and within subpopulation are small enough at time $t=1$, after a sufficient number of generations, and using equation A3-5 and decomposing $\mathbf{v}$ as $\mathbf{v}=\sum_i x_i \cdot \mathbf{e}_i$, where the $x_i$'s are scalars that can be computed and $\mathbf{e}_i$ are eigenvectors of $\mathbf{A}$, we can approximate equation A6-4 as:

$$ (A6-5) \qquad \mathbf{Q} \approx \sum_{t=1}^{\infty} \gamma^t \sum_i \lambda_i^{(t-1)} x_i \mathbf{e}_i $$

This is the same as the second part of equation 4.10 in Rousset's book, page 56 (Rousset, 2004). It is worthy of note that such an approximation is invalid in populations with poor levels of genetic diversity in the first generations.

We can also compute $\mathbf{Q}$ at equilibrium. For this we set equation A6-2 as:

$$ \Leftrightarrow \qquad \mathbf{Q} = \gamma \mathbf{A} \mathbf{Q} + \gamma \mathbf{v} $$

$$ \Leftrightarrow \qquad \mathbf{Q}(\mathbf{I} - \gamma \mathbf{A}) = \gamma \mathbf{v} $$

$$ (A6-6) \qquad \mathbf{Q} = \gamma (\mathbf{I} - \gamma \mathbf{A})^{-1} \mathbf{v} $$

We can use equation A3-6 to obtain:

$$ (A6-7) \qquad \mathbf{Q} = \gamma \sum_i \frac{1}{1 - \gamma \lambda_i} \cdot x_i \cdot \mathbf{e}_i $$

This equation corresponds to the first part of equation 4.10 given in Rousset's book.

We can also express $\mathbf{Q}$ as a function of probability of pairwise coalescence at time $t$. If we define a vector $\mathbf{C}_t$ of such probabilities within individuals and between individuals (to stick to our framework with two hierarchies), we can write:

$$ (A6-8) \qquad \mathbf{C}_t = \begin{pmatrix} c_{I(t)} \\ c_{S(t)} \end{pmatrix} $$
In equation A6-8, C(t) and C(S,t) are the probabilities, at time t, that two alleles of one individual (I) or of different individuals in the subpopulation (S), respectively, all randomly chosen, had coalesced somewhere in the past. At equilibrium, or after a lot of generations, identities will correspond to the sum of all coalescent events that occurred in the past, and if no mutation ever occurred, and hence:

$$ (A6-9) \qquad Q = \sum_{t=1}^{\infty} C_t \gamma^t $$

Combining equations A6-5 and A6-9 provides the following equality:

$$ \sum_{t=1}^{\infty} \gamma^t \sum_{i} \lambda_i^{(t-1)} x_i e_i \approx \sum_{t=1}^{\infty} \gamma^t C_t $$

This means that:

$$ (A6-10) \qquad C_t \approx \sum_i \lambda_i^{(t-1)} x_i e_i $$

This equation meets with equation 4.11 page 56 in Rousset's book (Rousset, 2004).

The mean coalescent time between two alleles in hierarchy J T(t) (J=I or S for the example treated in the present paper) at time t, can be computed as the sum of the products of the time of each event of coalescence by the probability of coalescence at that time for these two alleles of J (Rousset, 2004) (page 59), in vector format:

$$ (A6-11) \qquad T_n = \sum_{t=1}^{n} t C_t $$

Please note that in Rousset's book or other papers n=∞.

Using the result of equation A6-10 we get:

$$ (A6-12) \qquad T_n = \sum_{t=1}^{n} t \sum_{i} \lambda_i^{(t-1)} x_i e_i $$

<figure><img src="image_99.png" /></figure>

$$ T_n = \sum_{t=1}^{n} t \lambda_1^{(t-1)} x_1 e_1 + \sum_{t=1}^{n} t \lambda_2^{(t-1)} x_2 e_2 + \sum_{t=1}^{n} t \lambda_3^{(t-1)} x_3 e_3 + \dots + \sum_{t=1}^{n} t \lambda_i^{(t-1)} x_i e_i $$

for the general case of any squared transition matrices (for the present case this sum stops at λ₂).

The eigenpairs and scalars are constant through time and we can focus on the different sums, Sᵢ of each eigenpair of order i:

$$ S_i = \sum_{t=1}^{n} t \lambda_i^{(t-1)} = 1\lambda_i^0 + 2\lambda_i^1 + 3\lambda_i^2 + \dots + n\lambda_i^{(n-1)} $$

<figure><img src="image_100.png" /></figure>

$$ \lambda_i S_i = \lambda_i^1 + 2\lambda_i^2 + \dots + (t-1)\lambda_i^{t-1} + n\lambda_i^n $$

We can then set:

$$ S_i - \lambda_i S_i = 1 + \lambda_i^1 + \lambda_i^2 + \dots + \lambda_i^{(n-1)} - n\lambda_i^n $$

<figure><img src="image_101.png" /></figure>

$$ S_i(1 - \lambda_i) = 1 + \lambda_i^1 + \lambda_i^2 + \dots + \lambda_i^{(n-1)} - n\lambda_i^n $$

<figure><img src="image_102.png" /></figure>

$$ S_i(1 - \lambda_i) = S_i' - n\lambda_i^n $$

<figure><img src="image_103.png" /></figure>

$$ S_i = \frac{S_i' - n\lambda_i^n}{1 - \lambda_i} $$

where

$$ S_i' = 1 + \lambda_i^1 + \lambda_i^2 + \dots + \lambda_i^{(n-1)} $$

<figure><img src="image_104.png" /></figure>

$$ \lambda_i S_i' = \lambda_i^1 + \lambda_i^2 + \dots + \lambda_i^{(n-1)} + \lambda_i^n $$
We can again use the fact that:

<figure><img src="image_105.png" /></figure>

$$S'_{i} - \lambda_{i} S'_{i} = 1 - \lambda_{i}^{n}$$

<figure><img src="image_106.png" /></figure>

$$S'_{i}(1 - \lambda_{i}) = 1 - \lambda_{i}^{n}$$

$$S'_{i} = \frac{1 - \lambda_{i}^{n}}{1 - \lambda_{i}}$$

We can now replace this $S'_i$ in $S_i$ to obtain:

<figure><img src="image_107.png" /></figure>

$$S_{i} = \frac{\frac{1 - \lambda_{i}^{n}}{1 - \lambda_{i}} - t\lambda_{i}^{n}}{1 - \lambda_{i}}$$

$$S_i = \sum_{t=1}^{n} t \lambda_i^{(t-1)} = \frac{1 - \lambda_i^n}{(1 - \lambda_i)^2} - n \frac{\lambda_i^n}{1 - \lambda_i} = \frac{1}{1 - \lambda_i} \left[ \frac{1}{1 - \lambda_i} - \lambda_i^n \left( \frac{1}{1 - \lambda_i} + n \right) \right]$$

Now, using this $S_i$ in equation A6-12 yields:

$$(A6-13) \qquad T_n = \sum_i \frac{1}{1-\lambda_i} \left[ \frac{1}{1-\lambda_i} - \lambda_i^n \left( \frac{1}{1-\lambda_i} + n \right) \right] x_i e_i$$

Here, simplifying equation A6-13 is possible, but at the expense of another approximation. In the case of an isolated dioecious subpopulation, numerical applications suggested that if $n$ big (i.e. $n>400$ generations) or if the subpopulation is big enough ($N>4$) and $n>10$, then equation A6-13 can be approximated as:

$$(A6-14) \qquad T \approx \sum_i \frac{1}{(1-\lambda_i)^2} x_i e_i$$

For a dioecious population with random mating and even sex-ratio we can write Equation A6-2 (see also A6-3) as: $Q_t = \gamma(A \cdot Q_{t-1} + v)$.

Eigenpairs of matrix $A$ are of the form (see A2-4 or Scripts 1,2 and 4):

$$(A6-15) \qquad \left\{ \begin{array}{l} \lambda_1 = \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} \\ \mathbf{e}_1 = \begin{pmatrix} 1 \\ \lambda_1 \end{pmatrix} \\ \\ \lambda_2 = \frac{1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} \\ \mathbf{e}_2 = \begin{pmatrix} 1 \\ \lambda_2 \end{pmatrix} \end{array} \right.$$

Vector $v$ is composed of a combination of eigenvectors $e_1$ and $e_2$:

<figure><img src="image_108.png" /></figure>

$$v = \left(\frac{0}{2N}\right) = x_1 \left(\frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}\right) + x_2 \left(\frac{1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}\right)$$

<figure><img src="image_109.png" /></figure>

$$\left\{ \begin{array}{l} 0 = x_1 + x_2 \\ \frac{1}{2N} = x_1 \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} + x_2 \frac{1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} \end{array} \right.$$
<figure><img src="image_110.png" /></figure>

$$ \left\{ \begin{aligned} \frac{1}{2N} &= x_1 \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2} - 1 + \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} \\ \frac{x_2}{2N} &= -x_1 \frac{1}{\sqrt{1 + \left(\frac{1}{N}\right)^2}} \end{aligned} \right. $$

<figure><img src="image_111.png" /></figure>

$$ \left\{ \begin{aligned} \frac{x_2}{2N} &= -x_1 \sqrt{1 + \left(\frac{1}{N}\right)^2} = x_1(\lambda_1 - \lambda_2) \\ \frac{1}{2N} &= x_1 \sqrt{1 + \left(\frac{1}{N}\right)^2} = \frac{1}{2N(\lambda_1 - \lambda_2)} \end{aligned} \right. $$

(A6-16)

$$ \left\{ \begin{aligned} x_1 &= \frac{1}{2N\sqrt{1+\left(\frac{1}{N}\right)^2}} = \frac{1}{2N(\lambda_1-\lambda_2)} \\ x_2 &= -\frac{1}{2N\sqrt{1+\left(\frac{1}{N}\right)^2}} = -\frac{1}{2N(\lambda_1-\lambda_2)} \end{aligned} \right. $$

If we combine equations A6-14, and A6-16, we get:

<figure><img src="image_112.png" /></figure>

$$ T = \begin{pmatrix} T_1 \\ T_S \end{pmatrix} \approx \frac{1}{(1-\lambda_1)^2} \frac{1}{2N(\lambda_1-\lambda_2)} \begin{pmatrix} 1 \\ \lambda_1 \end{pmatrix} - \frac{1}{(1-\lambda_2)^2} \frac{1}{2N(\lambda_1-\lambda_2)} \begin{pmatrix} 1 \\ \lambda_2 \end{pmatrix} $$

<figure><img src="image_113.png" /></figure>

$$ \left\{ \begin{aligned} T_1 &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{1}{(1 - \lambda_1)^2} - \frac{1}{(1 - \lambda_2)^2} \right] \\ T_S &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{\lambda_1}{(1 - \lambda_1)^2} - \frac{\lambda_2}{(1 - \lambda_2)^2} \right] \end{aligned} \right. $$

<figure><img src="image_114.png" /></figure>

$$ \left\{ \begin{aligned} T_1 &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{(1 - \lambda_2)^2 - (1 - \lambda_1)^2}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \\ T_S &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{\lambda_1(1 - \lambda_2)^2 - \lambda_2(1 - \lambda_1)^2}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \end{aligned} \right. $$

<figure><img src="image_115.png" /></figure>

$$ \left\{ \begin{aligned} T_1 &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{1 + \lambda_2^2 - 2\lambda_2 - 1 - \lambda_1^2 + 2\lambda_1}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \\ T_S &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{\lambda_1(1 + \lambda_2^2 - 2\lambda_2) - \lambda_2(1 + \lambda_1^2 - 2\lambda_1)}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \end{aligned} \right. $$

<figure><img src="image_116.png" /></figure>

$$ \left\{ \begin{aligned} T_1 &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{\lambda_2^2 - \lambda_1^2 - 2\lambda_2 + 2\lambda_1}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \\ T_S &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{\lambda_1 + \lambda_2^2\lambda_1 - 2\lambda_1\lambda_2 - \lambda_2 - \lambda_1^2\lambda_2 + 2\lambda_1\lambda_2}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \end{aligned} \right. $$

<figure><img src="image_117.png" /></figure>

$$ \left\{ \begin{aligned} T_1 &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{(\lambda_2 - \lambda_1)(\lambda_2 + \lambda_1) + 2(\lambda_1 - \lambda_2)}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \\ T_S &\approx \frac{1}{2N(\lambda_1 - \lambda_2)} \left[ \frac{\lambda_1 - \lambda_2 + \lambda_1\lambda_2(\lambda_2 - \lambda_1)}{(1 - \lambda_1)^2(1 - \lambda_2)^2} \right] \end{aligned} \right. $$
<figure><img src="image_118.png" /></figure>

$$ (A6-17) \qquad \left\{ \begin{array}{l} T_I \approx \frac{(\lambda_1 - \lambda_2)}{2N(\lambda_1 - \lambda_2)} \left[ \frac{2 - \lambda_1 - \lambda_2}{(1 - \lambda_1)^2 (1 - \lambda_2)^2} \right] \\ T_S \approx \frac{\lambda_1 - \lambda_2}{2N(\lambda_1 - \lambda_2)} \left[ \frac{1 - \lambda_1 \lambda_2}{(1 - \lambda_1)^2 (1 - \lambda_2)^2} \right] \end{array} \right. $$ $$ \left\{ \begin{array}{l} T_I \approx \frac{1}{2N} \left[ \frac{2 - \lambda_1 - \lambda_2}{(1 - \lambda_1)^2 (1 - \lambda_2)^2} \right] \\ T_S \approx \frac{1}{2N} \left[ \frac{1 - \lambda_1 \lambda_2}{(1 - \lambda_1)^2 (1 - \lambda_2)^2} \right] \end{array} \right. $$

If we use A6-15 in A6-17, we obtain:

<figure><img src="image_119.png" /></figure>

$$ \left\{ \begin{array}{l} T_I \approx \frac{1}{2N} \left[ \frac{2 - \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} - \frac{1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}}{\left(1 - \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}\right)^2 \left(1 - \frac{1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}\right)^2} \right] \\ T_S \approx \frac{1}{2N} \left[ \frac{1 - \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} \frac{1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}}{\left(1 - \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}\right)^2 \left(1 - \frac{1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2}\right)^2} \right] \end{array} \right. $$

<figure><img src="image_120.png" /></figure>

$$ \left\{ \begin{array}{l} T_I \approx \frac{1}{2N} \left[ \frac{2 - \frac{1 - \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}}{2} + 1 - \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}}{\left(1 + \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}\right)^2 \left(1 + \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}\right)^2} \right] \\ T_S \approx \frac{1}{2N} \left[ \frac{1 - \frac{\left(1 - \frac{1}{N}\right)^2 - 1 - \left(\frac{1}{N}\right)^2}{4}}{\left(1 + \frac{1}{N} - \sqrt{1 + \left(\frac{1}{N}\right)^2}\right)^2 \left(1 + \frac{1}{N} + \sqrt{1 + \left(\frac{1}{N}\right)^2}\right)^2} \right] \end{array} \right. $$
$$ egin{cases} T_I \approx \frac{1}{2N} \frac{2 - \left(1 - \frac{1}{N}\right)}{\left(\frac{\left(1 + \frac{1}{N}\right)^2}{4} - 1 - \left(\frac{1}{N}\right)^2\right)^2} \\ T_S \approx \frac{1}{2N} \left[ \frac{1 - \frac{1 + \left(\frac{1}{N}\right)^2 - 2\frac{1}{N} - 1 - \left(\frac{1}{N}\right)^2}{4}}{\left(\frac{\left(1 + \frac{1}{N}\right)^2}{4} - 1 - \left(\frac{1}{N}\right)^2\right)^2} \right] \end{cases} $$

$$ egin{cases} T_I \approx \frac{1}{2N} \frac{1 + \frac{1}{N}}{\left(\frac{1 + \left(\frac{1}{N}\right)^2 + 2\frac{1}{N} - 1 - \left(\frac{1}{N}\right)^2}{4}\right)^2} \\ T_S \approx \frac{1}{2N} \left[ \frac{1 + \frac{1}{2N}}{\left(\frac{1 + \left(\frac{1}{N}\right)^2 + 2\frac{1}{N} - 1 - \left(\frac{1}{N}\right)^2}{4}\right)^2} \right] \end{cases} $$

$$ egin{cases} T_I \approx \frac{1}{2N} \frac{1 + \frac{1}{N}}{\left(\frac{1}{2N}\right)^2} \\ T_S \approx \frac{1}{2N} \frac{1 + \frac{1}{2N}}{\left(\frac{1}{2N}\right)^2} \end{cases} $$

$$ egin{cases} T_I \approx 2N \left(1 + \frac{1}{N}\right) \\ T_S \approx 2N \left(1 + \frac{1}{2N}\right) \end{cases} $$

$$ (A6-18) \qquad \begin{cases} T_I \approx 2(N + 1) \\ T_S \approx 2N + 1 \end{cases} $$

This result is the same as in Balloux's paper (Balloux, 2004) (equation 15) with an even sex ratio.
For a panmictic population of size $N_e$:
<figure><img src="image_121.png" /></figure>

$$ egin{cases} Q_{1t} = \frac{1}{N_e} \left( \frac{1}{2} + \frac{1}{2} Q_{1t-1} \right) + \left( 1 - \frac{1}{N_e} \right) Q_{S_{t-1}} \\ Q_{S_t} = \frac{1}{N_e} \left( \frac{1}{2} + \frac{1}{2} Q_{1t-1} \right) + \left( 1 - \frac{1}{N_e} \right) Q_{S_{t-1}} \end{cases} $$

<figure><img src="image_122.png" /></figure>

$$ egin{cases} Q_{1t} = \frac{1}{2N_e} Q_{S_{t-1}} + \left(1 - \frac{1}{N_e}\right) Q_{S_{t-1}} + \frac{1}{2N_e} \\ Q_{S_t} = \frac{1}{2N_e} Q_{S_{t-1}} + \left(1 - \frac{1}{N_e}\right) Q_{S_{t-1}} + \frac{1}{2N_e} \end{cases} $$

$$ egin{cases} Q_{1t} = \left(1 - \frac{1}{2N_e}\right) Q_{S_{t-1}} + \frac{1}{2N_e} \\ Q_{S_t} = \left(1 - \frac{1}{2N_e}\right) Q_{S_{t-1}} + \frac{1}{2N_e} \end{cases} $$

The corresponding transition matrix has the following eigenpair:

(A6-19)

$$ egin{cases} \lambda_{e1} = 1 - \frac{1}{2N_e} \\ \lambda_{e2} = 0 \\ \mathbf{e}_{e1} = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \\ \mathbf{e}_{e2} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \end{cases} $$

Hence:

$$ \mathbf{v}_e = \begin{pmatrix} \frac{1}{2N_e} \\ \frac{1}{2N_e} \end{pmatrix} = x_{e1} \begin{pmatrix} 1 \\ 1 \end{pmatrix} + x_{e2} \begin{pmatrix} 1 \\ 0 \end{pmatrix} $$

<figure><img src="image_123.png" /></figure>

$$ egin{cases} \frac{1}{2N_e} = x_{e1} + x_{e2} \\ \frac{1}{2N_e} = x_{e1} \end{cases} $$

<figure><img src="image_124.png" /></figure>

(A6-20)

$$ egin{cases} x_{e1} = \frac{1}{2N_e} \\ x_{e2} = 0 \end{cases} $$

If we apply equation A6-14 with the values of equations A6-19 and A6-20), we obtain:

$$ egin{cases} T_{1e} \approx \frac{1}{\left(1 - 1 - \frac{1}{2N_e}\right)^2} \frac{1}{2N_e} \\ T_{Se} \approx \frac{1}{\left(1 - 1 - \frac{1}{2N_e}\right)^2} \frac{1}{2N_e} \end{cases} $$

<figure><img src="image_125.png" /></figure>

(A6-21)

$$ egin{cases} T_{1e} \approx 2N_e \\ T_{Se} \approx 2N_e \end{cases} $$
From there, it is easy to understand that the coalescent effective population size can then be defined as in equation 17 of Balloux et al. (2003), i.e.:

$$ (A6-22) \qquad N_e \approx \frac{1}{2} \bar{T} $$

where $\bar{T}$ is the weighted average of the different $T_i$'s, here:

$$ (A6-23) \qquad \bar{T} = \frac{1}{N} T_1 + \left(1 - \frac{1}{N}\right) T_S $$

The weights in fact correspond to the probabilities to sample two genes from the considered hierarchy: within one individual, from two different individuals within the same sub-population, from two different sub-populations from the same archipelago, etc...

In our context, for a dioecious and isolated population of size $N$ with an even sex-ratio, combining equations A6-18, A6-22 and A6-23 leads to:

<figure><img src="image_126.png" /></figure>

$$ N_e \approx \frac{1}{2} \left[ \frac{1}{N} 2(N + 1) + \left(1 - \frac{1}{N}\right) (2N + 1) \right] $$

<figure><img src="image_127.png" /></figure>

$$ N_e \approx \frac{1}{2} \left[ 2 \left( 1 + \frac{1}{N} \right) + 2N + 1 - 2 - \frac{1}{N} \right] $$

<figure><img src="image_128.png" /></figure>

$$ N_e \approx \frac{1}{2} \left[ 2 + \frac{2}{N} + 2N - 1 - \frac{1}{N} \right] $$

<figure><img src="image_129.png" /></figure>

$$ N_e \approx \frac{1}{2} \left[ 1 + \frac{1}{N} + 2N \right] $$

$$ (A6-24) \qquad N_e \approx N + \frac{1}{2} + \frac{1}{2N} $$

Equation A6-24 is exactly the same as equation 10 in Balloux (2004). To give a biological meaning to this result, it corresponds to the census size (or number of breeders) plus half an individual that would have been coalescent through random selfing in a WF population, plus one coalescent individual that would occur in a WF population (which may sound redundant with the second).

## Appendix 7: Matrix method to compute eigenvalue effective population size in a dioecious population

Let $Q_{I(t)}$ and $Q_{S(t)}$ be the probabilities of identity between two alleles at time $t$ within individuals and between individuals in a dioecious random mating population. We can then use Equations 7 and 8 in the main text:

<figure><img src="image_130.png" /></figure>

$$ (A7-1) \qquad \left\{ \begin{array}{l} Q_{I(t)} = Q_{S(t-1)} \\ Q_{S(t)} = Q_{I(t-1)} \frac{1}{8} \left( \frac{1}{N_f} + \frac{1}{N_m} \right) + Q_{S(t-1)} \left( 1 - \frac{1}{4N_f} - \frac{1}{4N_m} \right) + \frac{1}{8} \left( \frac{1}{N_f} + \frac{1}{N_m} \right) \end{array} \right. $$

Equation A9-1 has transition matrix (see appendices 1-4):
$$ (A7-2) \qquad \mathbf{A} = \left( \frac{1}{2} \left( \frac{1}{4N_f} + \frac{1}{4N_m} \right)^0 + \left( 1 - \frac{1}{4N_f} - \frac{1}{4N_m} \right)^1 \right) $$

To save time we used wxMaxima to find the leading eigenvalue of A (see Script 1):

$$ \lambda_1 = \frac{\sqrt{(16N_f^2 + 1)N_m^2 + 2N_f N_m + N_f^2 + (4N_f - 1)N_m - N_f}}{8N_f N_m} $$

<figure><img src="image_131.png" /></figure>

$$ \lambda_1 = \frac{1}{2} - \frac{N}{8N_f N_m} + \frac{\sqrt{(4N_f N_m)^2 + N_f^2 + N_m^2 + 2N_f N_m}}{8N_f N_m} $$

<figure><img src="image_132.png" /></figure>

$$ \lambda_1 = \frac{1}{2} - \frac{N}{8N_f N_m} + \frac{1}{2} \sqrt{1 + \left( \frac{N}{4N_f N_m} \right)^2} $$

which is the same as equation 11 in the main text (QED).

## Appendix 8: Derivatives and Taylor-MacLaurin's expansion series

### Basic notions about derivative functions

Readers acquainted with derivatives can skip this first section.

The derivative of a function $f(x)$ describes the orientation and speed of variation of this function, measured between two points separated by a distance $\Delta x$ that tends to 0:

$$ f'(x) = \lim_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x)}{\Delta x} $$

For the present paper, we will need to compute the derivative of several functions. For instance for the function $f(x)=x^n$, then:

$$ f'(y) = \lim_{\Delta x \to 0} \frac{(x + \Delta x)^n - x^n}{\Delta x} $$

For any $n$, beginning with $n=2$ or 3, it is easy to show that:

$$ f'(x) = \lim_{\Delta x \to 0} nx^{(n-1)} + \Delta x \times g(x) $$

where $g(x)$ is a function of $x$ with one term in $x^{m \cdot n}$, one term in $\Delta x^{n-2}$ and other terms in $x\Delta x$, so that the limit when $\Delta x \to 0$ necessary is $nx^{n-1}$. Then $f'(x)=nx^{n-1}$. It is easy to see that the derivative of a sum of functions is simply the sum of derivatives of the different functions of this sum.

Next, we need to compute the derivative of $f(g(x))$ or more correctly $(f \circ g)(x)$. For this, it will be easier to change of notation:

$$ (f \circ g)'(x) = \lim_{\Delta x \to 0} \frac{f(g(x) + \Delta x) - f(g(x))}{\Delta x} $$

Let $f(g(x) + \Delta x) - f(g(x)) = \Delta u$, and $\Delta v = g(x + \Delta x) - g(x)$, then:

$$ (f \circ g)'(x) = \lim_{\Delta x \to 0} \frac{\Delta u}{\Delta x} $$

<figure><img src="image_133.png" /></figure>

$$ (f \circ g)'(x) = \lim_{\Delta x \to 0} \frac{\Delta u \Delta v}{\Delta v \Delta x} $$

<figure><img src="image_134.png" /></figure>

$$ (f \circ g)'(x) = g'(x) \lim_{\Delta x \to 0} \frac{\Delta u}{\Delta v} $$

Since it is easy to see that when $\Delta x \to 0$, then $\Delta v \to 0$, we can rewrite:
$$ (f \circ g)'(x) = f'(u)g'(x) $$

or

$$ (f \circ g)'(x) = f'(g(x))g'(x) $$

We then have the necessary tools for Taylor-MacLaurin's expansion series.

### Taylor-MacLaurin's expansion series

In the neighborhood of a, any infinitely differentiable function f(x), writes:

$$ f(x) = f(a) + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \cdots + \frac{f^{(n)}(a)}{n!}(x-a)^n + \varepsilon $$

Indeed, let f be a derivable function of variable x so that:

$$ f(x) = a_0 + a_1(x-a) + a_2(x-a)^2 + a_3(x-a)^3 + a_4(x-a)^4 + \cdots $$

If we derivate f, we get:

$$ f'(x) = a_1 + 2a_2(x-a) + 3a_3(x-a)^2 + 4a_4(x-a)^3 + \cdots $$

$$ f''(x) = 2a_2 + 3 \times 2 \times a_3(x-a) + 4 \times 3 \times a_4(x-a)^2 + \cdots $$

$$ f'''(x) = 3 \times 2 \times a_3 + 4 \times 3 \times 2 \times a_4(x-a) + \cdots $$

If $x \to a$, then $(x-a) \to 0$ and:

$$ f'(x) = a_1 $$

$$ f''(x) = 2a_2 $$

$$ f'''(x) = 3 \times 2 \times a_3 $$

$$ f''''(x) = 4 \times 3 \times 2 \times a_3 $$

$$ f^{(n)}(x) = n! \times a_n $$

where $f^{(n)}$ is the nth derivative of f.

We can thus set that, in the neighborhood of a:

$$ a_n = \frac{f^{(n)}(x)}{n!} $$

From there we can rewrite $f(x)$ in the neighborhood of a:

$$ f(x) = \frac{f(a)}{0!} + \frac{f'(a)}{1!}(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \frac{f'''(a)}{3!}(x-a)^3 + \cdots $$

If $a \to 0$, $f(x)$ writes (QED):

$$ f(x) = f(0) + f'(0)(x) + \frac{f''(0)}{2!}(x)^2 + \frac{f'''(0)}{3!}(x)^3 + \cdots $$

This method will offer very good approximations of cumbersome functions of small variables (e.g. $1/(2N)$ or $1/N^2$).

### Examples of Taylor-MacLaurin's expansion series

We need to find a proxy for $\sqrt{1 + \bar{x}}$ and for $1/(1-\bar{x})$, when X is small.

For $f(x) = \sqrt{x}$:

$$ f'(x) = \lim_{\Delta x \to 0} \frac{\sqrt{x + \Delta x} - \sqrt{x}}{\Delta x} $$

<figure><img src="image_135.png" /></figure>

$$ f'(x) = \lim_{\Delta x \to 0} \frac{[\sqrt{x + \Delta x} - \sqrt{x}][\sqrt{x + \Delta x} + \sqrt{x}]}{\Delta x[\sqrt{x + \Delta x} + \sqrt{x}]} $$

<figure><img src="image_136.png" /></figure>

$$ f'(x) = \lim_{\Delta x \to 0} \frac{x + \Delta x + \sqrt{x + \Delta x}\sqrt{x} - \sqrt{x}\sqrt{x + \Delta x} - x}{\Delta x[\sqrt{x + \Delta x} + \sqrt{x}]} $$
<figure><img src="image_137.png" /></figure>

$$f'(x) = \lim_{\Delta x \to 0} \frac{\Delta x}{\Delta x[\sqrt{x + \Delta x} + \sqrt{x}]}$$

<figure><img src="image_138.png" /></figure>

$$f'(x) = \lim_{\Delta x \to 0} \frac{1}{\sqrt{x + \Delta x} + \sqrt{x}}$$

<figure><img src="image_139.png" /></figure>

$$f'(x) = \frac{1}{2\sqrt{x}}$$

Let $g(X)$ be:

$$g(X) = \sqrt{1 + X}$$

Then:

$$g'(X) = (\sqrt{1 + X})' \times 1$$

<figure><img src="image_140.png" /></figure>

$$g'(X) = \frac{1}{2\sqrt{1 + X}}$$

For the function $1/x$:

$$(\frac{1}{x})' = \lim_{\Delta x \to 0} \frac{(\frac{1}{x + \Delta x}) - \frac{1}{x}}{\Delta x}$$

<figure><img src="image_141.png" /></figure>

$$(\frac{1}{x})' = \lim_{\Delta x \to 0} \frac{\frac{x - x - \Delta x}{x(x + \Delta x)}}{\Delta x}$$

<figure><img src="image_142.png" /></figure>

$$(\frac{1}{x})' = \lim_{\Delta x \to 0} \frac{-\Delta x}{x(x + \Delta x)\Delta x}$$

<figure><img src="image_143.png" /></figure>

$$(\frac{1}{x})' = \lim_{\Delta x \to 0} \frac{-1}{x(x + \Delta x)}$$

<figure><img src="image_144.png" /></figure>

$$(\frac{1}{x})' = \frac{-1}{x^2}$$

Using the same approach as for $g'(X)$:

$$g''(X) = - \frac{1}{2^2 (\sqrt{1 + X})^3}$$

If we now use Taylor for $g(X) = \sqrt{1 + X}$ in the neighborhood of $a$:

$$g(X) \approx \frac{g(a)}{0!} (X - a)^0 + \frac{g'(a)}{1!} (X - a)^1 + \frac{g''(a)}{2!} (X - a)^2 + \dots$$

<figure><img src="image_145.png" /></figure>

$$g(X) \approx \frac{\sqrt{1 + a}}{0!} (X - a)^0 + \frac{1}{2\sqrt{1 + a}} \frac{1}{1!} (X - a)^1 + \frac{1}{2^2(\sqrt{1 + a})^3} \frac{1}{2!} (X - a)^2 + \dots$$

When $a \to 0$, we get:

$$g(X) = \sqrt{1 + X} = 1 + \frac{1}{2}X - \frac{1}{8}X^2 + \dots$$

The same result can be obtained with the command "taylor(sqrt(1+X),X,0,3)" in wxMaxima (Vodopivec, 2017).
$$g(X) = \sqrt{1 + X} = 1 + \frac{1}{2}X - \frac{1}{8}X^2 + \frac{1}{16}X^3 + \dots$$

Now, if $X \ll 1$, we can approximate this expression as:

$$\sqrt{1 + X} \approx 1 + \frac{1}{2}X - \frac{1}{8}X^2$$

The first and second derivatives of $1/(1-X)$ are $1/(1-X)^2$ and $2(1-X)/(1-X)^4$. We can use Taylor-MacLaurin again and write: $1/(1-X)=1+X+X^2+X^3+\dots$ (note that the same result would have been obtained with Maxima (Vodopivec, 2017) typing "taylor(1/(1-X),X,0,3)").

## Appendix 9: Finding the root of $N_{e\_Eq3}-N_{e\_Eq13}$

Since Equation 15 gave an almost perfect estimate of equation 13, we studied the sign of $\Delta N_e=N_{e\_Eq3}-N_{e\_Eq15}$ instead of $N_{e\_Eq3}-N_{e\_Eq13}$ for the sake of simplicity.

$$\Delta N_e = \frac{4N_f N_m}{N} + \frac{1}{2} + \frac{1}{2N} - \frac{4N_f N_m}{N} - \frac{1}{2} - \frac{1}{4} \frac{N}{4N_f N_m}$$

⇒

$$\Delta N_e = \frac{1}{2N} - \frac{1}{4} \frac{N}{4N_f N_m}$$

->

$$\Delta N_e \propto \frac{1}{N} - \frac{N}{8N_f N_m}$$

->

$$\Delta N_e \propto \frac{8N_f N_m - N^2}{8N_f N_m N}$$

->

$$\Delta N_e \propto 8N_f N_m - (N_f + N_m)^2$$

->

$$\Delta N_e \propto 8N_f N_m - N_f^2 - N_m^2 - 2N_f N_m$$

->

$$\Delta N_e \propto 6N_f N_m - N_f^2 - N_m^2$$

We can divide the right term by $N_f^2$, then, noting $SR=N_m/N_f$:

$$\Delta N_e \propto 6SR - 1 - SR^2$$

->

$$\Delta N_e \propto SR^2 - 6SR + 1$$

->

$$\Delta N_e \propto SR^2 - 2 \times 3SR + 9 - 8$$

->

$$\Delta N_e \propto (SR - 3)^2 - 8$$

We need to find the two roots of the right term, which must satisfy:

$$(SR - 3)^2 = 8$$

⇒

$$SR - 3 = \pm \sqrt{8}$$

⇒

$$SR = 3 \pm \sqrt{8}$$

⇒

$$\begin{cases} SR_1 = 3 + \sqrt{8} \approx 5.8284 \\ SR_2 = 3 - \sqrt{8} \approx 0.1716 \end{cases}$$

A sex-ratio above 1 is not relevant here, since an excess of females would lead to the same result as in populations with an even sex-ratio. Consequently, $SR_2$ is the only relevant root. From there, it can be seen
that Balloux's equation will provide an over-estimate when $SR>SR_2$, an under-estimate when $SR<SR_2$ and will be exact when $SR=SR_2=3-2\sqrt{2}$.

## Appendix 10: Balloux's like method to compute $F_{IS}$ based $N_e$

Let $Q_t$ and $Q_{t-1}$ be the probabilities to sample twice the same allele in one individual and between individuals from the same population, respectively, $u$ the mutation rate, then, for a dioecious population with an even sex-ratio and random mating, we can set the following recurrences between generation $t$ and $t-1$ (see equations 7 and 8 with $N_f=N_m$):

$$ (A10-1) \qquad \left\{ \begin{array}{l} Q_{I(t)} = (1-u)^2 Q_{S(t-1)} \\ Q_{S(t)} = (1-u)^2 \left[ \frac{1}{N} \left( \frac{1}{2} + \frac{1}{2} Q_{I(t-1)} \right) + \left( 1 - \frac{1}{N} \right) Q_{S(t-1)} \right] \\ \\ Q_{I(t)} = (1-u)^2 Q_{S(t-1)} \\ Q_{S(t)} = (1-u)^2 \frac{1}{2N} Q_{I(t-1)} + (1-u)^2 \left( 1 - \frac{1}{N} \right) Q_{S(t-1)} + (1-u)^2 \frac{1}{2N} \end{array} \right. $$

Let $Q_t$ and $Q_{t-1}$ be the vectors defining genetic identities at generations $t$ and $t-1$, $A$ the transition matrix and $v$ the vector of residuals, then:

$$ (A10-2) \qquad \left\{ \begin{array}{l} Q_t = \begin{pmatrix} Q_{I(t)} \\ Q_{S(t)} \end{pmatrix} \\ A = \begin{pmatrix} 0 & (1-u)^2 \\ (1-u)^2 \frac{1}{2N} & (1-u)^2 \left(1-\frac{1}{N}\right) \end{pmatrix} \\ v = \begin{pmatrix} 0 \\ (1-u)^2 \frac{1}{2N} \end{pmatrix} \end{array} \right. $$

and $Q_t=A.Q_{t-1}+v$.

At equilibrium, we can write that the vector of genetic identities $Q$ writes $Q=(I-A)^{-1}.V$, where $I=\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$ is the identity matrix (see appendix 5).

To solve this equation, and get $Q_t$ and $Q_{t-1}$ at equilibrium, we used wxMaxima 17.10.1 (Vodopivec, 2017) as detailed in the section wxMaxima scripts, Script 2. Taking into account that $u<<1$, we obtained:

$$ (A10-3) \qquad F_{IS} = \frac{Q_I - Q_S}{1 - Q_S} \approx -\frac{1}{2N+1} $$

The same results can be obtained with classic algebra, without the use of matrix computations, but it is much faster this way. This is also the same results as equation 8 in Balloux's paper. It is worth mentioning here that equation A10-3 can also theoretically give access to the census size of individuals in the population ($N$) or, more precisely, to the exact number of adult parents of the individuals in the population, that some may call the effective number of breeders:

$$ (A10-4) \qquad N = -\frac{1+F_{IS}}{2F_{IS}} $$

If we go back to equation 16 of the main manuscript, we can compute the eigenvalue effective population size as:

$$ (A10-5) \qquad N_e \approx N + \frac{1}{2} + \frac{1}{4N} $$

If we combine equations A9-5 with A9-4, we obtain:
$$N_e \approx -\frac{1 + F_{IS}}{2F_{IS}} + \frac{1}{2} - \frac{1}{4\frac{1 + F_{IS}}{2F_{IS}}}$$

⇔

$$(A10-6) \qquad N_e \approx -\frac{1}{2F_{IS}} - \frac{F_{IS}}{2(1+F_{IS})}$$

Now, with a stronger approximation, $N_e \approx N+1/2$, which, combined with equation A9-4 yields Pudovkin et al.'s equation 3 (see equation 4 of the present manuscript). We can also notice that A9-6 is the average of equations 4 (Pudovkin et al. second equation) and 6 (Balloux).

## Appendix 11: Equilibrium value for $F_{IS}$ in a dioecious population (general case)

For this we need to use equation A6-1 and add a mutation rate $u$ so that equation A6-1 becomes:

$$(A10-7) \qquad \left\{ \begin{aligned} Q_{I(t)} &= Q_{S(t-1)}(1-u)^2 \\ Q_{S(t)} &= Q_{I(t-1)} \frac{1}{8} \left( \frac{1}{N_f} + \frac{1}{N_m} \right) (1-u)^2 + Q_{S(t-1)} \left( 1 - \frac{1}{4N_f} - \frac{1}{4N_m} \right) (1-u)^2 \\ &\quad + \frac{1}{8} \left( \frac{1}{N_f} + \frac{1}{N_m} \right) (1-u)^2 \end{aligned} \right.$$

The transition matrix and the associated vectors of this equation are:

$$\left\{ \begin{aligned} \mathbf{Q}_t &= \begin{pmatrix} Q_{I(t)} \\ Q_{S(t)} \end{pmatrix} \\ \mathbf{A} &= \begin{pmatrix} 0 & (1-u)^2 \\ \frac{1}{2}\left(\frac{1}{4N_f} + \frac{1}{4N_m}\right)(1-u)^2 & \left(1 - \frac{1}{4N_f} - \frac{1}{4N_m}\right)(1-u)^2 \end{pmatrix} \\ \mathbf{V} &= \begin{pmatrix} 0 \\ \frac{1}{8}\left(\frac{1}{N_f} + \frac{1}{N_m}\right)(1-u)^2 \end{pmatrix} \end{aligned} \right.$$

and equation A7-1 can be rewritten as $\mathbf{Q}_t = \mathbf{A}\mathbf{Q}_{t-1} + \mathbf{V}$.

At equilibrium, the vector of genetic identities satisfies the equation $\mathbf{Q}=(\mathbf{I}-\mathbf{A})^{-1}\mathbf{V}$, where $\mathbf{I}=\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$ is the identity matrix (see appendix 5).

To solve this equation, and get $Q_I$ and $Q_S$ at equilibrium, we used wxMaxima 17.10.1 (Vodopivec, 2017), as detailed in the section wxMaxima scripts (Script 3) and obtained:

$$F_{IS} - \frac{(N_m + N_f)u^2 - 2(N_m + N_f)u + N_m + N_f}{(N_m + N_f)u^2 - 2(N_m + N_f)u + N_m + N_f + 8N_m N_f}$$

⇔

$$(A10-8) \qquad F_{IS} = - \frac{(N_m + N_f)(1-u)^2}{(N_m + N_f)(1-u)^2 + 8N_m N_f}$$

Terms in $u$ are small in front of 1 so that equation A7-2 can be simplified as:

$$(A10-9) \qquad F_{IS} \approx - \frac{N_m + N_f}{N_m + N_f + 8N_m N_f}$$

## Appendix 12: Effective population size of an isolated monogamous population

We will use the same notations as in other sections. Monogamy implies an even sex ratio in the pool of adults that are involved in a mating. For the identity within individuals, the recurrence stays the same as in polygamous populations. The recursion for the identity between individuals can be determined by conditioning on the ancestry of the sampled pair in the previous generation. One possibility is that the two
sampled individuals are sibs, i.e., they share the same parents, which is true with probability 1/(N/2). In this case, with probability 1/2, the two alleles will have come from the same parent, in which case they are equally likely to be derived from a single parental allele or from both parental alleles. In the former case, the sampled alleles are necessarily IBD, whereas in the latter case, the probability that they are IBD is $Q_{I(t-1)}$. Alternatively, with probability 1/2, each sampled allele may have come from a different parent, in which case the probability that they are IBD is $Q_{S(t-1)}$. The second possibility, which has probability 1 - 1/(N/2), is that the two sampled individuals are not sibs, in which case the probability that the sampled alleles are IBD is $Q_{S(t-1)}$. We can thus set the following recurrence:

$$ \left\{ \begin{aligned} Q_{S(t)} &= \frac{1}{\frac{1}{2}N} \left[ \frac{1}{2} \left( \frac{1}{2} + \frac{1}{2} Q_{I(t-1)} \right) + \frac{1}{2} Q_{S(t-1)} \right] + \left( 1 - \frac{1}{\frac{1}{2}N} \right) Q_{S(t-1)} \\ Q_{I(t)} &= Q_{S(t-1)} \end{aligned} \right. $$

<figure><img src="image_146.png" /></figure>

$$ \left\{ \begin{aligned} Q_{S(t)} &= Q_{I(t-1)} \frac{1}{2N} + Q_{S(t-1)} \left( \frac{1}{N} + 1 - \frac{1}{\frac{1}{2}N} \right) + \frac{1}{2N} \\ Q_{I(t)} &= Q_{S(t-1)} \end{aligned} \right. $$

<figure><img src="image_147.png" /></figure>

$$ \left\{ \begin{aligned} Q_{S(t)} &= Q_{I(t-1)} \frac{1}{2N} + Q_{S(t-1)} \left( 1 - \frac{1}{N} \right) + \frac{1}{2N} \\ Q_{I(t)} &= Q_{S(t-1)} \end{aligned} \right. $$

Using Maxima, it is easy to compute the leading eigenvalue of the corresponding transition matrix as:

$$ \lambda_1 = \frac{1}{2} - \frac{1}{2N} + \frac{1}{2} \sqrt{1 + \left(\frac{1}{N}\right)^2} $$

For X small, Taylor-MacLaurin of $\sqrt{1 + X} \approx 1 + \frac{1}{2}X$, hence:

$$ \lambda_1 \approx \frac{1}{2} - \frac{1}{2N} + \frac{1}{2} \left[ 1 + \frac{1}{2} \left( \frac{1}{N} \right)^2 \right] $$

<figure><img src="image_148.png" /></figure>

$$ \lambda_1 \approx 1 - \frac{1}{2N} \left( 1 - \frac{1}{2N} \right) $$

The eigenvalue effective population size is:

$$ N_e \approx \frac{1}{2(1 - \lambda_1)} $$

<figure><img src="image_149.png" /></figure>

$$ N_e \approx \frac{1}{2 \left( 1 - 1 + \frac{1}{2N} \left( 1 - \frac{1}{2N} \right) \right)} $$

<figure><img src="image_150.png" /></figure>

$$ N_e \approx \frac{1}{\frac{1}{N} \left( 1 - \frac{1}{2N} \right)} $$

<figure><img src="image_151.png" /></figure>

$$ N_e \approx \frac{N}{1 - \frac{1}{2N}} $$

<figure><img src="image_152.png" /></figure>

Using Taylor-MacLaurin again leads to:
$$ \frac{1}{1 - \frac{1}{2N}} \approx 1 + \frac{1}{2N} + \left(\frac{1}{2N}\right)^2 $$

Then:

$$ N_e \approx N \left[ 1 + \frac{1}{2N} + \left( \frac{1}{2N} \right)^2 \right] $$

<figure><img src="image_153.png" /></figure>

$$ N_e \approx N + \frac{1}{2} + \frac{1}{2N} $$

This result is exactly the same as for a dioecious pangamic population with even sex ratio.

## wxMaxima scripts

### Script 1: Computing the eigenvalues of Matrix A (equation A6-2)

(%i1) A: matrix([0,1], [(1/N_f+1/N_m)/8,1-1/(4·N_f)-1/(4·N_m)]);

$$ (A) \begin{pmatrix} 0 & 1 \\ \frac{1}{N_m} + \frac{1}{N_f} & -\frac{1}{4N_m} - \frac{1}{4N_f} + 1 \end{pmatrix} $$

(%i2) eigenvalues(A);

$$ (\%02) \quad [[-\frac{\sqrt{(16N_f^2+1)N_m^2+2N_fN_m+N_f^2+(1-4N_f)N_m+N_f}}{8N_fN_m}, \frac{\sqrt{(16N_f^2+1)N_m^2+2N_fN_m+N_f^2+(4N_f-1)N_m-N_f}}{8N_fN_m}], [1,1]] $$

### Script 2: recomputing N_e and FIS in dioecious populations with an even sex ratio

(%i1) A: matrix([0,(1-u)^2], [(1-u)^2/(2·N),(1-u)^2·(1-1/N)]);

$$ (A) \begin{pmatrix} 0 & (1-u)^2 \\ \frac{(1-u)^2}{2N} & (1 - \frac{1}{N})(1-u)^2 \end{pmatrix} $$

(%i2) I: matrix([1,0], [0,1]);

$$ (I) \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} $$

(%i3) V: matrix([0], [(1-u)^2/(2·N)]);

$$ (V) \begin{pmatrix} 0 \\ \frac{(1-u)^2}{2N} \end{pmatrix} $$

(%i4) Q:invert(I-A).V;

$$ (Q) \begin{pmatrix} \frac{(1-u)^4}{2N} & \frac{-(1-u)^4}{2N} - \frac{1}{N}(1-u)^2 + 1 \\ \frac{1-u^2}{2N} & \frac{-(1-u)^4}{2N} - \frac{1}{N}(1-u)^2 + 1 \end{pmatrix} $$

(%i5) QI:(1-u)^4/(2·N·(-(1-u)^4/(2·N)-(1-1/N)·(1-u)^2+1));

$$ (QI) \frac{(1-u)^4}{2N} \frac{-(1-u)^4}{2N} - \frac{1}{N}(1-u)^2 + 1 $$

(%i6) QS:(1-u)^2/(2·N·(-(1-u)^4/(2·N)-(1-1/N)·(1-u)^2+1));
$$ (\text{QS}) \frac{(1-u)^2}{2N\left(-\frac{(1-u)^4}{2N}-\left(1-\frac{1}{N}\right)(1-u)^2+1\right)} $$

$$ (\%i7) \text{ FIS:} \frac{\text{QI-SQ}}{(1-\text{QS})}; $$

$$ (\text{FIS}) \frac{\frac{(1-u)^4}{2N\left(-\frac{(1-u)^4}{2N}-\left(1-\frac{1}{N}\right)(1-u)^2+1\right)} - \frac{(1-u)^2}{2N\left(-\frac{(1-u)^4}{2N}-\left(1-\frac{1}{N}\right)(1-u)^2+1\right)}{1-\frac{(1-u)^2}{2N\left(-\frac{(1-u)^4}{2N}-\left(1-\frac{1}{N}\right)(1-u)^2+1\right)}} $$

$$ (\%i9) \text{ FIS2:} \text{ratsimp(FIS)}; $$

$$ (\text{FIS2}) - \frac{u^2-2u+1}{u^2-2u+2N+1} $$

$$ (\%i10) \text{ eigenvalues(A)}; $$

$$ (\%i10) \left[ -\frac{\sqrt{N^2+1}(u^2-2u+1)+(1-N)u^2+(2N-2)u-N+1}{2N}, \frac{\sqrt{N^2+1}(u^2-2u+1)+(N-1)u^2+(2-2N)u+N-1}{2N} \right], [1,1] $$

$$ (\%i11) \lambda_1: \text{ratsubst}(0, u, (\sqrt{N^2+1} \cdot (u^2-2u+1) + (N-1) \cdot u^2 + (2-2N) \cdot u + N-1) / (2 \cdot N)); $$

$$ (\lambda1) \frac{\sqrt{N^2+1}+N-1}{2N} $$

## Script 3: Computing FIS in a dioecious population (general case)

$$ (\%i2) \text{ A: matrix}([0, (1-u)^2], [(1-u)^2 \cdot (1/(4 \cdot N_f) + 1/(4 \cdot N_m)) / 2, (1-1/(4 \cdot N_f) - 1/(4 \cdot N_m)) \cdot (1-u)^2]); $$

$$ (A) \begin{pmatrix} 0 & (1-u)^2 \\ \frac{\left(\frac{1}{4N_m} + \frac{1}{4N_f}\right)(1-u)^2}{2} & \left(-\frac{1}{4N_m} - \frac{1}{4N_f} + 1\right)(1-u)^2 \end{pmatrix} $$

$$ (\%i3) \text{ V: matrix}([0], [(1-u)^2 \cdot (1/N_f + 1/N_m) / 8]); $$

$$ (V) \left( \frac{0}{\frac{1}{N_m} + \frac{1}{N_f}} (1-u)^2 \right) $$

$$ (\%i4) \text{ I: matrix}([1,0], [0,1]); $$

$$ (I) \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} $$

$$ (\%i5) \text{ Q: invert}(I-A).V; $$

$$ (Q) \begin{pmatrix} \frac{\left(\frac{1}{N_m} + \frac{1}{N_f}\right)(1-u)^4}{8\left(-\frac{\left(\frac{1}{4N_m} + \frac{1}{4N_f}\right)(1-u)^4}{2} - \left(-\frac{1}{4N_m} - \frac{1}{4N_f} + 1\right)(1-u)^2 + 1\right)} & \frac{\left(\frac{1}{N_m} + \frac{1}{N_f}\right)(1-u)^4}{8\left(-\frac{\left(\frac{1}{4N_m} + \frac{1}{4N_f}\right)(1-u)^4}{2} - \left(-\frac{1}{4N_m} - \frac{1}{4N_f} + 1\right)(1-u)^2 + 1\right)} \\ \frac{\left(\frac{1}{N_m} + \frac{1}{N_f}\right)(1-u)^2}{8\left(-\frac{\left(\frac{1}{4N_m} + \frac{1}{4N_f}\right)(1-u)^4}{2} - \left(-\frac{1}{4N_m} - \frac{1}{4N_f} + 1\right)(1-u)^2 + 1\right)} & \frac{\left(\frac{1}{N_m} + \frac{1}{N_f}\right)(1-u)^4}{8\left(-\frac{\left(\frac{1}{4N_m} + \frac{1}{4N_f}\right)(1-u)^4}{2} - \left(-\frac{1}{4N_m} - \frac{1}{4N_f} + 1\right)(1-u)^2 + 1\right)} \end{pmatrix} $$

$$ (\%i6) \text{ ratsimp}(\%); $$

$$ (\%o6) \frac{\left( - \frac{(N_m + N_f)u^4 + (-4N_m - 4N_f)u^3 + (6N_m + 6N_f)u^2 + (-4N_m - 4N_f)u + N_m + N_f}{(N_m + N_f)u^4 + (-4N_m - 4N_f)u^3 + ((8N_f + 4)N_m + 4N_f)u^2 - 16N_fN_mu - N_m - N_f} \right)}{\left( N_m + N_f \right)u^4 + (-4N_m - 4N_f)u^3 + ((8N_f + 4)N_m + 4N_f)u^2 - 16N_fN_mu - N_m - N_f} $$

$$ (\%i7) \text{ QI:} - \frac{((N_m + N_f)u^4 + (-4 \cdot N_m - 4 \cdot N_f)u^3 + (6 \cdot N_m + 6 \cdot N_f)u^2 + (-4 \cdot N_m - 4 \cdot N_f)u + N_m + N_f)}{((N_m + N_f)u^4 + (-4 \cdot N_m - 4 \cdot N_f)u^3 + ((8 \cdot N_f + 4) \cdot N_m + 4 \cdot N_f)u^2 - 16 \cdot N_f N_m u - N_m - N_f)} $$

$$ \text{QI:} - \frac{((N_m + N_f)u^4 + (-4 \cdot N_m - 4 \cdot N_f)u^3 + (6 \cdot N_m + 6 \cdot N_f)u^2 + (-4 \cdot N_m - 4 \cdot N_f)u + N_m + N_f)}{((N_m + N_f)u^4 + (-4 \cdot N_m - 4 \cdot N_f)u^3 + ((8 \cdot N_f + 4) \cdot N_m + 4 \cdot N_f)u^2 - 16 \cdot N_f N_m u - N_m - N_f)} $$
$$ (QI) \frac{-(N_m+N_f)u^4-(-4N_m-4N_f)u^3-(6N_m+6N_f)u^2-(-4N_m-4N_f)u-N_m-N_f}{(N_m+N_f)u^4+(-4N_m-4N_f)u^3+((8N_f+4)N_m+4N_f)u^2-16N_fN_mu-N_m-N_f} $$

(%i8)

$$ QS: -((N_m+N_f)u^4+(-2 \cdot N_m-2 \cdot N_f)u+N_m+N_f)/((N_m+N_f)u^4+(-4 \cdot N_m-4 \cdot N_f)u^3+((8 \cdot N_f+4) \cdot N_m+4 \cdot N_f)u^2-16 \cdot N_f N_m u-N_m-N_f) $$

$$ (QS) \frac{-(N_m+N_f)u^2-(-2N_m-2N_f)u-N_m-N_f}{(N_m+N_f)u^4+(-4N_m-4N_f)u^3+((8N_f+4)N_m+4N_f)u^2-16N_fN_mu-N_m-N_f} $$

(%i9) FIS: (QI-QS)/(1-QS);

(FIS)

$$ \frac{-(N_m+N_f)u^4-(-4 \times N_m-4 \times N_f)u^3-(6 \times N_m+6 \times N_f)u^2-(-4 \times N_m-4 \times N_f)u-N_m-N_f}{(N_m+N_f)u^4+(-4 \times N_m-4 \times N_f)u^3+((8 \times N_f+4) \times N_m+4 \times N_f)u^2-16 \times N_f \times N_m u-N_m-N_f} \frac{-(N_m+N_f)u^2-(-2 \times N_m-2 \times N_f)u-N_m-N_f}{1-\frac{-(N_m+N_f)u^2-(-2 \times N_m-2 \times N_f)u-N_m-N_f}{(N_m+N_f)u^4+(-4 \times N_m-4 \times N_f)u^3+((8 \times N_f+4) \times N_m+4 \times N_f)u^2-16 \times N_f \times N_m u-N_m-N_f}} $$

(%i10) FIS2: ratsimp(FIS);

$$ (FIS2) \frac{(N_m+N_f)u^2+(-2N_m-2N_f)u+N_m+N_f}{(N_m+N_f)u^4+(-2N_m-2N_f)u+(8N_f+1)N_m+N_f} $$

## Script 4: Computing eigenpairs in a dioecious population (even sex-ratio)

(%i1) A: matrix([0,1], [1/(2·N),1-1/N]);

$$ (A) \begin{pmatrix} 0 & 1 \\ \frac{1}{2N} & 1-\frac{1}{N} \end{pmatrix} $$

(%i2) eigenvectors(A);

$$ (%i2) \text{eigenvectors}(A); \text{[}[[[-\frac{\sqrt{N^2+1}-N+1}{2N}, \frac{\sqrt{N^2+1}-N+1}{2N}], [1,1]], \text{[[}[1, -\frac{\sqrt{N^2+1}-N+1}{2N}]\text{]]}, \text{[[}[1, \frac{\sqrt{N^2+1}-N-1}{2N}]\text{]]}]] $$

## Script 5: Computing partial Q in a dioecious population (even sex-ratio)

(%i7) D: matrix([\lambda_1^1(t-1),0], [0,\lambda_2^1(t-1)]);

$$ (D) \begin{pmatrix} \lambda_1^{t-1} & 0 \\ 0 & \lambda_2^{t-1} \end{pmatrix} $$

(%i8) P: matrix([1,1], [\lambda_1,\lambda_2]);

$$ (P) \begin{pmatrix} 1 & 1 \\ \lambda_1 & \lambda_2 \end{pmatrix} $$

--> ;

--> Q_1: matrix([Q_I_1], [Q_S_1]);

$$ (Q_1) \begin{pmatrix} Q_{I_1} \\ Q_{S_1} \end{pmatrix} $$

(%i10) Q_partial: v^1(t-1)·P.D invert(P).Q_1;

$$ (Q_partial) \begin{pmatrix} \gamma^{t-1} \left( \lambda_1^{t-1} \left( \frac{Q_{I_1} \lambda_2}{\lambda_2 - \lambda_1} - \frac{Q_{S_1}}{\lambda_2 - \lambda_1} \right) + \lambda_2^{t-1} \left( \frac{Q_{S_1}}{\lambda_2 - \lambda_1} - \frac{Q_{I_1} \lambda_1}{\lambda_2 - \lambda_1} \right) \right) \\ \gamma^{t-1} \left( \lambda_1^t \left( \frac{Q_{I_1} \lambda_2}{\lambda_2 - \lambda_1} - \frac{Q_{S_1}}{\lambda_2 - \lambda_1} \right) + \lambda_2^t \left( \frac{Q_{S_1}}{\lambda_2 - \lambda_1} - \frac{Q_{I_1} \lambda_1}{\lambda_2 - \lambda_1} \right) \right) \end{pmatrix} $$