Linkage disequlibrium under recurrent bottlenecks

Genetic variation at single neutral loci has been investigated in great detail for population models under different demographic processes, such as population expansions, single bottlenecks, or genetic hitchhiking caused by nearby selective sweeps (see for example (Eriksson et al., 2008) for a review of such models). All of of these models, either with non-overlapping generations, such as the Wright-Fisher model (Fisher, 1930/1999; Wright, 1931), or with continuous reproduction and mortality, such as the Moran model Moran (1958), have in common the underlying assumption of random mating.

Real biological populations exhibit abundance fluctuations on both short and long time scales, caused by e. g. environmental and ecological changes. Such size fluctuations in the form of repeated bottlenecks are characteristic of populations expanding into new territories: examples include the human out-of-Africa scenario (Liu et al., 2006; Ramachandran et al., 2005), the accompanying expansion of the parasite Plasmodicum falciparum causing severe malaria (Tanabe et al., 2010), and the recolonization by the marine snail Littorina saxatilis of Sweden’s west coast archipelago (Johannesson, 2003). It is common practice to accommodate such fluctuations in the theory by using an effective population size instead of the census population size (see (Ewens, 1982) for a review of different measures of the effective population size).

Recent research has highlighted the importance of two competing time scales in the context of effective population-size approximations: the coalescent time scale (which is inversely proportional to the population size and reflects the time to the most recent common ancestor (MRCA) in populations with constant population size), and the time scale of the population-size fluctuations. Clearly, relatively slow demographic fluctuations can be ignored and the effective population size can be approximated by the initial population size Sjödin et al. (2005). In the opposite case of very frequent demographic fluctuations, it has been argued (Wright, 1938; Crow and Kimura, 1970) that genetic variation of a population with varying population size is well described in terms of a population with constant population size $N_{\rm eff}$, given by the harmonic average of the population size $N_{\tau}$ in generation $\tau$:

$$N_{\rm eff}=\lim_{T\rightarrow\infty}\Bigl{(}\frac{1}{T}\sum_{\tau=0}^{T}\frac{1}{N_{\tau}}\Bigr{)}^{-1}\,\,.$$

(1)

See (Sjödin et al., 2005; Jagers and Sagitov, 2004; Wakeley and Sargsyan, 2009) for recent developments of this concept.

Thus, for both fast and slow demographic fluctuations, the statistical properties of gene genealogies agree with those of the constant population-size model. By contrast, when both time scales are of the same order, it has been shown (Nordborg and Krone, 2003; Kaj and Krone, 2003; Sjödin et al., 2005; Eriksson et al., 2010) that the distribution of total branch lengths in samples of one-locus gene genealogies does not in general agree with that predicted by the standard coalescent approximation. Especially, when subject to periodic fluctuations, the total branch length of gene genealogies may exhibit a maximum for periods matching the coalescent time scale (Eriksson et al., 2010). In summary, the effect of population-size fluctuations upon genetic variation of a single locus is well understood.

But how do population-size fluctuations affect multi-locus patterns of genetic variation on the same chromosome? Such patterns are influenced by recombination. Genetic recombination introduces a new time scale which is inversely proportional to the probability of recombination $r$ between a pair of loci in a single generation.

Genetic recombination plays an important role in shaping empirically observed multi-locus patterns of genetic variation in biological populations. Measures of linkage disequilibrium (LD) quantify the degree of association of genetic variation at pairs of loci on the same chromosome. Common measures of LD, such as $\hat{r}^{2}$ Hill and Robertson (1968), and its approximation $\sigma^{2}_{d}$ Ohta and Kimura (1971); McVean (2002), depend upon the frequencies of alleles at two loci. These measures are thus closely related to the covariance of the times (i. e. the number of generations) to the MRCA of the underlying gene genealogies (McVean, 2002).

In order to illustrate this concept, we show in Fig. 1 our results of the Wright-Fisher dynamics (grey lines) of a population experiencing recurrent bottlenecks, with random durations and separations between bottlenecks. Panels a and b show results for two different sets of parameters. Details are given in Section II. The single-locus properties of a similar model were studied by Sjödin et al. (2005) and also by Eriksson et al. (2010). Each grey line in Fig. 1 shows the covariance of the times to the MRCA for two chromosomes at two loci, as a function of genetic distance along the chromosomes, corresponding to a single realisation of the sequence of bottlenecks, $\mathcal{D}$ (averaged over all pairs of loci the same distance apart). The red lines are the averages of the covariances within each panel. In panel a, the bottlenecks happen frequently and have a short duration; in this case, the single-locus properties are expected to be in good agreement with those of a population with the constant effective population size, given by the harmonic time average, Eq. (1), of the population size $N_{\tau}$. For such populations, the coalescent approximation predicts (Griffiths, 1981; Hudson, 1983, 1990)

$$\text{cov}[t_{a(ij)},t_{b(ij)}]=x_{\rm eff}^{2}\frac{Rx_{\rm eff}+18}{(Rx_{\rm eff})^{2}+13Rx_{\rm eff}+18},$$

(2)

where $t_{a(ij)}$ and $t_{b(ij)}$ are the times to the MRCA of two loci (called $a$ and $b$), in a sample of two chromosomes (denoted by $i$ and $j$). Moreover, $R=2N_{0}r$ is the scaled recombination rate, and $x_{\rm eff}=N_{\rm eff}/N_{0}$ is the effective population size relative to $N_{0}$, the population size at the present time. Units of time are chosen so that $\tau=\lfloor tN_{0}\rfloor$, and $\lfloor tN_{0}\rfloor$ is the largest integer not larger than $tN_{0}$. Effective population-size approximations, according to Eq. (2), are shown as dashed lines in Fig. 1. In Fig. 1a we observe good agreement between Eq. (2) and the average covariance of the times to the MRCA obtained numerically. However, in Fig. 1b, Eq. (2) agrees with our numerically obtained covariance only for short genetic distances. For large genetic distances, by contrast, the numerical results show that the covariance decreases much more slowly than expected according to Eq. (2). Thus, the results shown in Fig. 1b imply long-range LD in two-locus gene genealogies.

The examples shown in Fig. 1 open many questions in relation to multi-locus gene-genealogies. What are the conditions for the effective population-size approximation to be valid in the multi-locus case? Why does it fail when these conditions are not met? How significant are deviations of the exact result from the effective population-size approximation? Why does long-range linkage appear in some cases? How large are fluctuations around the covariance of the coalescent times, averaged over an ensemble of gene-genealogies and over different demographic histories? What is the significance of the fluctuations around such averages for data analysis?

The aim of this paper is to provide answers to the above questions, by computing the covariance of the times to the MRCA under a model of recurrent bottlenecks introduced in Section II. Our analysis enables us to qualitatively and quantitatively determine the effects of fluctuating population size on the two-locus statistics in terms of the time scales of population-size fluctuations, of coalescence, and of recombination. Using both a numerical and an analytical approach, we estimate the range of validity of the coalescent effective population-size approximation for the two-locus case. We find that the coalescent effective population-size approximation inevitably fails for large recombination rates: the failure is sometimes minor (as in the case shown in Fig. 1a) and sometimes significant (as in the case shown in Fig. 1b). By taking different limits of the parameters of the model, we provide both qualitative and quantitative understanding of how the constant effective population-size approximation may fail. Finally, when bottlenecks are severe, we show that gene genealogies can be approximated by those of a constant-sized population with simultaneous multiple pairwise coalescent events (the so-called Xi-coalescent described by Schweinsberg (2000) and Möhle and Sagitov (2001)). Last but not least, we demonstrate that $\sigma^{2}_{d}$ is surprisingly little affected by the long-range covariances.

We use a Wright-Fisher model of a population of chromosomes, to trace the ancestry of $L$ loci on a pair of chromosomes backwards in time, until the most recent common ancestor of each locus is found. In each generation $\tau$, a set of parents is chosen randomly from the previous generation. Genetic recombination of a pair of chromosomes occurs independently with probability $r$ between each adjacent pair of loci.

To investigate the role of population-size fluctuations, we consider a model of recurrent bottlenecks in which the population size can take one of two values, $N_{0}$ or $xN_{0}$. Here $xN_{0}$ is the population size in the bottlenecks, and $0<x<1$. The probability of changing the population size, going one generation back in time, depends on the current population size. The switching probability is $p$ and $q$ when the current population size is $N_{0}$ and $xN_{0}$, respectively. Hence, the expected durations of the high and the low population-size phases are $1/p$ and $1/q$ generations, respectively. The population size in the first generation is taken to be $N_{0}$. We note that for the one-locus case, such a model has been investigated by Sjödin et al. (2005) and also by Eriksson et al. (2010).

Fig. 2 illustrates population-size fluctuations in this recurrent bottleneck model (panels a and b) and examples of gene genealogies of two loci (called $a$ and $b$) in a sample of two chromosomes (panels c and d). In panels c and d, generations with low population size ($xN_{0}$) are marked with yellow, otherwise the population size is $N_{0}$. Each chromosome is represented by a pair of lines (red and blue lines correspond to loci $a$ and $b$, respectively). In the generations where a common ancestor is found for a pair of ancestral lines, or a recombination between two loci occurs, the chromosomes are represented by circles instead of lines. The MRCA of a locus is shown as a filled circle. In some cases recombination causes the ancestry of one locus to become associated with a chromosome that lacks direct descendants in the sample (grey circles). The ancestries of such segments of DNA are irrelevant to the gene genealogy of the sample, and these ancestral lines are not traced further.

As Fig. 1 shows, the covariance of the times to the MRCA at two loci (called $a$ and $b$) depends on the particular sequence of bottlenecks in the history of the population. Let $\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]$ denote the covariance conditional on the demographic history $\mathcal{D}$. In the following, we average the conditional covariance over random demographic histories to obtain (see red lines in Fig. 1):

	$\displaystyle\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle$	$\displaystyle=\langle\langle t_{a(ij)}t_{b(ij)}|\mathcal{D}\rangle-\langle t_{a(ij)}|\mathcal{D}\rangle\langle t_{b(ij)}|\mathcal{D}\rangle\rangle$
		$\displaystyle=\langle t_{a(ij)}t_{b(ij)}\rangle-\langle\langle t_{a(ij)}|\mathcal{D}\rangle^{2}\rangle\,\,,$		(3)

where $\langle\dots|\mathcal{D}\rangle$ denotes the expectation conditional on the particular demographic history $\mathcal{D}$. In the second equality we have used that the expected times to the MRCA are the same for both loci. Note that the averaged conditional covariance is not the same as the unconditional covariance of the times to the MRCA for the full process (given by $\langle t_{a(ij)}t_{b(ij)}\rangle-\langle t_{a(ij)}\rangle^{2}$).

We now derive an approximate expression for $\langle{\text{cov}}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle$ using the coalescent approximation (Kingman, 1982), valid in the limit of large population sizes. Thus, in the following calculations, we assume $N_{0}\gg 1$, and $xN_{0}\gg 1$. As is usual in coalescent calculations, it is convenient to scale the rates $p$ and $q$ in the Wright-Fisher model by a suitable representative of the population size. Defining the rates $\lambda=pN_{0}$, and $\lambda_{\text{B}}=qxN_{0}$ allows us to express the probability for two lines to coalesce between two consecutive bottlenecks as $(1+\lambda)^{-1}$. Correspondingly, the probability for two lines to coalesce during a single bottleneck is given by $(1+\lambda_{\rm B})^{-1}$. The units of time are chosen so that $\tau=\lfloor tN_{0}\rfloor$, we take $R=2N_{0}r$ as the scaled recombination rate as mentioned in the introduction. In these units we denote the population size at time $t$ by $N(t)$, indicating that in the limit of $N_{0}\rightarrow\infty$, $t$ becomes a continuous variable.

We find the following expression for the term $\langle\langle t_{a(ij)}|\mathcal{D}\rangle\langle t_{b(ij)}|\mathcal{D}\rangle\rangle\equiv\langle\langle t_{a(ij)}|\mathcal{D}\rangle^{2}\rangle$, occurring in Eq. (III):

$$\langle\langle t_{a(ij)}|\mathcal{D}\rangle^{2}\rangle=\frac{\lambda_{\rm B}(2x\lambda+\lambda_{\rm B}+3)+x\lambda(x\lambda+x+2)+2}{(\lambda+\lambda_{\rm B}+1)(\lambda+\lambda_{\rm B}+2)}\,\,.$$

(4)

Details of the derivation of this result are summarised in Appendix A.

In order to evaluate the remaining term in Eq. (III), $\langle t_{a(ij)}t_{b(ij)}\rangle$, we adapt the method described by Eriksson and Mehlig (2004) for calculating the covariance of the times to the MRCA at two loci, to the present model of recurrent bottlenecks. As can be seen in Fig. 2c, d, there are only a small number of possible combinations of ancestral lines in gene genealogies of two loci for two chromosomes. Thus, we can write down a Markov process for how states of the ancestral lines change along the gene genealogy. The corresponding graph is shown in Fig. 3, where the vertices represent states (combinations of ancestral lines), and the edges represent transitions between the states (the transition rates $w_{ji}$ from $i$ to $j$ are shown along the edges from $i$ to $j$). The vertices labeled by a prime correspond to states in a bottleneck.

Following Eriksson and Mehlig (2004) we note that the expectation $\langle t_{a(ij)}t_{b(ij)}\rangle$ is determined by a six-dimensional sub-graph of the graph shown in Fig. 3, consisting of the vertices $1,2,3,1^{\prime},2^{\prime}$, and $3^{\prime}$. Let $\bf{M}$ be the corresponding $6\times 6$ transition matrix. Its off-diagonal elements are given by the transition rates, $w_{ji}$, from state $i$ to state $j$. The diagonal elements $M_{ii}$ are equal to the negative sum of the rates of all edges leaving node $i$ in the graph, i. e. $M_{ii}=-\sum_{j\neq i}M_{ji}$. At the present time ($t=0$), the system is in state $1$ in the graph. This is represented by the vector

$\displaystyle\mbox{{\boldsymbol{$v$}}}=\begin{bmatrix}1&0&0&0&0&0\end{bmatrix}^{\sf T},$

(5)

where $\sf T$ denotes the transpose. It is convenient to combine the transition rates $w_{ji}$ into vectors and matrices as follows:

	$u$	$\displaystyle=\begin{bmatrix}1&0&0&x^{-1}&0&0\end{bmatrix},$
	$\displaystyle\bf{Q}$	$\displaystyle=\begin{bmatrix}0&2&2&0&0&0\\ 0&0&0&0&2x^{-1}&2x^{-1}\end{bmatrix},$
	$\displaystyle\bf{K}$	$\displaystyle=\begin{bmatrix}-(\lambda+1)&\lambda_{\text{B}}x^{-1}\\ \lambda&-(\lambda_{\text{B}}+1)x^{-1}\end{bmatrix},$		(6)
	$c$	$\displaystyle=\begin{bmatrix}1&x^{-1}\end{bmatrix}\,.$

In terms of these vectors and matrices we can write Eriksson and Mehlig (2004)

$\displaystyle\langle t_{a(ij)}t_{b(ij)}\rangle=\int_{0}^{\infty}dt_{1}\,t_{1}^{2}\,{\mbox{{\boldsymbol{$u$}}}}\,{\rm e}^{{\bf M}\,t_{1}}{\mbox{{\boldsymbol{$v$}}}}+\int_{0}^{\infty}dt_{1}\int_{t_{1}}^{\infty}dt_{2}\,t_{1}t_{2}\,{\mbox{{\boldsymbol{$c$}}}}\,{\rm e}^{{\bf K}\left(t_{2}-t_{1}\right)}\,{\bf Q}\,{\rm e}^{{\bf M}\,t_{1}}{\mbox{{\boldsymbol{$v$}}}},$

(7)

where $t_{1}$ and $t_{2}$ are the first and second coalescent events, respectively (i. e. $\min(t_{a(ij)},t_{b(ij)})$ and $\max(t_{a(ij)},t_{b(ij)})$). The first term corresponds to a common MRCA of loci $a$ and $b$ (i. e. a transition from states $1$ or $1^{\prime}$ to state $5$), and the second to different MRCA (transitions from states $2$ or $3$ to state $4$, or from states $2^{\prime}$ or $3^{\prime}$ to state $4^{\prime}$, followed by a transition to state $5$). The matrix $\bf{K}$ describes the change in the population size due to recurrent bottlenecks, whereas the vector $c$ contains the coalescent rates in each population-size regime. Both $\bf{M}$ and $\bf{K}$ have negative real eigenvalues. Hence, the integrals in Eq. (7) can be evaluated in terms of matrix inverses Eriksson and Mehlig (2004):

$\displaystyle\langle t_{a(ij)}t_{b(ij)}\rangle=2\,{\mbox{{\boldsymbol{$u$}}}}\,(-{\bf M})^{-3}\,{\mbox{{\boldsymbol{$v$}}}}+2\,{\mbox{{\boldsymbol{$c$}}}}\,(-{\bf K})^{-2}\left\{(-{\bf K})^{-1}\,{\bf Q}+{\bf Q}\,(-{\bf M})^{-1}\right\}(-{\bf M})^{-2}\,{\mbox{{\boldsymbol{$v$}}}}.$

(8)

Combining Eqs. (4) and (8) yields

$\displaystyle\langle\mbox{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle=\frac{R^{3}C_{3}+R^{2}C_{2}+RC_{1}+C_{0}}{R^{4}D_{4}+R^{3}D_{3}+R^{2}D_{2}+RD_{1}+D_{0}}\,\,,$

(9)

where the coefficients $C_{i}$ and $D_{i}$ are functions of parameters $x$, $\lambda$ and $\lambda_{\rm B}$ given in Appendix C.

We now discuss three special limits of this result. First, when the time to the first bottleneck is much longer than the expected times to the MRCA for two chromosomes at two loci in the large population-size regime (i. e. $\lambda\ll 1$), the subsequent bottlenecks are irrelevant to the mean covariance. In this case, Eq. (9) reduces to Eq. (2), the expression valid in the case of constant population size with $x_{\text{eff}}=1$.

Second, assume that the population-size fluctuations are rapid compared to all other processes (this case is described by taking the limit $\lambda\rightarrow\infty$ and $\lambda_{\rm B}\rightarrow\infty$ in such a way that the ratio $\lambda/\lambda_{\rm B}$ is kept constant). Keeping only the leading order of $\lambda$ and $\lambda_{\text{B}}$ in the numerator and the denominator of (9) yields Eq. (2), with $x_{\text{eff}}=(x\lambda+\lambda_{\text{B}})/(\lambda+\lambda_{\text{B}})$. This is the harmonic time average (1) of $N_{\tau}/N_{0}$ in the recurrent bottleneck model.

Third, we consider the case of severe bottlenecks ($x\ll 1$). Because the expected duration of a bottleneck is given by $x\lambda_{\text{B}}^{-1}$, this regime implies that bottlenecks are typically of short duration. Such demographic histories can occur during range expansions, where small groups of animals repeatedly colonize new areas (examples of this kind are given in the introduction). This case can be treated analytically by taking the leading order of $x^{-1}$ in the numerator and denominator of Eq. (9). We find:

$\displaystyle\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle\approx\frac{R^{2}A_{2}+RA_{1}+A_{0}}{R^{2}B_{2}+RB_{1}+B_{0}},$

(10)

where $A_{i}$ and $B_{i}$ are functions of parameters $\lambda$ and $\lambda_{\rm B}$ given in Appendix C. Note that this function reaches a plateau for large values of $R$:

	$\displaystyle\langle\mbox{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle$	$\displaystyle\approx\frac{A_{2}}{B_{2}}$		(11)
	$\displaystyle=$	$\displaystyle\frac{2\lambda_{\rm B}(1+\lambda_{\rm B})\lambda}{(1+\lambda_{\rm B}+\lambda)(2+\lambda_{\rm B}+\lambda)(9(2+\lambda)+\lambda_{\rm B}(27+8\lambda+\lambda_{\rm B}(10+\lambda_{\rm B}+\lambda)))}\,.$

This expression implies long-range linkage disequilibrium since the right hand side does not depend upon $R$.

In order to further illustrate the role of the relevant time scales of genetic drift and recombination in shaping LD, we compare the full coalescent result, Eq. (9), and the different limiting cases considered in Section III, to the average covariance calculated from the Wright-Fisher simulations. The comparisons are shown in Fig. 4.

The parameters used in Fig. 4a correspond to rapid population-size fluctuations (the second example described in the previous section). As can be seen in Fig. 4a, the agreement between the numerical result (red line) and the approximation (2) (dashed line) is good for a wide range of recombination rates. A small disagreement appears at large values of $R$, more precisely at $R\approx 100$. This discrepancy is expected, as at such large recombination rates, the population-size fluctuations are no longer rapid compared to the process of recombination. In summary, in the case of rapid population-size fluctuations, the results of the Wright-Fisher simulations are well approximated by Eq. (2) as long as genetic recombination does not occur too frequently.

Fig. 4b shows results for parameters corresponding to severe bottlenecks (the third example described in the previous section). As we noted already in the introduction, this case exhibits long-range linkage disequilibrium which cannot be accurately described by the effective population-size approximation (2). We observe that the average covariance curve can be divided into three regions where the curve behaves qualitatively differently. First, for very small recombination rates, the covariance can be approximated using the harmonic average effective population size. This is expected since, in this region, the population-size fluctuations are fast compared to the process of recombination. Second, the constant effective population-size approximation breaks down when $R\approx\lambda=10$. For recombination rates larger than this value, by contrast, recombination occurs frequently between bottlenecks but rarely within, and multiple coalescent events during bottlenecks become frequent. Since the bottlenecks are short in this regime, this may lead to long-range LD. In this regime, the covariance is approximated by Eq. (10). We observe that, for a large range of recombination rates, Eq. (10) is in excellent agreement with both the full coalescent result, Eq. (9), and the simulations. Third, we see that the agreement between Eq. (10) and the full coalescent result breaks down when recombination events can no longer be ignored in the bottlenecks (i. e. when $R$ is of the same order as the rate of leaving a bottleneck, $\lambda_{\text{B}}x^{-1}$, or larger). For still larger recombination rates, only the full coalescent result agrees with the Wright-Fisher simulations. The slight deviations between the Wright-Fisher simulations and the coalescent result (9), visible in Fig. 4, are discussed in the concluding section.

In the previous sections we have shown how sustained population-size fluctuations in the form of recurrent bottlenecks give rise to long-range linkage disequilibrium, measured by the covariance (III) of the times to the MRCA. An important question is how such population-size fluctuations affect more common measures of linkage disequilibrium such as, for example, $\hat{r}^{2}$, and its approximation, $\sigma^{2}_{d}$ (Ohta and Kimura, 1971; McVean, 2002). In this section we discuss the effect of recurrent bottlenecks upon $\sigma^{2}_{d}$. While the covariance of the times to the MRCA is obtained by comparing two chromosomes, the measure $\sigma_{d}^{2}$ is computed in a large sample of chromosomes. In the following we consider the limit of infinite sample size. In this limit, and in terms of the expectations and covariances conditional on the demographic history $\mathcal{D}$, Eq. (9) in McVean (2002) becomes:

$\displaystyle\sigma^{2}_{d}=\left\langle\frac{\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]-2\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ik)}]+\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(kl)}]}{\langle t_{a(ij)}|\mathcal{D}\rangle^{2}+\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(kl)}]}\right\rangle\,.$

(12)

As before, $a$ and $b$ denote two loci, and $i$, $j$, $k$, and $l$ refer to four different chromosomes in a large sample. The main properties of this measure are determined by how the numerator depends on the recombination rate McVean (2002). In order to simplify the analysis, we therefore focus on the expected value of the numerator, i. e.

$$\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle-2\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ik)}]\rangle+\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(kl)}]\rangle.$$

The covariance $\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle$ is given by Eq. (9). The covariances $\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ik)}]\rangle$ and $\langle\text{cov}_{\mathcal{D}}[t_{a(ij)},t_{b(kl)}]\rangle$ can be calculated in the same way as Eq. (9) was obtained, but starting from different initial conditions (McVean, 2002; Eriksson and Mehlig, 2004). In our Markov representation, this corresponds to taking ${\mbox{{\boldsymbol{$v$}}}}=[0\ 1\ 0\ 0\ 0\ 0]^{T}$ and ${\mbox{{\boldsymbol{$v$}}}}=[0\ 0\ 1\ 0\ 0\ 0]^{T}$, respectively, in Eq. (7).

Fig. 5 shows the relation between the covariances $\langle{\rm cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ij)}]\rangle$ (blue lines), $\langle{\rm cov}_{\mathcal{D}}[t_{a(ij)},t_{b(ik)}]\rangle$ (red lines), and $\langle{\rm cov}_{\mathcal{D}}[t_{a(ij)},t_{b(kl)}]\rangle$ (green lines). As can be seen, in a region of low values of $R$ each covariance is well approximated by the corresponding constant effective population-size approximation, as explained in Section III. Note that in the case shown in panel b, all three covariances exhibit a plateau at the same level, in approximately the same range of $R$ values. Thus, the linear combination of covariances in the numerator of $\sigma^{2}_{d}$ cancel an enhancement relative to the effective population-size approximation. This is the reason why the plateau, present in each covariance, does not show in the numerator of $\sigma^{2}_{d}$. In other words, the information about the sustained population-size fluctuations is not preserved in $\sigma_{d}^{2}$.

The aim of this paper was to provide an understanding on how sustained random population-size fluctuations influence gene-history correlations. Our conclusions are based on both analytical and numerical calculations, which we find to agree well. Using the particular population-size model, depicted in Fig. 2a, we have derived an exact result for the covariance of the times to the MRCA of two loci, Eq. (9). We have discussed three particular limits of our result. First, if the expected times to the MRCA of two loci are less than the expected time to the most recent bottleneck, our model reduces to the constant population-size model with the effective population size equal to the population size at the present time. Second, if the population size fluctuates much faster than the remaining two processes (coalescence and recombination), the coalescent effective population-size approximation again works well, but with the effective population size given by the harmonic average (1). These two cases are consistent with earlier findings of investigations of single-locus properties of populations which exhibit population-size variations (Sjödin et al., 2005; Eriksson et al., 2010). In the third case, bottlenecks are severe (large difference in population sizes in and between bottlenecks) and the time between bottlenecks is assumed to be brief. In this case, the result of the simulations depends on the relation between the time scale of recombination and the time scale of population-size changes. When recombination is the slowest process (i. e. when the recombination rate is sufficiently small), the coalescent effective population-size approximation (with the effective population size given by the harmonic average (1)) is a good approximation (this is essentially the same as in the second case). Conversely, when recombination is the fastest process, the covariance is approximately $1/R$ (same as in the first case). Finally, when recombination is intermediate, slow enough that recombination events during bottlenecks are rare, but fast enough to decorrelate gene genealogies between bottlenecks, the covariance exhibits a plateau. The plateau corresponds to an enhanced covariance of the times to the MRCA, with respect to that expected in the effective constant population-size case. Thus, in this case, pairs of distant loci are expected to exhibit a high degree of linkage.

These conclusions rely on analysing covariances averaged over different demographic histories. This raises the question how typical such averages are. In other words, how large are the fluctuations around the average? Fig. 1 shows that in the case of severe bottlenecks, the fluctuations around the mean covariance are much higher than, for example, in the case of fast population-size fluctuations. This can be explained as follows: in the case of severe bottlenecks, the times to the MRCA of both loci are determined by the time to the bottleneck that hosts two pairwise mergers (for very strong bottlenecks, this is simply the time to the first bottleneck).

The coalescent approximations employed in this paper assume large population sizes. While we generally find very good agreement between the coalescent approximations and the Wright-Fisher dynamics, we observe some deviations, in particular for severe bottlenecks and large recombination rates. The coalescent approximations assume that the time between two recombination events is exponentially distributed. However, this is only accurate in the limit of small recombination rates, $r$. Because the relevant parameter is $R=2N_{0}r$, by increasing $N_{0}$ we reach a better agreement between the corresponding analytical and numerical results in the range of values of $R$ shown in Figs. 1 and 4. We have made additional simulations to confirm that this is the case (not shown).

It is worth mentioning that the result for the case of severe bottlenecks ($x\ll 1$), Eq. (10), can be understood in terms of the so-called Xi-coalescent approximation. Xi-coalescents form a broad family of gene-genealogical models allowing for simultaneous multiple pairwise coalescent events (mergers). The Kingman coalescent is a special case, allowing only for pairwise mergers. See (Schweinsberg, 2000; Möhle and Sagitov, 2001) for detailed descriptions of the family of Xi-coalescents. In the case of severe bottlenecks ($x\ll 1$), coalescent events during a bottleneck may appear as simultaneous multiple mergers (as pointed out also by Birkner et al. (2009)). We show in Appendix C how the Markov process with simultaneous multiple mergers is obtained in this case, and compute the corresponding transition rates $w_{ji}$ . This provides an alternative way of deriving Eq. (10). More importantly, it provides insight into why the plateau forms. It turns out that the plateau arises as a direct consequence of simultaneous multiple mergers. We note that this means that long-range gene-history correlations are also expected in other situations where simultaneous multiple mergers are important. Examples are populations with strongly skewed reproduction laws, or populations subject to selective sweeps.

We conclude with the observation that $\sigma^{2}_{d}$, a measure of LD, fails to show the plateaus present in its constituent covariances (this was observed already in Eriksson and Mehlig (2004) for the case of a single, recent bottleneck). Because of the close link between $\sigma^{2}_{d}$ and $\hat{r}^{2}$, a common measure of LD McVean (2002), this casts doubt on the suitability of such measures for characterising LD (another example is the measure $HR^{2}$ (Sabatti and Risch, 2002)), in populations that may have been subject to recent population bottlenecks and range expansions. A more accurate approach, especially for detecting long-range LD, may be to estimate the covariance of the times to the MRCA directly. For example, simulations show that the covariance of the number of mutations in small windows (e. g. a few hundred nucleotides long) can be used to estimate the covariance of the times to the MRCA Eriksson and Mehlig (2004). However, it remains to investigate which observables are most suitable for detecting long-range dependencies in the underlying gene genealogies for more general demographic histories.

Acknowledgements. Support by Swedish Research Council grants, the Göran Gustafsson stiftelse, and by the Centre for Theoretical Biology at the University of Gothenburg are gratefully acknowledged. AE was supported by a Philip Leverhume Award and a Biotechnology and Biological Sciences Research Council grant (BB/H005854/1).