# Demographic inference using genetic data from a single individual: separating population size variation from population structure

Olivier Mazet, Willy Rodríguez, Lounès Chikhi

► To cite this version:

Olivier Mazet, Willy Rodríguez, Lounès Chikhi. Demographic inference using genetic data from a single individual: separating population size variation from population structure. Theoretical Population Biology, 2015, 104, pp.46 - 58. 10.1016/j.tpb.2015.06.003. hal-01087849

HAL Id: hal-01087849

https://hal.science/hal-01087849

Submitted on 27 Nov 2014

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>
# Supplementary Materials

Mazet Olivier, Rodríguez Valcarce Willy, Chikhi Lounès

## 1 Proof of the lemma

**Lemma 1.1** Given set of $n_L$ independent observations $(t_1, t_2, ..., t_{n_L})$, the log-likelihood function $\log(\mathbb{L}_{SSPSC})$ has no critical points in $\mathbb{R}^2$.

**Proof** Given $n_L$ independent values $t = (t_1, t_2, ..., t_{n_L})$, the likelihood is:

$$ \mathbb{L}_{SSPSC}(\alpha, T) = \prod_{i=1}^{n_L} \mathbb{L}_{t_i}(T, \alpha), $$

and taking the log:

$$ \log(\mathbb{L}_{SSPSC}(\alpha, T)) = \sum_{i=1}^{n_L} \log(\mathbb{L}_{t_i}(\alpha, T)), \quad (1) $$

where

$$ \mathbb{L}_{t_i}(T, \alpha) = \frac{1}{\alpha} e^{-T - \frac{1}{\alpha}(t_i - T)} \mathbb{I}_{T \le t_i} + e^{-t_i} \mathbb{I}_{T > t_i} \quad (2) $$

First note that:

* For $T > t_i$, $\mathbb{L}_{t_i}(\alpha, T) = e^{-t_i}$ is **constant** (with respect to $(\alpha, T)$).

* If $\alpha \neq 1$ ($\alpha = 1$ means there is no change in the population's size) then $\mathbb{L}_{t_i}(T, \alpha)$ has a discontinuity at $T = t_i$.

As we are interested in the case $\alpha \neq 1$, the log-likelihood function has discontinuities at each $t_i$, $i = 1... n_L$.

For $i \in \{0, 1, ..., n_L + 1\}$, let $C_i = \{(\alpha, T) \in \mathbb{R}^+ \times \mathbb{R}^+, t_i < T < t_{i+1}\}$ with $t_0 = 0$ and $t_{n_L+1} = +\infty$.

Now let be $C = \bigcup_{i=0}^{n_L} C_i$ (Figure 1). We can see that $\mathbb{L}_{SSPSC}(\alpha, T)$ is continuously differentiable in the interior of $C$.
Demographic inference using genetic data from a single individual:
separating population size variation from population structure

Mazet Olivier*, Rodríguez Valcarce Willy*, and Chikhi Lounès§,†,‡

*Université de Toulouse, Institut National des Sciences Appliquées, Institut de Mathématiques de Toulouse, F-31077 Toulouse, France

§CNRS, Université Paul Sabatier, ENFA, UMR 5174 EDB (Laboratoire Évolution & Diversité Biologique), F-31062 Toulouse, France

†Université de Toulouse, UPS, EDB, F-31062 Toulouse, France

‡Instituto Gulbenkian de Ciência, P-2780-156 Oeiras, Portugal
<figure><img src="image_1.png" /><figcaption>Figure 1: Cutting the parameter space, taking into account the discontinuities of the log-likelihood function in $\mathbb{R}^2$</figcaption></figure>
Running Head: Structure vs Population Size Change

Key Words: symmetric island model, population size change, maximum likelihood estimation, demographic history, coalescence time

Corresponding Author:

Chikhi Lounès

CNRS, Université Paul Sabatier, Laboratoire Evolution & Diversité Biologique, Bâtiment
4R1, 118 route de Narbonne, 31062 Toulouse cedex 9, France.

lounes.chikhi@univ-tlse3.fr
Given that we do not consider negative values for $\alpha$ or $T$, we split the parameter space into two subsets $C$ and $\mathbb{R}_+^2 \setminus C$. If $(\alpha, T) \in C$, taking the log and the derivative with respect to $T$ in equation (2) gives:

$$ \frac{\partial}{\partial T} \mathbb{L}_{t_i}(\alpha, T) = \begin{cases} -1 + \frac{1}{\alpha} & \text{if } T < t_i \\ 0 & \text{otherwise.} \end{cases} $$

As we can see, if $\alpha < 1$ then $\frac{\partial}{\partial T} \mathbb{L}_{t_i}(\alpha, T) > 0$ for all $i$ and if $\alpha > 1$ then $\frac{\partial}{\partial T} \mathbb{L}_{t_i}(\alpha, T) < 0$ for all $i$. A consequence, $\nabla \mathbb{L}_{SSPSC}(\alpha, T)$ will never be zero in the interior of $C$ if $\alpha \neq 1$. $\square$

This fact suggests that the min and max values of $\mathbb{L}_{SSPSC}$ (if they exist) have the form $(\alpha, t_i)$.

Let's find the critical points of $\mathbb{L}_{SSPSC}$ over the lines $(\alpha, t_i)$ with $t_i \in \{t_1, t_2, ..., t_{n_L}\}$. When we fix the value of $T$ the function becomes a function of the single variable $\alpha$. If $T = t_a$ for $a \in \{1, ..., n_L\}$ it follows from (1) that:

$$ \log(\mathbb{L}_{SSPSC}(\alpha, t_a)) = \log\left(\frac{1}{\alpha}\right) - t_a - \sum_{i=1}^{n_L} \left[\log\left(\frac{1}{\alpha}\right) - t_a - \frac{1}{\alpha}(t_i - t_a)\right] \mathbb{I}_{t_a < t_i} - \sum_{i=1}^{n_L} t_i \mathbb{I}_{t_a > t_i} $$

Denoting $K = \sum_{i=1}^{n_L} \mathbb{I}_{t_i < t_a}$, we then have:

$$ \log(\mathbb{L}_{SSPSC}(\alpha, t_a)) = (K+1) \left( \log\left(\frac{1}{\alpha}\right) - t_a \right) - \frac{1}{\alpha} \sum_{i=1}^{n_L} (t_i - t_a) \mathbb{I}_{t_a < t_i} - \sum_{i=1}^{n_L} t_i \mathbb{I}_{t_a > t_i} $$

Let us find the zeros of the derivative in $\alpha$:

$$ \begin{align*} \log(\mathbb{L}_{SSPSC}(\alpha, t_a)) = 0 &\Leftrightarrow -(K+1)\frac{1}{\alpha} + \frac{1}{\alpha^2}\sum_{i=1}^{n_L}(t_i - t_a)\mathbb{I}_{t_a < t_i} = 0 \\ &\Leftrightarrow -\alpha(K+1) + \sum_{i=1}^{n_L}(t_i - t_a)\mathbb{I}_{t_a < t_i} = 0 \\ &\Leftrightarrow \frac{\sum_{i=1}^{n_L} t_i \mathbb{I}_{t_a < t_i} - Kt_a}{K+1} = 0 \end{align*} $$

Hence, the maximum value of the log-likelihood function (if it exists) is of the form:
Abstract

The rapid development of sequencing technologies represents new opportunities for pop-
ulation genetics research. It is expected that genomic data will increase our ability to re-
construct the history of populations. While this increase in genetic information will likely
help biologists and anthropologists to reconstruct the demographic history of populations,
it also represents new challenges. Recent work has shown that structured populations gen-
erate signals of population size change. As a consequence it is often difficult to determine
whether demographic events such as expansions or contractions (bottlenecks) inferred from
genetic data are real or due to the fact that populations are structured in nature. Given
that few inferential methods allow us to account for that structure, and that genomic data
will necessarily increase the precision of parameter estimates, it is important to develop new
approaches. In the present study we analyse two demographic models. The first is a model
of instantaneous population size change whereas the second is the classical symmetric island
model. We (i) re-derive the distribution of coalescence times under the two models for a sam-
ple of size two, (ii) use a maximum likelihood approach to estimate the parameters of these
models (iii) validate this estimation procedure under a wide array of parameter combina-
tions, (iv) implement and validate a model choice procedure by using a Kolmogorov-Smirnov
test. Altogether we show that it is possible to estimate parameters under several models and
perform efficient model choice using genetic data from a single diploid individual.
$$m_a = \left( \frac{\sum_{i=1}^{n_L} t_i \mathbb{I}_{t_a < t_i} - Kt_a}{K+1}, t_a \right), a \in \{1, 2, \dots, n_L\}.$$

We then take $(\hat{\alpha}, \hat{T}) = \operatorname{argmax}_{a \in \{1, \dots, n_L\}} \{\log(\mathbb{L}_{SSPSC}(m_a))\}$ as the Maximum Likelihood Estimation.

**Remark:** Note that the function $\mathbb{L}_{SSPSC}$ does not have an upper bound. Let's say $(T_1, T_2, \dots, T_{n_L})$ are the $n_L$ observations of $T_2$ sorted from the lower value to the higher value.

For $T = T_{n_L}$ we have from (1) and (2):

$$\mathbb{L}_{SSPSC}(\alpha, T_{n_L}) = \log\left(\frac{1}{\alpha}\right) - T_{n_L} - \sum_{i=1}^{n_L-1} T_i$$

which clearly goes to $+\infty$ as $\alpha$ goes to zero. The same problem is present for every $T \in [T_{n_L-1}, T_{n_L}]$

Nevertheless, in practice, our $m_a$ is a good estimate for the parameters.
20

INTRODUCTION

The sheer amount of genomic data that is becoming available for many organisms with the rapid development of sequencing technologies represents new opportunities for population genetics research. It is hoped that genomic data will increase our ability to reconstruct the history of populations (LI and DURBIN 2011) and detect, identify and quantify selection (VITTI et al. 2013). While this increase in genetic information will likely help biologists and anthropologists to reconstruct the demographic history of populations, it also exposes old challenges in the field of population genetics. In particular, it becomes increasingly necessary to understand how genetic data observed in present-day populations are influenced by a variety of factors such as population size changes, population structure and gene flow (NIELSEN and BEAUMONT 2009). Indeed, the use of genomic data does not lead to an improvement of statistical inference. If the model assumed to make statistical inference is fundamentally mis-specified, then increasing the amount of data will lead to increased precision for perhaps misleading if not meaningless parameters and will not reveal new insights (NIELSEN and BEAUMONT 2009; CHIKHI et al. 2010; HELLER et al. 2013).

For instance, several recent studies have shown that the genealogy of genes sampled from a deme in an island model is similar to that of genes sampled from a non structured isolated population submitted to a demographic bottleneck (PETER et al. 2010; CHIKHI et al. 2010; HELLER et al. 2013). As a consequence, using a model of population size change for a spatially structured population may falsely lead to the inference of major population size changes (STÄDLER et al. 2009; PETER et al. 2010; CHIKHI et al. 2010; HELLER et al. 2013; PAZ-VINAS et al. 2013). Conversely, assuming a structured model to estimate rates of gene flow when a population has been submitted to a population size change, may also generate misleading conclusions, even though the latter case has been much less documented. More generally, previous studies have shown that spatial processes can mimic selection (CURRAT et al. 2006), population size changes (CHIKHI et al. 2010; HELLER et al. 2013) or that changes in gene flow patterns can mimic changes in population size (WAKELEY 1999; BROQUET et al.
## 2 Proportion of rejections

<figure><img src="image_2.png" /><figcaption>Figure 2: Proportion of rejected data sets simulated under the SSPSC model. Left panels: the reference model is the StSI model (the wrong model). Right panels: the reference model is the SSPSC (the right model), i.e. the model under which the data were simulated. Note that for the abscissa we used $2n_L$ instead of $n_L$ because in order to perform the KS test it is necessary to first estimate the parameters using $n_L$ loci and then an independent set of $n_L$ values of $T_2$.</figcaption></figure>
2010). The fact that such dissimilar processes can generate similar coalescent trees poses exciting challenges (NIELSEN and BEAUMONT 2009). One key issue here is that it may be crucial to identify the kind of model (or family of models) that should be used before estimating and interpreting parameters.

One solution to this problem is to identify the "best" model among a set of competing models. This research program has been facilitated by the development of approximate Bayesian computation (ABC) methods (PRITCHARD et al. 2000; BEAUMONT et al. 2002; CORNUET et al. 2008; BEAUMONT 2010). For instance, using an ABC approach, PETER et al. (2010) showed that data sets produced under population structure can be discriminated from those produced under a population size change by using up to two hundred microsatellite loci genotyped for 25 individual. In some cases, relatively few loci may be sufficient to identify the most likely model (SOUSA et al. 2012; PETER et al. 2010), but in others, tens or hundreds of loci may be necessary (PETER et al. 2010). ABC approaches are thus potentially very powerful but it may still be important to improve our understanding of the coalescent under structured models.

In the present study we are interested in describing the properties of the coalescent under two models of population size change and population structure, respectively, and in devising a new statistical test and estimation procedures. More specifically we re-derive the full distribution of $T_2$, the time to the most recent common ancestor for a sample of size two for a model of sudden population size change and for the *n*-island model. We then use a maximum likelihood approach to estimate the parameters of interest for each model (timing and ratio of population size change former and rate of gene flow and deme size for the latter). We develop a statistical test that identifies data sets generated under the two models. Finally, we discuss how these results may apply to genomic data and how they could be extended to real data sets (since $T_2$ is not usually known) and other demographic models.

In particular we discuss how our results are relevant in the context of the PSMC (Pairwise Sequentially Markovian Coalescent) method (LI and DURBIN 2011), which has been now
<figure><img src="image_3.png" /><figcaption>Figure 3: Proportion of rejected data sets simulated under the StSI model. Left panels: the reference model is the SSPSC (the wrong model). Right panels: the reference model is the StSI model (the right model), i.e. the model under which the data were simulated. Note that for the abscissa we used $2n_L$ instead of $n_L$ because in order to perform the KS test it is necessary to first estimate the parameters using $n_L$ loci and then an independent set of $n_L$ values of $T_2$.</figcaption></figure>
74 extensively used on genomic data and also uses a sample size of two.
# 3 Estimation accuracy

All the scripts used in the simulations are available from http://github.com/willyrv/SSPSC_vs_StSI

Figures 4, 5, 6, 7, 8 and 9 show the results of the estimation of $\alpha$ for all the values of $T$ and all the sample sizes, under the SSPSC model. The values where $\alpha = (2, 4, 10, 20, 50, 100)$, $T = (0.1, 0.2, 0.5, 1, 2, 5)$ and $n_L = (20, 50, 100, 200, 500, 1000, 10000)$

<figure><img src="image_4.png" /><figcaption>Figure 4: Estimations of $\alpha$ (real value $\alpha = 2.0$)</figcaption></figure>
76 Demographic models

77 Population size change: We consider a simple model of population size change, where N(t) represents the population size (N, in units of genes or haploid genomes) as a function of time (t) expressed in generations scaled by N, the population size, and where t = 0 is the present, and positive values represent the past (Figure 1 (a)). More specifically we assume a sudden change in population size at time T in the past, where N changes instantaneously by a factor α. This can be summarized as N(t) = N(0) = N₀ for t ∈ [0, T[, N(t) = N(T) = αN₀ for t ∈ [T, +∞[. If α > 1 the population went through a bottleneck (Figure 1) whereas if α < 1 it expanded. Since N represents the population size in terms of haploid genomes, the number of individuals will therefore be N/2 for diploid species. Note also that for a population of constant size the expected coalescence time of two genes is N generations, which therefore corresponds to t = 1. We call this model the SSPSC, which stands for Single Step Population Size Change.

[Figure 1 about here.]

90 *Structured population:* Here we consider the classical symmetric *n*-island model (WRIGHT 1931), see Figure 1 (b), where we have a set of *n* islands (or demes) of constant size *N*, interconnected by gene flow with a migration rate *m*, where *M* = *Nm* is the number of immigrants (genes) in each island every generation. The whole metapopulation size is therefore *nN* (this is the total number of genes, not the effective size). Again, *N* is the number of haploid genomes, and *N*/2 the number of diploid individuals. We call this model the StSI, which stands for Structured Symmetrical Island model.
<figure><img src="image_5.png" /><figcaption>Figure 5: Estimations of α (real value α = 4.0)</figcaption></figure>
97 The distribution of coalescence times: qualitative and quantitative analyses

In this section we used previous results (HERBOTS 1994; DONNELLY and TAVARE 1995) to derive the distribution of coalescent times for the two models of interest. We show that even though they are different, these distributions can be similar under an indefinitely large number of parameter values (Figures 2 and 3). Moreover we show that even when the distributions are distinguishable, their first moments may not be. In particular, we show that the first two moments (mean and variance) are near identical for a large number of parameter combinations. Before doing that we start by providing a simple intuitive rationale explaining why and how a model of population structure can be mistaken for a model of population size change. This intuitive approach is important because it allows us to understand how the parameters of the two models ((T, α) and (M, n), respectively) are linked.

*Intuitive and qualitative rationale:* We start by taking two genes sampled in the present-day population under the Single Step Population Size Change (SSPSC) model. If we assume that $\alpha > 1$ (population bottleneck from an ancient population of size $N_1$ to a current population of size $N_0$, with $N_1 = \alpha N_0$) the probability that the two genes coalesce will vary with time as a function of $N_0$, $N_1$ and $T$. If $T$ is very small, then most genes will coalesce at a rate determined by $N_1$, whereas if $T$ is very large the coalescence rate will be mostly determined by $N_0$. If we now take two genes sampled from the same island in the Structured Symmetrical Island (StSI) model, we can also see that their coalescence rate will depend on $N$, the size of the island and on $m$, the migration rate. If $m$ is very low, the coalescence rate should mostly depend on $N$. If $m$ is high, the two genes may see their lineages in different islands before they coalesce. As a consequence the coalescence rate will depend on the whole set of islands and therefore on the product $nN$, where $n$ is the total number of islands.

This intuitive description suggests that there is an intrinsic relationship between $T$ and $1/M$, and between $\alpha$ and $n$. The reason why structured populations exhibit signals of
<figure><img src="image_6.png" /><figcaption>Figure 6: Estimations of α (real value α = 10.0)</figcaption></figure>
bottlenecks is because in the recent past the coalescence rate depends on the local island size $N$, whereas in a more distant past it depends on $nN$. In other words, it is as if the population size had been reduced by a factor of $n$. As we will see this rationale is only qualitatively correct, but it suggests that if we want to distinguish them it may be necessary to derive the full distribution of the coalescence times under the two models. We shall denote these coalescence times $T_2^{SSPSC}$ and $T_2^{StSI}$, respectively.

### Derivation of the distribution of coalescence times:

**The distribution of $T_2^{SSPSC}$:** The generalisation of the coalescent in populations of variable size was first rigorously treated in DONNELLY and TAVARE (1995)), and is clearly exposed in TAVARÉ (2004)). If we denote by $\lambda(t)$ the ratio $N(t)/N(0)$ where $t$ is the time scaled by the number of genes (i.e. units of coalescence time, corresponding to $\lfloor N(0)t \rfloor$ generations), we can compute the probability density function (pdf) $f_{T_2}^{SSPSC}(t)$ of the coalescence time $T_2^{SSPSC}$ of two genes sampled in the present-day population. Indeed, the probability that two genes will coalesce at a time greater than $t$ is

$$ \mathbb{P}(T_2^{SSPSC} > t) = e^{-\int_0^t \frac{1}{\lambda(x)} dx}, \qquad (1) $$

where

$$ \lambda(x) = \mathbb{I}_{[0,T]}(x) + \alpha \mathbb{I}_{[T, +\infty]}(x), $$

and $\mathbb{I}_{[a,b]}(x)$ is the Kronecker index such that

$$ \mathbb{I}_{[a,b]}(x) = \begin{cases} 1 & \text{for } x \in [a, b[ \\ 0 & \text{otherwise.} \end{cases} $$

Given that the pdf is

$$ f_{T_2}^{SSPSC}(t) = (1 - \mathbb{P}(T_2^{SSPSC} > t))' $$
<figure><img src="image_7.png" /><figcaption>Figure 7: Estimations of α (real value α = 20.0)</figcaption></figure>
Equation (1) can be rewritten as

$$ \mathbb{P}(T_2^{SSPSC} > t) = e^{-t} \mathbb{I}_{[0,T]} + e^{-T-\frac{1}{\alpha}(t-T)} \mathbb{I}_{[T, +\infty[} $$

This leads to the following pdf

$$ f_{T_2}^{SSPSC}(t) = e^{-t} \mathbb{I}_{[0,T]}(t) + \frac{1}{\alpha} e^{-T-\frac{1}{\alpha}(t-T)} \mathbb{I}_{[T, +\infty[}(t). \quad (2) $$

**The distribution of $T_2^{StSI}$:** Following HERBOTS (1994), an easy way to derive the distribution of the coalescence time $T_2^{StSI}$ of two genes for our structured model, is to compute the probability that two genes are identical by descent when they are sampled from the same or from different populations. These two probabilities are respectively denoted by $p_s(\theta)$ and $p_d(\theta)$, where $\theta = 2uN$ is the scaled mutation rate, $u$ being the per locus mutation rate.

Indeed, using a classical scaling argument (see for instance TAVARÉ (2004), page 34), we can note that

$$ p_s(\theta) = \mathbb{E}(e^{-\theta T_2^{StSI}}) $$

In other words $p_s(\theta)$ is the Laplace transform of $T_2^{StSI}$.

We can compute this probability as follows. Taking two genes from the same island and going back in time, there are three events that may occur: a coalescence event (with rate 1), a mutation event (with rate $\theta$) and a migration event (with rate $M$). Taking now two genes from different islands, they cannot coalesce and therefore only a mutation or a migration event may occur. Migration events can then bring the lineages in the same island with probability $\frac{1}{n-1}$, and in different islands with probability $\frac{n-2}{n-1}$. We thus obtain the following coupled equations:

$$ p_s(\theta) = \frac{1}{1+M+\theta} + \frac{M}{1+M+\theta}p_d(\theta), $$

and

$$ p_d(\theta) = \frac{M/(n-1)}{M+\theta}p_s(\theta) + \frac{M(n-2)/(n-1)}{M+\theta}p_d(\theta). $$
<figure><img src="image_8.png" /><figcaption>Figure 8: Estimations of α (real value α = 50.0)</figcaption></figure>
By solving them, we obtain

$$p_s(\theta) = \frac{\theta + \gamma}{D} \quad \text{and} \quad p_d(\theta) = \frac{\gamma}{D}$$

with

$$\gamma = \frac{M}{n-1} \quad \text{and} \quad D = \theta^2 + \theta(1 + n\gamma) + \gamma.$$

We can then obtain the full distribution through the Laplace transform formula, if we note that

$$p_s(\theta) = \frac{\theta + \gamma}{(\theta + \alpha)(\theta + \beta)} = \frac{a}{\theta + \alpha} + \frac{1 - a}{\theta + \beta}$$

with

$$a = \frac{\gamma - \alpha}{\beta - \alpha} = \frac{1}{2} + \frac{1 + (n - 2)\gamma}{2\sqrt{\Delta}},$$

150 where

$$\alpha = \frac{1}{2} \left( 1 + n\gamma + \sqrt{\Delta} \right)$$

and

$$\beta = \frac{1}{2} \left( 1 + n\gamma - \sqrt{\Delta} \right).$$

Noting now that for any $\theta$ and any $\alpha$ we have

$$\int_{0}^{+\infty} e^{-\alpha s} e^{-\theta s} ds = \frac{1}{\theta + \alpha},$$

151 it is straightforward to see that the pdf of $T_2^{StSI}$ is an exponential mixture:

$$f_{T_2}^{StSI}(t) = ae^{-\alpha t} + (1-a)e^{-\beta t}. \qquad (3)$$

152 First moments: Equations 2 and 3 are different hence showing that it is in principle possible to identify genetic data produced under the two demographic models of interest. The two 153 equations can be used to derive the expectation and variance of the two random variables 154 of interest, $T_2^{SSPSC}$ and $T_2^{StSI}$. Their analytic values can be easily expressed as functions of 155 the model parameters: 156
<figure><img src="image_9.png" /><figcaption>Figure 9: Estimations of α (real value α = 100.0)</figcaption></figure>
$$
\begin{align*}
\mathbb{E}(T_2^{\text{SSPSC}}) &= 1 + e^{-T}(\alpha - 1), \\
\text{Var}(T_2^{\text{SSPSC}}) &= 1 + 2Te^{-T}(\alpha - 1) + 2\alpha e^{-T}(\alpha - 1) - (\alpha - 1)^2 e^{-2T}, \\
\mathbb{E}(T_2^{\text{StSI}}) &= n, \\
\text{Var}(T_2^{\text{StSI}}) &= n^2 + \frac{2(n-1)^2}{M}.
\end{align*}
$$

It is interesting to note that the expected time in the StSI model is $n$ and does not depend on the migration rate (DURRETT 2008). The variance is however, and as expected, a function of both $n$ and $M$. For the SSPSC model, the expected coalescence time is a function of both $T$ and $\alpha$. We note that it is close to 1 when $T$ is very large and to $\alpha$ when $T$ is close to zero. Indeed, when the population size change is very ancient, even if $\alpha$ is very large the expected coalescence time will mostly depend on the present-day population size, $N_0$. Similarly, when $T$ is small it will mostly depend on $N_1$. The relationship that we mentioned above between $n$ and $\alpha$ (and between $M$ and $1/T$) can be seen by noting that when $T$ is close to zero (and $M$ is large), the expectations under the two models are $\alpha$ and $n$, and the variances are $\text{Var}(T_2^{\text{SSPSC}}) \approx 1 + 2\alpha(\alpha - 1) - (\alpha - 1)^2 = \alpha^2$ and $\text{Var}(T_2^{\text{StSI}}) \approx n^2$. This exemplifies the intuitive rationale presented above. This relationship is approximate and will be explored below, but can be illustrated in more general terms by identifying scenarios with similar moments.

As figure 2 shows, the two models provide near-identical pairs of values for $(\mathbb{E}(T_2), \text{Var}(T_2))$ for “well chosen” parameters $(T, \alpha)$ and $(M, n)$. Here by setting $T$ to 0.1 (and $M$ to 9, i.e. $1/M \approx 0.11$) whereas $\alpha$ and $n$ were allowed to vary from 1 to 100, and from 2 to 100, respectively, we see that the two models exhibit very similar behaviours. We also plotted a second example obtained by setting $M$ to 0.5 and $T$ to 1.09, and varying $n$ and $\alpha$ as above. These examples illustrate how $n$ and $\alpha$ (respectively, $M$ and $1/T$) are intimately related.

[Figure 2 about here.]

The near-identical values obtained for the expectation and variance under the two models
Figures 10, 11, 12, 13, 14 and 15 show the results of the estimation of $T$ for all the values of $\alpha$ and all the sample sizes, under the SSPSC model. The values where $\alpha = (2, 4, 10, 20, 50, 100)$, $T = (0.1, 0.2, 0.5, 1, 2, 5)$ and $n_L = (20, 50, 100, 200, 500, 1000, 10000)$

<figure><img src="image_10.png" /><figcaption>Figure 10: Estimations of $T$ (real value $T = 0.1$)</figcaption></figure>
explains why it may be difficult to separate models of population size change from models of
population structure when the number of independent genetic markers is limited. However,
the differences between the distributions of coalescence times under the two models suggest
that we can probably go further and identify one model from another. For instance, figure
3 shows that even in cases where the first two moments are near-identical ($T = 0.1$ and
$\alpha = 10$ versus $M = 7$ and $n = 9$), it should be theoretically possible to distinguish them.
This is exactly what we aim to do in the next section. In practice, we will assume that
we have a sample of $n_L$ independent $T_2$ values (corresponding to $n_L$ independent loci) and
will use these $T_2$ values to (i) estimate the parameter values that best explain this empirical
distribution under the two models of interest, (ii) use a statistical test to compare the
empirical distribution with the expected distribution for the ML estimates and reject (or
not) one or both of the models. For simplicity, and to make it easier to read, we will often
use the term loci in the rest of the manuscript when we want to mention the number of
independent $T_2$ values.

[Figure 3 about here.]
<figure><img src="image_11.png" /><figcaption>Figure 11: Estimations of T (real value T = 0.2)</figcaption></figure>
193 Model choice and parameter estimation

General principle and parameter combinations: Given a sample (t₁, ..., tₙₐ) of nₐ independent
observations of the random variable T₂, we propose a parameter estimation procedure and
a goodness-of-fit test to determine whether the observed distribution of the T₂ values is
significantly different from that expected from the theoretical T₂<sup>SSPSC</sup> or T₂<sup>StSI</sup> distributions.
This sample can be seen as a set of T₂ values obtained or estimated from nₐ independent loci.
We took a Maximum Likelihood (ML) approach to estimate the parameters (T, α) and (M, n)
under the hypothesis that the nₐ-sample was generated under the T₂<sup>SSPSC</sup> and the T₂<sup>StSI</sup>
distributions, respectively (see Supplementary materials for the details of the estimation
procedure). The ML estimates (T̂, α̂) and (M̂, n̂) were then used to define T₂<sup>SSPSC</sup> or T₂<sup>StSI</sup>
reference distributions. The Kolmogorov-Smirnov (KS) test which allows to compare a
sample with a reference distribution was then used to determine whether the observed nₐ
sample could have been generated by the respective demographic models. In other words
this allowed us to reject (or not) the hypothesis that the (t₁, ..., tₙₐ) sample was a realization
of the reference distributions (T₂<sup>StSI</sup> or T₂<sup>SSPSC</sup>). Note that the estimation procedure and
the KS test were performed on independent sets of T₂ values. We thus simulated twice as
many T₂ values as needed (2nₐ instead of nₐ). With real data that would require that half
of the loci be used to estimate (T̂, α̂) and (M̂, n̂), whereas the other half would be used to
perform the KS test.

We expect that if the estimation procedure is accurate and if the *KS* test is performing
well we should reject the SSPSC (respectively, the StSI) model when the data were simu-
lated under the StSI (resp., the SSPSC) model. On the contrary we should not reject data
simulated under the SSPSC (resp., the StSI) model when they were indeed simulated under
that model. To validate our approach we used (t₁, ..., t₂nₐ) data sampled from the two T₂
distributions and quantified how the estimation procedure and the *KS* test performed. In
order to do that, we varied the parameter values ((T, α) and (M, n)) for various 2nₐ values
<figure><img src="image_12.png" /><figcaption>Figure 12: Estimations of T (real value T = 0.5)</figcaption></figure>
as follows. For $T$ and $\alpha$ we used all 36 pairwise combinations between these two sets of values (0.1, 0.2, 0.5, 1, 2, 5), and (2, 4, 10, 20, 50, 100), respectively. For $M$ and $n$ we used all the 48 combinations between the following values (0.1, 0.2, 0.5, 1, 5, 10, 20, 50) and (2, 4, 10, 20, 50, 100), respectively. For $2n_L$ we used the following values (40, 100, 200, 400, 1000, 2000, 20000). Altogether we tested 588 combinations of parameters and number of loci. For each $2n_L$ value and for each parameter combination $(T, \alpha)$ (or $(M, n)$) we realized 100 independent repetitions of the following process. We first simulated a sample of $2n_L$ values using the pdfs of the SSPSC (resp. StSI) model with $(T, \alpha)$ (resp. $(M, n)$). We then used the first $n_L$ values to obtain the ML estimates $(\hat{T}, \hat{\alpha})$ for the SSPSC model and $(\hat{M}, \hat{n})$ for the StSI model. Then, we performed a KS test using a 0.05 threshold on the second half of the simulated data (i.e. $n_L$ values) with each of the theoretical distributions defined by the estimated parameters. Finally, after having repeated this process 100 times we recorded all estimated parameters and counted the number of times we rejected the SSPSC and StSI models for each parameter combination and each $2n_L$ value.

**Maximum likelihood estimation in the SSPSC case:** We know from section the pdf of the coalescence time in the SSPSC model of two genes. We can thus write the likelihood function for any couple of parameters $(\alpha, T)$, given one observation $t_i$ as:

$$ \mathbb{L}_{t_i}(\alpha, T) = \frac{1}{\alpha} e^{-T - \frac{1}{\alpha}(t_i - T)} \mathbb{I}_{[0, t_i]}(T) + e^{-t_i} \mathbb{I}_{[t_i, +\infty]}(T). $$

Given $n_L$ independent values $t = (t_1, t_2, ..., t_{n_L})$, the likelihood is:

$$ \mathbb{L}_{SSPSC}(\alpha, T) = \prod_{i=1}^{n_L} \mathbb{L}_{t_i}(T, \alpha), $$

and taking the log it gives:

$$ \log(\mathbb{L}_{SSPSC}(\alpha, T)) = \sum_{i=1}^{n_L} \log(\mathbb{L}_{t_i}(\alpha, T)). $$
<figure><img src="image_13.png" /><figcaption>Figure 13: Estimations of T (real value T = 1.0)</figcaption></figure>
**Lemma 0.1** Given a set of $n_L$ independent observations $\{t_1, t_2, ..., t_{n_L}\}$, the log-likelihood function $\log(\mathbb{L}_{SSPSC})$ has no critical points in $\mathbb{R}^2$.

For the proof and some comments, see Supplementary Materials.

As a consequence of this lemma, we take $(\hat{\alpha}, \hat{T}) = \operatorname*{argmax}_{a \in \{1, ..., n_L\}} \{\log(\mathbb{L}_{SSPSC}(m_a))\}$ as the Maximum Likelihood estimates, where

$$m_a = \left( \frac{\sum_{i=1}^{n_L} t_i \mathbb{I}_{t_a < t_i} - Kt_a}{K+1}, t_a \right), a \in \{1, 2, ..., n_L\}.$$

with

$$K = \sum_{i=1}^{n_L} \mathbb{I}_{t_i < t_a}$$

*Maximum likelihood estimation in the StSI case:* Under the StSI model the expression of the critical points is not analytically derived. We know from section the *pdf* of coalescence times for two genes. Given $n_L$ independent values $t = (t_1, t_2, ..., t_{n_L})$ we can compute the log-likelihood function for any set of parameters $(n, M)$ as:

$$\log(\mathbb{L}_{StSI}(n, M)) = \sum_{i=1}^{n_L} \log(ae^{-\alpha t} + (1-a)e^{-\beta t})$$

We used the Nelder-Mead method (NELDER and MEAD 1965) implementation of *scipy* (JONES et al. 2001) to find numerically an approximation to the maximum of the likelihood function. This method returns a pair of real numbers $(\hat{n}, \hat{M})$. Since $n$ should be an integer we kept either $\lfloor \hat{n} \rfloor$ or $\lfloor \hat{n} \rfloor + 1$, depending on which had the largest log-likelihood value.
<figure><img src="image_14.png" /><figcaption>Figure 14: Estimations of T (real value T = 2.0)</figcaption></figure>
Figure 4 shows, for various values of $n_L$, the results of the estimation of $\alpha$ (panels (a), (c), and (e)), for simulations assuming $\alpha = 10$ and $T = (0.1, 1, 2)$, respectively; see Supplementary Material for the other values) and the estimation of $n$ (panels (b), (d), and (f)) for simulations with $n = 10$ and $M = (10, 1, 0.5)$, respectively; see Supplementary Material for the other values, corresponding to 26 figures and 168 panels). The first thing to notice is that both $\alpha$ and $n$ are increasingly well estimated as $n_L$ increases. This is what we expect since $n_L$ represents the amount of information (the number of $T_2$ values or independent loci.) The second thing to note is that the two parameters are very well estimated when we use 10,000 values of $T_2$. This is particularly obvious for $n$ compared to $\alpha$, probably because $n$ must be an integer, whereas $\alpha$ is allowed to vary continuously. For instance, for most simulations we find the exact $n$ value (without error) as soon as we have more than 1000 loci. However, we should be careful in drawing very general rules. Indeed, when fewer $T_2$ values are available (i.e. fewer independent loci), the estimation precision of both parameters depends also on $T$ and $M$, respectively. Interestingly, the estimation of $\alpha$ and $n$ are remarkable even when these parameters are small. This means that even “mild” bottlenecks may be very well quantified (see for instance the Supplementary materials for $\alpha = 2$, $T$ values between 0.1 and 1 when we use only 1000 loci). We should also note that when the bottleneck is very old ($T = 5$) the estimation of the parameters is rather poor and only starts to be reasonable and unbiased for $n_L = 10,000$. This is not surprising since the expected $T_{MRCA}$ is 1. Under the SSPSC model most genes will have coalesced by $t = 5$, and should therefore exhibit $T_2$ values sampled from a stationary population (i.e. $\alpha = 1$). As the number of loci increases, a small proportion will not have coalesced yet and will then provide information on $\alpha$. The expected proportion of genes that have coalesced by $T = 5$ is 0.993.

Figure 5 shows for various values of $n_L$ the results of the estimation of $T$ (panels (a), (c), and (e)), for simulations assuming $T = 0.2$ and $\alpha = (2, 20, 100)$, respectively; see Supplementary Material for the other values) and the estimation of $M$ (panels (b), (d), and (f)), for
<figure><img src="image_15.png" /><figcaption>Figure 15: Estimations of T (real value T = 5.0)</figcaption></figure>
simulations with $M = 20$ and $n = (2, 20, 100)$, respectively; see Supplementary Material for the other values). As expected again, the estimates are getting better as $n_L$ increases. For the values shown here we can see that $T$, the age of the bottleneck is very well estimated even when $\alpha = 2$ (for $n_L = 10,000$). In other words, even a limited bottleneck can be very precisely dated. For stronger bottlenecks fewer loci (between 500 and 1000) are needed to still reach a high precision. This is particularly striking given that studies suggest that it is hard to identify bottlenecks with low $\alpha$ values (GIROD et al. 2011). Interestingly, the panels (b), (d) and (f) seem to suggest that it may be more difficult to estimate $M$ than $T$. As we noted above this observation should be taken with care. Indeed, $T$ and $M$ are not equivalent in the same way as $\alpha$ and $n$. This is why we chose to represent a value of $M$ such that $M = 1/T$, and why one should be cautious in drawing general conclusions here. Altogether this and the previous figure show that it is possible to estimate with a high precision the parameters of the two models by using only 500 or 1000 loci from a single diploid individual. There are also parameter combinations for which much fewer loci could be sufficient (between 50 and 100).

In Figure 6 we show some results of the KS test for the two cases (See the Supplementary Materials for the other parameter combinations). In the left-hand panels ((a), (c), and (e)) the data were simulated under the SSPSC model and we used the StSI model as a reference (i.e. we ask whether we can reject the hypothesis that genetic data were generated under a structured model when they were actually generated under a model of population size change). In the right-hand panels ((b), (d) and (f)) the same data were compared using the SSPSC model as reference and we computed how often we rejected them using a 5% rejection threshold. The left-hand panels exhibit several important features. The first is that, with the exception of $T_2 = 5$ we were able to reject the wrong hypothesis in 100% of the cases when we used 10,000 independent $T_2$ values.

This shows that our estimation procedure (as we saw above in figures 4 and 5) and the KS test are very powerful. The second feature is that for $T = 5$, the test performs badly
Figures 16, 17, 18, 19, 21, 22, 23 and 24 show the results of the estimation of M for all the values of n and all the sample sizes, under the StSI model. The values where $M = (0.1, 0.2, 0.5, 1, 5, 10, 20, 50)$, $n = (2, 4, 10, 20, 50, 100)$ and $n_L = (20, 50, 100, 200, 500, 1000, 10000)$

(0.1, 0.2, 0.5, 1, 5, 10, 20, 50) and (2, 4, 10, 20, 50, 100)

<figure><img src="image_16.png" /><figcaption>Figure 16: Estimations of M (real value $M = 0.1$)</figcaption></figure>
whatever the number of independent loci (at least up to 10,000). This is expected since the
expected $T_{MRCA}$ of two genes is $t = 1$, and 99.3% of the loci will have coalesced by $t = 5$.
This means that out of the 10,000, only c.a. 70 loci are actually informative regarding the
pre-bottleneck population size. Another important feature of the left-hand panels is that
the best results are generally obtained for $T = 1, 0.5$ and $2$, whichever the value of $\alpha$. This
is in agreement with GIROD et al. (2011) in that very recent population size changes are
difficult to detect and quantify. The observation is valid for ancient population size changes
as well. The right-hand panels are nearly identical, whichever $\alpha$ value we used (see also
Supplementary Materials), and whichever number of $T_2$ values we use. They all show that
the *KS* test always rejects a rather constant proportion of data sets. This proportion varies
between 3 and 15%, with a global average of 8.9%. Altogether our *KS* test seems to be
conservative. This is expected because for low $n_L$ values the estimation of the parameters
will tend to be poor. Since the *KS* test uses a reference distribution based on the estimated
rather than the true values, it will reject the simulated data more often than the expected
value of 5%.

Figure 7 is similar to Figure 6 but the data were simulated under the StSI model and
the *KS* test was performed first using the SSPSC model as a reference ((a), (c), (e)) and
then using the StSI model as a reference ((b), (d), (f)). The left-hand panels ((a), (c), and
(e)) show results when we ask whether we can reject the hypothesis that genetic data were
generated under a population size change model when they were actually generated under
a model of population structure. In the right-hand panels ((b), (d), and (f)) we computed
how often we rejected the hypothesis that genetic data were generated under the StSI model
when they were indeed generated under that model of population structure. Altogether,
the left-hand panels suggest that the results are generally best when $M = (0.1, 0.2, 1)$, but
that we get very good results for most values of *M* when we have 10,000 loci and can reject
the SSPSC when they were actually generated under the StSI model. The right-hand panels
show, as in Figure 6, that for all the values of $n_L$ and *n* we reject a rather constant proportion
<figure><img src="image_17.png" /><figcaption>Figure 17: Estimations of M (real value M = 0.2)</figcaption></figure>
of data sets (between 5 and 10%). Altogether the two previous figures (figures 6 and 7) show
that it is possible to identify the model under which genetic data were generated by using
genetic data from a single diploid individual.

Figure 8 is divided in four panels showing the relationships between $T$ and $M$ (panels
(b) and (d), for various values of $\alpha$ and $n$) and between $\alpha$ and $n$ (panels (a) and (c),
for various values of $T$ and $M$). In each of the panels we simulated data under a model for
specific parameter values represented on the x-axis, and estimated parameters from the other
model, and represented the estimated value on the y-axis. Since we were interested in the
relationship between parameters (not in the quality of the estimation, see above), we used the
largest $n_L$ value and plotted the average of 100 independent estimation procedures. In panel
(a) we simulated a population size change (SSPSC) for various $T$ values (represented each by
a different symbol) and several values of $\alpha$ on the x-axis. We then plotted the estimated value
of $\hat{n}$ for each case (i.e. when we assume that the data were generated under the StSI model).
We find a striking linear relationship between these two parameters conditional on a fixed $T$
value. For instance, a population bottleneck by a factor 50 that happened $N_0$ generations ago
($T = 1$) is equivalent to a structured population with $\hat{n} \approx 22$ islands (and $\hat{M} \approx 0.71$). Panel
(c) is similar and shows how data simulated under a structured population generates specific
parameters of population bottlenecks. Panels (b) and (d) show the relationship between $T$
and $M$. We have plotted as a reference the curve corresponding to $y = 1/x$. As noted above
and shown on this graph, this relationship is only approximate and depends on the value
of $\alpha$ and $n$. Altogether, this figure exhibits the profound relationships between the model
parameters. They show that the qualitative relationships between $\alpha$ and $n$, and between $T$
and $1/M$ discussed above are indeed real but only correct up to a correcting factor. Still this
allows us to identify profound relationships between population structure and population
size change.

[Figure 4 about here.]
<figure><img src="image_18.png" /><figcaption>Figure 18: Estimations of M (real value M = 0.5)</figcaption></figure>
359

[Figure 5 about here.]

360

[Figure 6 about here.]

361

[Figure 7 about here.]

362

[Figure 8 about here.]
<figure><img src="image_19.png" /><figcaption>Figure 19: Estimations of M (real value M = 1.0)</figcaption></figure>
In this study we have analysed the distribution of coalescence times under two simple demographic models. We have shown that even though these demographic models are strikingly different (Figure 1) there is always a way to find parameter values for which both models will have the same first two moments (Figure 2). We have also shown that there are intrinsic relationships between the parameters of the two models (Figure 8). However, and this is a crucial point, we showed that the distributions were different and could therefore be distinguished. Using these distributions we developed a ML estimation procedure for the parameters of both models ($\hat{T}, \hat{\alpha}$) and ($\hat{M}, \hat{n}$) and showed that the estimates are accurate, given enough genetic markers. Finally, we showed that by applying a simple KS test we were able to identify the model under which specific data sets were generated. In other words, we were able to determine whether a bottleneck signal detected in a particular data set could actually be caused by population structure using genetic data from a single individual. The fact that a single individual provides enough information to estimate demographic parameters is in itself striking (see in particular the landmark paper by LI and DURBIN (2011)), but the fact that one individual (or rather sometimes as few as 500 or 1000 loci from that one individual) potentially provides us with the ability to identify the best of two (or more) models is even more remarkable. The PSMC (pairwise sequentially markovian coalescent) method developed by LI and DURBIN (2011) reconstructs a theoretical demographic history characterized by population size changes, assuming a single non structured population. Our study goes further and shows that it could be possible to test whether the signal identified by the PSMC is due to actual population size changes or to population structure. However, models putting together these two scenarios have been proposed. In WAKELEY (1999), a model considering a structured population who went through a bottleneck in the past is developed, allowing to estimate the *time* and the *ratio* of the bottleneck. Moreover, WAKELEY (1999) discussed the idea that, in structured populations, changes in the migration rate and/or the number of islands (demes) can shape the observed data in the same way that
<figure><img src="image_20.png" /><figcaption>Figure 20: Estimations of M (real value M = 2.0)</figcaption></figure>
effective population size changes do. Hence, we think that our work could be helpful to the
aim of setting these two scenarios apart in order to detect (for example) false bottleneck
signals. Nevertheless, while our study provides several new results, there are still several
important issues that need to be discussed and much progress that can still be made.

T₂ and molecular data

The first thing to note is that we assume, throughout our study, that we have access to the coalescence times $T_2$. In real data sets, this is never the case and the $T_2$ are rarely estimated from molecular data. While this may seem as a limitation, we should note that the recent method of LI and DURBIN (2011), that uses the genome sequence of a single individual to infer the demographic history of the population it was sampled from, actually estimates the distribution of $T_2$ values. In its current implementation the PSMC software does not output this distribution but it could be modified to do it. Note however, that the PSMC should only be able to provide a discretized distribution in the form of a histogram with classes defined by the number of time periods for which population size estimates are computed. In any case, this suggests that it is in theory possible to use the theoretical work of Li and Durbin to generate $T_2$ distributions, which could then be used with our general approach. Moreover, it should be possible to use the theory developed here to compute, conditional of the $T_2$ distribution, the distribution of several measures of molecular polymorphism. For instance, under an infinite site mutation model it is in principle possible to compute the distribution of the number of differences between pairs of non recombining sequences for the two demographic models analysed here. Similarly, assuming a single stepwise mutation model it should be possible to compute the distribution of the number of repeat differences between two alleles conditional on $T_2$. However, we must add here that while it is easy to use these distributions to simulate genetic data (rather than $T_2$ values) it is not straightforward to then use the genetic data to estimate the models' parameters and apply a test to identify
<figure><img src="image_21.png" /><figcaption>Figure 21: Estimations of M (real value M = 5.0)</figcaption></figure>
415
416
417

the model under which the data were simulated. This would probably require a Khi-2 test
since discrete rather continuous distributions would be compared. This is an issue that would
deserve a full and independent study.

418

Demographic models

419
420
421
422
423
424
425
426
427
428
429
430
431

In our study we limited ourselves to two simple but widely used models. It would thus be
important to determine the extent to which our approach could apply to other demographic
models. The n-island or StSI model is a widely used model and it was justified to use it
here. One of its strongest assumptions is that migration is identical between all demes.
This is likely to be problematic for species with limited vagility. In fact, for many species a
model where migration occurs between neighbouring populations such as the stepping-stone
is going to be more likely. At this stage it is unclear whether one could derive analytically
the *pdf* of *T*₂ for a stepping-stone model. The work by HERBOTS (1994) suggests that it
may be possible to compute it numerically by inversing the Laplace transform derived by
this author. This work has not been done to our knowledge and would still need to be done.
Interestingly, this author has also shown that it is in principle possible to derive analytically
the *pdf* of *T*₂ in the case of a two-island model with populations of different sizes. Again,
this still needs to be done.

432
433
434
435
436
437
438
439

The SSPSC model has also been widely used (ROGERS and HARPENDING 1992) and
represents a first step towards using more complex models of population size change. For
instance, the widely used method of BEAUMONT (1999) to detect, date and quantify pop-
ulation size changes (GOOSSENS et al. 2006; QUÉMÉRÉ et al. 2012; SALMONA et al. 2012)
assumes either an exponential or a linear population size change. It should be straightfor-
ward to compute the *pdf* of *T*₂ under these two models because, as we explained above, the
coalescent theory for populations with variable size has been very well studied (DONNELLY
and TAVARE 1995; TAVARÉ 2004) and it is possible to write the *pdf* of *T*₂ for any demo-
<figure><img src="image_22.png" /><figcaption>Figure 22: Estimations of M (real value M = 10.0)</figcaption></figure>
graphic history involving any type of population size changes. To go even further one could
in principle use the history reconstructed by the PSMC (LI and DURBIN 2011) as the ref-
erence model of population size change and compare that particular demographic history to
an n-island model using our results, and our general approach. At the same time, we should
note that for complex models of population size change, including relatively simple ones such
as the exponential model of BEAUMONT (1999), it is not straightforward to compute the
number of differences between non recombining sequences. Significant work would probably
be needed to apply the general approach outlined in our study to specific demographic mod-
els. But we believe that the possibilities opened by our study are rather wide and should
provide our community with new interesting problems to solve in the next few years.

Comparison with previous work and generality of our results

The present work is part of a set of studies aimed at understanding how population
structure can be mistaken for population size change and at determining whether studies
identifying population size change are mistaken or valid (CHIKHI et al. 2010; HELLER et al.
2013; PAZ-VINAS et al. 2013). It is also part of a wider set of studies that have recognised in
the last decade the importance of population structure as potential factor biasing inference
of demographic (LEBLOIS et al. 2006; STÄDLER et al. 2009; PETER et al. 2010; CHIKHI
et al. 2010; HELLER et al. 2013; PAZ-VINAS et al. 2013) or selective processes (CURRAT
et al. 2006). Here we demonstrated that it is indeed possible to separate the SSPSC and
StSI models. While we believe that it is an important result, we also want to stress that
we should be cautious before extending these results to any set of models, particularly
given that we only use the information from T<sub>2</sub>. Much work is still needed to devise new
tests and estimation procedures for a wider set of demographic models and using more
genomic information, including recombination patterns as in the PSMC method (LI and
DURBIN 2011). Beyond the general approach outlined here we would like to mention the
<figure><img src="image_23.png" /><figcaption>Figure 23: Estimations of M (real value M = 20.0)</figcaption></figure>
study of PETER et al. (2010) who also managed to separate one structure and one PSC (Population Size Change) model. These authors used an ABC approach to separate a model of exponential PSC from a model of population structure similar to the StSI model. Their structured model differs from ours by the fact that it is not an equilibrium model. They assumed that the population was behaving like an n-island model in the recent past, until T generations in the past, but that before that time, the ancestral population from which all the demes derived was not structured. When T is very large their model is identical to the StSI, but otherwise it may be quite different. The fact that they managed to separate the two models using an ABC approach is promising as it suggests that there is indeed information in the genetic data for models beyond those that we studied here. We can therefore expect that our approach may be applied to a wider set of models. We should also add that in their study these authors use a much larger sample size (25 diploid individuals corresponding to 50 genes). They used a maximum of 200 microsatellites which corresponds therefore to 10,000 genotypes, a number very close to the maximum number we used here. Altogether our study provides new results and opens up new avenues of research for the distribution of coalescent times under complex models.

## Sampling and population expansions

Recent years have also seen an increasing recognition of the fact that the sampling scheme together with population structure may significantly influence demographic inference (WAKELEY 1999; STÄDLER et al. 2009; CHIKHI et al. 2010; HELLER et al. 2013). For instance, WAKELEY (1999) showed that in the n-island model genes sampled in different demes will exhibit a genealogical tree similar to that expected under a stationary Wright-Fisher model. Since our work was focused on $T_2$ we mostly presented our results under the assumption that the two genes of interest were sampled in the same deme. For diploids this is of course a most reasonable assumption. However, we should note that the results presented above also allow
<figure><img src="image_24.png" /><figcaption>Figure 24: Estimations of M (real value M = 50.0)</figcaption></figure>
us to express the distribution of $T_2$ when the genes are sampled in different demes. We did not explore this issue further here, but it would be important to study the results under such conditions. Interestingly we find that if we assume that the two genes are sampled in two distinct demes, we can detect population expansions rather than bottlenecks. This could happen if we considered a diploid individual whose parents came from different demes. In that case, considering the two genes sampled in the deme where the individual was sampled would be similar to sampling his two parental genes in two different demes. Interestingly, this has also been described by PETER et al. (2010) who found that when the 25 individuals were sampled in different demes, they would detect population size expansions rather than bottlenecks. Our results are therefore in agreement with theirs. Similarly HELLER et al. (2013) also found that some signals of population expansion could be detected under scattered sampling schemes.

## Conclusion: islands within individuals

To conclude, our results provide a general framework that can in theory be applied to whole families of models. We showed for the first time that genomic data from a single individual could be used to estimate parameters that have to our knowledge never been estimated. In particular we showed that we were able to estimate the number of islands (and the number of migrants) in the StSI model. This means that one can in principle use genomic data from non model organisms to determine how many islands make up the metapopulation from which one single individual was sampled. This is of course true as well for model organisms but it is particularly meaningful for species for which the number of individuals with genomic data is limited. Our ability to estimate *n* is one of the most powerful results of our study. While such estimates should not be taken at face value, they surely should be obtained across species for comparative analyses. Also, during the last decade there has been a major effort to use programs such as STRUCTURE (PRITCHARD et al. 2000) to estimate the number of
Figures 25, 26, 27, 28, 29 and 30 show the results of the estimation of *n* for all the values of *M* and all the sample sizes, under the StSI model. The values where $M = (0.1, 0.2, 0.5, 1, 5, 10, 20, 50)$, $n = (2, 4, 10, 20, 50, 100)$ and $n_L = (20, 50, 100, 200, 500, 1000, 10000)$

<figure><img src="image_25.png" /><figcaption>Figure 25: Estimations of *n* (real value *n* = 2.0)</figcaption></figure>
"subpopulations" within a particular sample. Our work suggests that we might in principle provide additional results with only one individual. It is important to stress though that the answer provided here is very different from those obtained with STRUCTURE and similar methods and programs (PRITCHARD et al. 2000; GUILLOT et al. 2005; CHEN et al. 2007; CORANDER et al. 2004). We do not aim at identifying the populations from which a set of individuals come. Rather we show that his/her genome informs us on the whole set of populations. In other words, even though we assume that there are *n* populations linked by gene flow, we show that each individual, is somehow a genomic patchwork from this (poorly) sampled metapopulation. We find these results reassuring, in an era where genomic data are used to confine individuals to one population and where division rather than connectivity is stressed.
<figure><img src="image_26.png" /><figcaption>Figure 26: Estimations of *n* (real value *n* = 4.0)</figcaption></figure>
ACKNOWLEDGEMENTS

527 We are grateful to S. Boitard and S. Grusea for numerous and fruitful discussions and
528 input throughout the development of this work. We thank I. Paz, I. Pais, J. Salmona,
529 J. Chave for discussion and helpful comments that improved the clarity of the text. We
530 thank J. Howard and the Fundação Calouste Gulbenkian for their continuous support. We
531 are grateful to the genotoul bioinformatics platform Toulouse Midi-Pyrenees for providing
532 computing and storage resources. This work was partly performed using HPC resources from
533 CALMIP (Grant 2012 - projects 43 and 44) from Toulouse, France. This study was funded
534 by the Fundação para a Ciência e Tecnologia (ref. PTDC/BIA- BIC/4476/2012), the Projets
535 Exploratoires Pluridisciplinaires (PEPS 2012 Bio-Maths-Info) project, the LABEX entitled
536 TULIP (ANR-10-LABX-41) as well as the Pôle de Recherche et d'Enseignement Supérieur
537 (PRES) and the Région Midi-Pyrénées, France.
<figure><img src="image_27.png" /><figcaption>Figure 27: Estimations of *n* (real value *n* = 10.0)</figcaption></figure>
538

539

BEAUMONT, M. A., 1999 Detecting population expansion and decline using microsatellites. Genetics 153(4): 2013-2029.

540

541

BEAUMONT, M. A., 2010 Approximate Bayesian computation in evolution and ecology. Annual review of ecology, evolution, and systematics 41: 379-406.

542

543

BEAUMONT, M. A., W. ZHANG, and D. J. BALDING, 2002 Approximate Bayesian computation in population genetics. Genetics 162(4): 2025-2035.

544

545

BROQUET, T., S. ANGELONE, J. JAQUIERY, P. JOLY, J.-P. LENA, T. LENGAGNE, S. PLENET, E. LUQUET, and N. PERRIN, 2010 Genetic bottlenecks driven by population disconnection. Conservation Biology 24(6): 1596-1605.

546

547

548

CHEN, C., E. DURAND, F. FORBES, and O. FRANÇOIS, 2007 Bayesian clustering algorithms ascertaining spatial population structure: a new computer program and a comparison study. Molecular Ecology Notes 7(5): 747-756.

549

550

551

CHIKHI, L., V. C. SOUSA, P. LUISI, B. GOOSSENS, and M. A. BEAUMONT, 2010 The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes. Genetics 186(3): 983-995.

552

553

554

CORANDER, J., P. WALDMANN, P. MARTTINEN, and M. J. SILLANPÄÄ, 2004 BAPS 2: enhanced possibilities for the analysis of genetic population structure. Bioinformatics 20(15): 2363-2369.

555

556

557

CORNUET, J.-M., F. SANTOS, M. A. BEAUMONT, C. P. ROBERT, J.-M. MARIN, D. J. BALDING, T. GUILLEMAUD, and A. ESTOUP, 2008 Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation. Bioinformatics 24(23): 2713-2719.

558

559

560

561

CURRAT, M., L. EXCOFFIER, W. MADDISON, S. P. OTTO, N. RAY, M. C. WHITLOCK, and S. YEAMAN, 2006 Comment on "Ongoing adaptive evolution of ASPM, a brain

562
<figure><img src="image_28.png" /><figcaption>Figure 28: Estimations of *n* (real value *n* = 20.0)</figcaption></figure>
563

size determinant in *Homo sapiens"* and "Microcephalin, a gene regulating brain size,
continues to evolve adaptively in humans". *Science* 313(5784): 172a-172a.

564

565

DONNELLY, P. and S. TAVARE, 1995 Coalescents and genealogical structure under neutrality. *Annual review of genetics* 29(1): 401-421.

566

567

DURRETT, R., 2008 *Probability models for DNA sequence evolution*. Springer.

568

GIROD, C., R. VITALIS, R. LEBLOIS, and H. FRÉVILLE, 2011 Inferring Population Decline and Expansion From Microsatellite Data: A Simulation-Based Evaluation of the Msvar Method. *Genetics* 188(1): 165-179.

569

570

571

GOOSSENS, B., L. CHIKHI, M. ANCRENAZ, I. LACKMAN-ANCRENAZ, P. ANDAU, and M. W. BRUFORD, 2006 Genetic signature of anthropogenic population collapse in orang-utans. *PLoS Biology* 4(2): e25.

572

573

574

GUILLOT, G., F. MORTIER, and A. ESTOUP, 2005 GENELAND: a computer package for landscape genetics. *Molecular Ecology Notes* 5(3): 712-715.

575

576

HELLER, R., L. CHIKHI, and H. R. SIEGISMUND, 2013 The confounding effect of population structure on Bayesian skyline plot inferences of demographic history. *PLoS One* 8(5): e62992.

577

578

579

HERBOTS, H. M. J. D., 1994 Stochastic models in population genetics: genealogy and genetic differentiation in structured populations. *Ph. D. thesis.*

580

581

JONES, E., T. OLIPHANT, P. PETERSON, and OTHERS, 2001 SciPy: Open source scientific tools for Python. [Online; accessed 2014-11-18].

582

583

LEBLOIS, R., A. ESTOUP, and R. STREIFF, 2006 Genetics of recent habitat contraction and reduction in population size: does isolation by distance matter? *Molecular Ecology* 15(12): 3601-3615.

584

585

586

LI, H. and R. DURBIN, 2011 Inference of human population history from individual whole-genome sequences. *Nature* 475(7357): 493-496.

587
<figure><img src="image_29.png" /><figcaption>Figure 29: Estimations of *n* (real value *n* = 50.0)</figcaption></figure>
588

NELDER, J. A. and R. MEAD, 1965 A simplex method for function minimization. The computer journal 7(4): 308–313.

589

590

NIELSEN, R. and M. A. BEAUMONT, 2009 Statistical inferences in phylogeography. Molecular Ecology 18(6): 1034–1047.

591

592

PAZ-VINAS, I., E. QUÉMÉRÉ, L. CHIKHI, G. LOOT, and S. BLANCHET, 2013 The demographic history of populations experiencing asymmetric gene flow: combining simulated and empirical data. Molecular ecology 22(12): 3279–3291.

593

594

595

PETER, B. M., D. WEGMANN, and L. EXCOFFIER, 2010 Distinguishing between population bottleneck and population subdivision by a Bayesian model choice procedure. Molecular Ecology 19(21): 4648–4660.

596

597

598

PRITCHARD, J. K., M. STEPHENS, and P. DONNELLY, 2000 Inference of population structure using multilocus genotype data. Genetics 155(2): 945–959.

599

600

QUÉMÉRÉ, E., X. AMELOT, J. PIERSON, B. CROUAU-ROY, and L. CHIKHI, 2012 Genetic data suggest a natural prehuman origin of open habitats in northern Madagascar and question the deforestation narrative in this region. Proceedings of the National Academy of Sciences 109(32): 13028–13033.

601

602

603

604

ROGERS, A. R. and H. HARPENDING, 1992 Population growth makes waves in the distribution of pairwise genetic differences. Molecular biology and evolution 9(3): 552–569.

605

606

SALMONA, J., M. SALAMOLARD, D. FOUILLOT, T. GHESTEMME, J. LAROSE, J.-F. CENTON, V. SOUSA, D. A. DAWSON, C. THEBAUD, and L. CHIKHI, 2012 Signature of a pre-human population decline in the critically endangered Reunion Island endemic forest bird Coracina newtoni. PloS one 7(8): e43524.

607

608

609

610

SOUSA, V. C., M. A. BEAUMONT, P. FERNANDES, M. M. COELHO, and L. CHIKHI, 2012, May)Population divergence with or without admixture: selecting models using an ABC approach. Heredity 108(5): 521–530.

611

612
<figure><img src="image_30.png" /><figcaption>Figure 30: Estimations of n (real value n = 100.0)</figcaption></figure>
613

614

615

STÄDLER, T., B. HAUBOLD, C. MERINO, W. STEPHAN, and P. PFAFFELHUBER, 2009 The impact of sampling schemes on the site frequency spectrum in nonequilibrium subdivided populations. Genetics 182(1): 205-216.

616

617

TAVARÉ, S., 2004 Part I: Ancestral inference in population genetics. In *Lectures on probability theory and statistics*, pp. 1-188. Springer.

618

619

VITTI, J. J., S. R. GROSSMAN, and P. C. SABETI, 2013 Detecting natural selection in genomic data. Annual review of genetics 47: 97-120.

620

621

WAKELEY, J., 1999 Nonequilibrium migration in human history. Genetics 153(4): 1863-1871.

622

WRIGHT, S., 1931 Evolution in Mendelian populations. Genetics 16(2): 97.
<figure><img src="image_2.png" /><figcaption>Figure 1: Demographic models. (a): Single step population size change (SSPSC) model. The x-axis represents *t*, the time to the past in units of generations scaled by the number of genes. At time *t* = *T*, the population size changes instantaneously from *N*₁ to *N*₀ by a factor *α*. The y-axis represents the population sizes in units of *N*₀ (i.e. *N*(*t*)/*N*(*0*)). (b): Structured symmetrical island (StSI) model for *n* = 5 islands. Each circle represents a deme of size *N*. All demes are connected to each other by symmetrical gene flow, represented by the edges. In this example the total number of genes is 5*N*.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Expected value and Variance of $T_2$ under the SSPSC and StSI models. This figure illustrates how both models can have the same pair of values ($E(T_2)$, $Var(T_2)$) for many sets of parameters. For the SSPSC model the time at which the population size change occured was fixed to $T = 0.1$ whereas $\alpha$ varied from 1 to 100 in one case, and $T = 1.09$, whereas $\alpha$ varied from 1 to 200 in the other case. For the StSI model the migration rate was fixed to $M = 9$ and $M = 0.5$, whereas $n$ varies from 2 to 100.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Density of $T_2$ under the SSPSC and StSI models. Two sets of parameter values (panels (a) and (b), respectively) were chosen on the basis that expectations and variances were close. Panel (a): Density for the SSPSC model with $T = 0.1$ and $\alpha = 10.94$, and for the StSI model with $M = 9$ and $n = 10$. For this set of parameters we have $E(T_2^{SSPSC}) = 9.994$, and $E(T_2^{StSI}) = 10$, $Var(T_2^{SSPSC}) = 118.7$ and $Var(T_2^{StSI}) = 118.0$. Panel (b): The same, but for $T = 1.09$ and $\alpha = 125.91$, and for $M = 0.5$ and $n = 43$. The corresponding expectations and variances are $E(T_2^{SSPSC}) = 42.997$, and $E(T_2^{StSI}) = 43$, $Var(T_2^{SSPSC}) = 8905$ and $Var(T_2^{StSI}) = 8905$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Estimation of $\alpha$ and $n$. Panels (a), (c) and (e): Estimation of $\alpha$ under the SSPSC model for different sample sizes and $T$ values. Simulations performed with $\alpha = 10$ and $T = (0.1, 1, 2)$. Panels (b), (d) and (f): Estimation of $n$ under the StSI model for different sample sizes and $M$ values. Simulations performed with $n = 10$ and $M = (10, 1, 0.5)$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Estimation of $T$ and $M$. Panels (a), (c), (e)): Estimation of $T$ under the SSPSC model for different sample sizes and values of $\alpha$. Simulations performed with $\alpha = (2, 20, 100)$ and $T = 0.2$. Panels (b), (d), (f): Estimation of $M$ under the StSI model for different sample sizes and values of $n$. Simulations performed with $n = (2, 20, 100)$ and $M = 5$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Proportion of rejected data sets simulated under the SSPSC model. Panels (a), (c) and (e): the reference model is the StSI model. Panels (b), (d), and (f): the reference model is the SSPSC, i.e. the model under which the data were simulated. Note that for the abscissa we used $2n_L$ instead of $n_L$ because in order to perform the KS test it is necessary to first estimate the parameters using $n_L$ loci and then an independent set of $n_L$ values of $T_2$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Proportion of rejected data sets simulated under the StSI model. Panels (a), (c), and (e): the reference model is the SSPSC. Panels (b), (d), and (f): the reference model is the StSI model, i.e. the model under which the data were simulated. Note that for the abscissa we used $2n_L$ instead of $n_L$ because in order to perform the KS test it is necessary to first estimate the parameters using $n_L$ loci and then an independent set of $n_L$ values of $T_2$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Relationships between parameters of the models</figcaption></figure>