Abstract
In this article we discuss the ancestry of sequences sampled from the coalescent with recombination with constant population size 2N. We have studied a number of variables based on simulations of sample histories, and some analytical results are derived. Consider the leftmost nucleotide in the sequences. We show that the number of nucleotides sharing a most recent common ancestor (MRCA) with the leftmost nucleotide is ≈log(1 + 4N Lr)/4Nr when two sequences are compared, where L denotes sequence length in nucleotides, and r the recombination rate between any two neighboring nucleotides per generation. For larger samples, the number of nucleotides sharing MRCA with the leftmost nucleotide decreases and becomes almost independent of 4N Lr. Further, we show that a segment of the sequences sharing a MRCA consists in mean of 3/8Nr nucleotides, when two sequences are compared, and that this decreases toward 1/4Nr nucleotides when the whole population is sampled. A measure of the correlation between the genealogies of two nucleotides on two sequences is introduced. We show analytically that even when the nucleotides are separated by a large genetic distance, but share MRCA, the genealogies will show only little correlation. This is surprising, because the time until the two nucleotides shared MRCA is reciprocal to the genetic distance. Using simulations, the mean time until all positions in the sample have found a MRCA increases logarithmically with increasing sequence length and is considerably lower than a theoretically predicted upper bound. On the basis of simulations, it turns out that important properties of the coalescent with recombinations of the whole population are reflected in the properties of a sample of low size.
UNDERSTANDING the genealogical relationship between sequences in a diploid population has been central to recent analyses of the dynamics of sequence evolution at the population level. The stochastic process generating the genealogical relationship between k sampled sequences from a population with constant size N and no recombination was first described by Watterson (1975) and further developed into the theory of the coalescent by Kingman (1982). The process of evolution of sequences subject to both coalescence and recombination in a population was first described by Hudson (1983). In Hudson’s approach the combined coalescent and recombination process is followed back in time until any position in the extant sequences has found a most recent common ancestor (MRCA). Distant positions will not necessarily share the same history, and the ancestral positions can be located on different sequences. However, the genealogies of distinct but linked positions are correlated: Positions far apart have ancestries almost independent of each other, whereas positions close to each other tend to have identical ancestry. Griffiths and Marjoram (1997) proved that the set of MRCAs to a sample is finite for any sample size. But, even positions sharing the same MRCA can have very different histories.
In this article we discuss the ancestry of a sample of k sequences subject to both coalescence and recombination. This is done mainly through simulations of sample histories. The combinatorial complexity of the coalescent with recombination makes exact results difficult to derive and, in most cases, restricted to samples of size 2. We measure the sequence length in expected number of recombinations per sequence per 2N generations, where N is population size. The population size is assumed to be constant from generation to generation.
Our results can be broadly divided into two parts. In the first part, we have focused on the structure of a single MRCA. Consider the MRCA at position 0 of the sequences. Call this ancestor MRCA_{k} (0), where k refers to sample size. If there is no recombination in the history of the sample, all positions q > 0 will share the same MRCA, i.e., MRCA_{k} (0) = MRCA_{k} (q). However, if recombination is present, only a subset of the positions q > 0 will share this MRCA. In the example in Figure 1 positions 0 ≤ q < ¼ and ½ ≤ q ≤ 1 share MRCA spread on two distinct segments, while ¼ ≤ q < ½ share MRCA.
Furthermore, we are interested in the following variables: (1) length of ancestral material that shares MRCA with position 0 [in the example in Figure 1 this amounts to (1  ½) + (¼  0) = ¾], (2) the number of segments into which the positions sharing MRCA with position 0 are partitioned (in the example in Figure 1 this is two segments).
In the second part, we have focused on the time back to MRCAs. The total branch length, G_{k}(q), and the height, T_{k}(q) (time until a MRCA), of the genealogy of a single nucleotide are distributed according to the coalescent without recombination. In contrast, the distribution of the time until all nucleotides have found a MRCA, T_{k} = max{T_{k}(q)q} and the distribution of G_{k} = max{G_{k}(q)q} depend on the total genetic length ρ. We have investigated the expected values of these two variables. Moreover, we discuss a notion of shared sequence ancestry that relates to the correlation between genealogies.
The recombination rate ρ = 4NLr has been varied from 0 to 50 in the simulations. The quantity r is the probability of a recombination event between any two neighboring positions in a sequence per generation, and L is number of nucleotides in a sequence. Let us assume that r is 10^{7} in the human genome and that the effective population size of humans is 10^{4}. With ρ = 50 we have L = ρ/(4Nr) = 50 · 10^{7}/(2 · 10^{4}) = 2.5 · 10^{4} nucleotides. Thus, our simulation results cover the ancestry of a sample of human DNA sequences of length up to 25,000 nucleotides.
THE COALESCENT WITH RECOMBINATION
The model of a population of sequences subject to recombination is the following: Each sequence is L nucleotides long and recombination is assumed to occur to the right of a nucleotide. The population is of constant size N and diploid, i.e., there are 2N sequences in the population.
A new generation is obtained from the present by (1) selecting with probability 1  r a single parent uniformly at random and (2) selecting with probability r two parents uniformly at random and recombining these. Each sequence in the next generation chooses one or two parents in this manner. The collection of these offspring forms the next generation. The process starts at the present and time increases as it goes backward.
This process is transformed into one of a continuous time and continuous sequence by letting N → ∞ and measuring time in 2N generations and by letting L → ∞ and r → 0, such that 4rLN → ρ. Here 2rLN is the expected number of recombinations per 2N sequences per generation. Sequence length is measured in expected number of recombinations per 2N sequences per generation; that is, the entire sequence length is ρ/2. Hudson (1983) showed that the waiting time until a sequence is created by a recombination event from two sequences is exponentially distributed with intensity parameter ρ_{0}/2. For the extant sequences, ρ_{0}/2 is simply the length of the sequences, i.e., ρ_{0} = ρ. For ancestral sequences, ρ_{0}/2 is the length of the interval spanned by regions that have ancestral material. Note that this interval can include regions with nonancestral material (cf. Figure 1). The recombination breakpoint is uniformly distributed within this material. The waiting time going backward in time until k sequences have only k  1 ancestors in the population is exponentially distributed with intensity parameter k(k  1)/2, and the two sequences that have a common ancestor at that time are uniformly distributed among different pairs. This was first realized by Watterson (1975), and later developed into the theory of the coalescent by Kingman (1982).
The coalescent with recombination has further been investigated by Hudson and Kaplan (1985), Kaplan and Hudson (1985), Griffiths and Marjoram (1996, 1997), and Wiuf and Hein (1997, 1999).
The genealogy of a sample of sequences can be simulated by going back in time, waiting for what occurs first, a recombination or a coalescence, and then performing the appropriate operation on the set of ancestral sequences. Recombination increases the number of sequences carrying ancestral material by one, but does not increase the total amount of ancestral material. A coalescence decreases the number of sequences with ancestral material by one. It can increase the amount of material where recombination can occur, because coalescence can trap some nonancestral material (Figure 1). When any position on the extant sequences has found a MRCA, not necessarily the same ancestor, all segments with ancestral material spliced together constitute one sequence. Above this point, coalescence cannot reduce the amount of ancestral material and all that occurs is redistribution of ancestral material on different sequences by recombination and coalescence. Because the rate of coalescence is quadratic in the number of sequences, and the rate of recombination is at most linear, all positions eventually find a MRCA.
RESULTS
In this section we present simulated and mathematical results related to the MRCAs of a sample on k sequences (sections 16). We used an algorithm described in Wiuf and Hein (1999) to simulate sample histories. For each value of k = 2, 3, 5, 10, 25, 50, and 100 we simulated 2000 sample histories with recombination rate ρ = 50.
1. Definitions: We define a number of mathematical quantities that relate to the coalescent with recombination and to the results that we derive and discuss below.
Assume a sample of size k is given, with k possibly infinite. Let MRCA_{k}(q) denote the MRCA to position q in the sample of k sequences. The time until the MRCA_{k}(q) is distributed according to the coalescent process without recombination because one position cannot be subject to recombination. Further, let A_{k}(q) = 1 if there is a shift from one MRCA to another MRCA in position q, and A_{k}(q) = 0 otherwise; and let B_{k}(q) = 1 if there is a recombination breakpoint in position q within ancestral material, and B_{k}(q) = 0 otherwise. A_{k} stands for ancestor and B_{k} for breakpoint. We have A_{k}(q) = 1 iff the MRCAs to the left and to the right of position q are different, i.e., if MRCA_{k}(q  ε) ≠ MRCA_{k}(q + ε), provided ε is small. The definitions are illustrated in Figure 2. Both quantities A_{k}(q) and B_{k}(q) depend on the ancestral history of positions local to q only and not on the entire sequence history. Note that if A_{k}(q) = 1 then also B_{k}(q) = 1, but not necessarily the other way around. All recombination events do not necessarily result in a shift from one MRCA to another MRCA. Moreover, the distributions of the A_{k}(q)’s and B_{k}(q)’s, q ≥ 0 are invariant under translations along the sequences. As an example, (A_{k}(0), A_{k}(q)) is distributed like (A_{k}(p), A_{k}(q + p)); the distribution depends on the relative distance between positions only (q = q  0 = q + p  p) and not on the actual positions. This makes the two processes stationary processes (Daley and VereJones 1988).
On the basis of the definitions of A_{k}(q) and B_{k}(q) we define
The number R_{k}(ρ) was first studied by Hudson and Kaplan (1985), and S_{k}(ρ) by Griffiths and Marjoram (1997). We call S_{k}(ρ) the number of segments carrying ancestral material in the set of MRCAs. In light of the above discussion, this can be slightly misleading but is kept for matters of convenience.
In what follows we denote the length of a sequence by R ≡ ρ/2. Let
If the coalescent process with recombination is studied on a grid of points equally spaced with distance ε (in contrast to a continuum of points), the definition of L_{k}(ρ) would be L_{k}(ρ) = ε R^{R}^{ε1}_{i}_{=1} 1{MRCA_{k}(0) = MRCA_{k} (iε)}; that is, the number of times the MRCA_{k}(0) is visited moving along the sequences multiplied by the distance between the points.
We call MRCA_{k}(0) (or the MRCA to position 0), the leftmost MRCA, and L_{k}(ρ) the amount of material sharing MRCA with position 0.
2. Segment length: Hudson and Kaplan (1985) showed that the number, R_{k}(ρ), of recombination events within ancestral material until all positions have found a MRCA has expectation
We are interested in the sequence length between successive recombination breakpoints and the length between successive shifts between MRCAs. The above equations give us the expected number of each kind, recombination events/breakpoints and shifts between MRCAs, in sequences of length R ≡ ρ/2.
Denote by r_{k} the length between q = 0 and the first recombination point along the sequences, and by s_{k} the length to the first shift from one MRCA to another MRCA measured from q = 0 (Figures 2 and 3). We here assume that sequences are potentially infinite so that there always is a first recombination event and a shift between MRCAs. Because of the stationary property of the process, it follows that the expected value of r_{k}, given a recombination event happened in position q = 0, is
Similarly, one obtains the expected value of s_{k}, given that there is a shift in position q = 0 from one MRCA to another MRCA, by
The two expressions (6) and (7) hold for k = ∞ as well, yielding
In Wiuf and Hein (1999) the expected length of r_{k} is calculated
3. Number of MRCAs: The number of different MRCAs is upward bounded by S_{k}(ρ), the number of segments carrying ancestral material in the set of MRCAs and hence bounded in expectation by
4. The leftmost MRCA: Consider now the MRCA to position q = 0, MRCA_{k}(0). In the case k = 2 we can calculate the expected amount of positions sharing MRCA with position 0, L_{2}(ρ) [see (3)], as a function of R ≡ ρ/2, the sequence length. We find (see appendix)
Further, we find the following lower bound on the variance of L_{2}(ρ) (see appendix):
Combining the expression for the expectation of L_{2}(ρ) with the lower bound on the variance, we find that the normalized variable 2L_{2}(ρ)/log(1 + ρ) has expected value ∼1, but has a variance that increases without bound for increasing sequence length. Thus, there is very large variation in the amount of positions sharing MRCA with position 0 when sample size is 2.
For larger sample sizes, k > 2, we have P{MRCA_{k} (0) = MRCA_{k}(q)} ≤ P{MRCA_{2}(0) = MRCA_{2}(q)} (see appendix), and hence that
This bound is very crude. We have simulated the length, L_{k}(ρ), for samples of different sizes and found that for large sample sizes, the expected length is a slowly growing function of ρ (Figure 5), and for k > 2 it almost becomes constant. The difference in expected length between samples of size 10 and 100 is <5% within the range ρ is varied. This indicates a quick convergence of expected length for increasing sample sizes and fixed ρ.
The set of positions sharing MRCA with position 0 is (potentially) partitioned into several segments (Figure 3). Figure 6 shows the expected number of such segments on the leftmost MRCA. For small sample sizes there are less segments than for large sample sizes. Moreover, as shown in the previous figure, the expected length of positions sharing MRCA with position 0, E[L_{k}(ρ)], is smaller for large sample sizes than for small sample sizes. This supports the conclusion that the material sharing MRCA with position 0 is chopped into more segments for large sample sizes than for small sample sizes, and that these segments tend to be shorter for large samples than for smaller samples.
The histograms in Figures 7 and 8 show the amount of trapped material on the leftmost MRCA. Trapped material is between the segments sharing MRCA with position 0 (Figure 3). Figure 3 supports the conclusion that as sample size increases so does the chance of finding trapped material on the leftmost MRCA, i.e., the chance that the ancestral material is located on different segments increases. We conclude, as we did in Figure 6, that the chance of several segments being on the leftmost MRCA is highest for large samples.
5. Shared sequence ancestry: Consider two positions. As the distance between the positions gets larger, their genealogical histories become less correlated. In general, the correlation between genealogies might be measured in different ways according to what aspects of the genealogy are of interest. Kaplan and Hudson (1985) found that the covariance between the total branch lengths, G_{k}(q_{1}) and G_{k}(q_{2}), of the genealogies in positions q_{1} and q_{2} is about k/(4(k  1)R), where R ≡ ρ/2 is the distance between q_{1} and q_{2}. A similar result will hold for the tree heights of the genealogies.
If the positions are completely linked, the positions ancestral to q_{1} and q_{2} are on the same ancestral sequence. When recombination is present, the ancestral positions to q_{1} and q_{2} will not necessarily share sequence, but might be on different sequences. The time they share ancestral sequences is a measure of the correlation between the two genealogies. We define and discuss a notion of shared sequence ancestry in this context.
Let a sample of size 2 be given. Fix two positions, q_{1} and q_{2}, on the sequences s_{1} and s_{2} with distance R (recombination rate ρ = 2R).
Denote an ancestral state to the sample by a list ((x_{i}_{1}, x_{i}_{2})i = 1,..., n), where n is the number of ancestral sequences, and (x_{i}_{1}, x_{i}_{2}) denotes an ancestral sequence. The variable x_{ij} is ^{*} if position q_{j} on the sequence represented by (x_{i}_{1}, x_{i}_{2}) is nonancestral to any position in the sample, 0 if ancestral to q_{j} on both s_{1} and s_{2}, 1 if ancestral to q_{j} on s_{1} only, or 2 if ancestral to q_{j} on s_{2} only. The definition is illustrated in Figure 9. A presentday sample is represented by ((1,1), (2,2)). We say that the positions share sequence ancestry whenever the ancestral state is ((x_{i}_{1}, x_{i}_{2})i = 1,..., n) = ((1,1), (2,2)) or = ((0,0)). This implies that the positions ancestral to q_{1} and q_{2} on s_{1} share an ancestor, and at the same time the ancestral positions to q_{1} and q_{2} on s_{2} share an ancestor, possibly the same.
If the positions are completely linked (R = 0) the ancestral state is ((1,1), (2,2)) until the positions find a MRCA and the state becomes ((0,0)). If the positions are less linked, the ancestral state might be different from ((1,1), (2,2)). As a measure of shared sequence ancestry we take the expectation of the time T_{S} spent in the state S = ((1,1), (2,2)) compared to the expectation of the time T_{j}, j = 1, 2 until a position finds a MRCA, i.e.,
Similarly, we can calculate the shared sequence ancestry given that the two positions find a MRCA at the same time, i.e., given MRCA_{2}(q_{1}) = MRCA_{2}(q_{2}) or T_{1} = T_{2}. We find
The value of E[T_{j}T_{1} = T_{2}] is ∼3/ρ for large ρ’s. The time spent in S is only about ^{1}/_{3} the total time, so that in general several events happen before the MRCA and not just a single coalescent event.
6. Tree heights and branch lengths: Let T_{k}(q) denote the time until a MRCA in position q in a sample of size k, and G_{k}(q) the total branch length of the genealogy in position q. Because one position cannot be subject to recombination these two variables depend on only the coalescent process and not the recombination process. The expectation of T_{k}(q) is
However, the distribution of T_{k} = max_{0≤}_{q}_{≤ρ/2} T_{k}(q) is highly dependent on the recombination rate ρ. The variable T_{k} is the time until all positions along the sequences have found a MRCA. Griffiths and Marjoram (1997) found that the expectation of T_{k} is bounded:
This bound is uniform in k, and linear in ρ. We have simulated T_{k} to see how good this bound is (Figure 10). E[T_{k}]  E[T_{k}(0)] seems to converge toward a logarithmic limit.
Similarly, the expectation of G_{k}(q) is
DISCUSSION
We have discussed properties of the ancestry of k sequences sampled from the coalescent process with recombination. We have done so mainly by simulations of sample histories. A number of variables derived from the genealogies were observed, and the expectations over 2000 simulations were calculated.
Each of these variables describes an aspect of the ancestry of a sample. One should in general be cautious in interpreting the behavior of a process from expectations only: The variation in the process is ignored when relying on expected values, and it is not guaranteed that the expected value represents a “typical outcome” of the variable.
However, a comparison of expected values for varying sample sizes, k, and varying recombination rates, ρ, gives an idea of how the variables depend on k and ρ and thereby an idea of the amount of information in the size of the sample.
In particular we were interested in the leftmost MRCA, i.e., the MRCA to position 0 in the sample. The length of material sharing MRCA with position 0 decreases with increasing sample size toward a limit (sample size ∞). At the same time, the number of segments on the leftmost MRCA increases with increasing sample sizes.
For samples of size 2, we discussed a concept of sharing sequence ancestry between two positions (or loci). It was shown that, even when the two positions share a MRCA, the proportion of time they share ancestral sequences is short. For an increasing recombination rate the positions share sequence ancestry in ⅓ of the time until a MRCA.
Finally, simulations indicated that the expectations of the variables T_{k}  T_{k}(q) and G_{k}  G_{k}(q) are bounded in k (sample size) by a logarithmic function of ρ. The variable T_{k} is the time until all positions in the sample have found a MRCA, and G_{k} is the maximum of the total branch length of the genealogy over all positions. The bounds revived theoretically are far higher than the simulated curves.
It is interesting that the structure of a MRCA to a sample of sequences converges very quickly toward a limit structure (in expectations). In all figures, the difference between the simulation results for samples of size 10 (in a few cases 25) and larger sample sizes is close to zero. This indicates that the structure of a MRCA of a sample of size 10 has identical structure to a MRCA of the whole population.
Moreover, the expectation of the waiting time until all positions in a sample of size 10 have found a MRCA is about the same as the expectation of the waiting time until all positions in the whole population have found a MRCA. This finding is very similar to a result about the coalescent without recombination: The distribution of the waiting time until a sample of size 10 has found a MRCA is almost distributed like the same waiting time until the whole population has found a MRCA.
The explanation for this seems to be the following: Consider a large sample. The time during which there are many ancestors to the sample is considerably smaller than the time during which there are a few ancestors only. The rate of recombination is kρ/2 if there are k sequences, and the rate of coalescence is k(k  1)/2. If k is much larger than ρ, most events will in the beginning be coalescence events. Thus, the time from the present until the whole sample has been reduced to a small number of ancestors by coalescence events will be distributed similarly to the time until a large sample is reduced to a small sample in the coalescent without recombination. It is, therefore, the size of the minor number of ancestors that determines the structure of the variables we have discussed.
However, it is surprising that the convergence in sample size k seems almost uniform in ρ. The reason for this might be that the range within which ρ has been varied is too narrow to detect the dependence on ρ.
As an example, consider a large sample of human DNA sequences. Assume that the probability r of a recombination event between two nucleotides per generation per sequence is 10^{7} and that the effective population size of the human population is 2N = 10^{4}. If the number of nucleotides is L = 10^{4} (typical gene length), then ρ = 4NLr = 2 · 10^{4} · 10^{4} · 10^{7} = 20 and sequence length is R = ρ/2 = 10. In this case, there are ∼15 different MRCA consisting of ∼21 segments in total [see (5) and Figure 4]. Each segment will on average be E[s_{k}A_{k}(0) = 1] · L/R = ^{1}/_{2} · 10^{4}/10 = 500 [see (7)] nucleotides long and each MRCA 700 nucleotides long (L/15 = 10^{4}/15). Focus now on nucleotide 500 in the sequences. The length of the sequences to the right of the nucleotide is R/2 = 5. From Figure 5 we find that the expected number of positions sharing MRCA with nucleotide 500 is about 0.75 · 10^{4}/10 = 750. Similarly, ∼750 nucleotides to the left of number 500 will share MRCA with nucleotide 500; in total, 1500 nucleotides.
The expected time back until all nucleotides have found a MRCA is ∼5 · 2N = 50,000 generations (Figure 10). Counting 1 generation as 20 years, this is about 1 million years ago, whereas a random spot has average time to the MRCA of 40,000 years.
APPENDIX
Numbers refer to sections in results.
2. Following Daley and VereJones (1988), the expressions E[r_{k}B_{k}(0) = 1] and E[s_{k}A_{k}(0) = 1] are understood in the following sense. Extend the sequences by a small interval I of length ∂ (in units of expected number of recombinations per sequence per 2N generations) to the left of the position 0. Let β_{k}(∂) denote the event in which there is a least one recombination event within the sequence interval I in the history of the sample of size k. Similarly, let α_{k}(∂) denote the event in which there is at least one shift in MRCAs within the interval I. Then
Consider now the case k = ∞. The chance that there is at least one recombination event in the history of the sample in any interval of the sequences is one, independent of the size of the interval considered. That is, P(β_{k}(∂)) = 1 for all ∂ > 0 and P(R_{∞}(η) ≥ 1) = 1 for all η > 0. Thus,
To prove E[s_{∞}A_{∞}(0) = 1] = ^{1}/_{2} we show that E[S_{∞} · (ρ)] = 1 + ρ. Unfortunately, S_{k}(ρ) does not converge toward S_{∞}(ρ) in any regular way. Label all sequences in an infinite sample by numbers. Let S_{k}^{*}(ρ) be the number of recombination breaks, q, in the history of the first k sequences such that MRCA_{∞}(q  ε) ≠ MRCA_{∞}(q + ε) provided ε is small. Clearly
4. By definition of L_{2}(ρ) we have L_{2}(ρ) ≤ ρ/2. Thus, by dominated convergence
To prove P{MRCA_{k}(0) = MRCA_{k}(q)} ≤ P{MRCA_{2}(0) = MRCA_{2}(q)} we proceed as follows: If MRCA_{k}(0) = MRCA_{k}(q), then the last event must be a coalescence between two sequences, both carrying ancestral positions to 0 and q. Let S denote the state consisting of two ancestral sequences and both positions 0 and q ancestral to the sample. Let s denote the time until state S is entered for the first time, and let F(·) denote the probability distribution of s. We have
5. The state space of the Markov chain consists of all possible ancestral configurations ((x_{i}_{1}, x_{i}_{2})i = 1,..., n) to the sample. The variables x_{ij} are defined in 5. The state ((0,0)) is absorbing. This is very similar to a Markov chain in Simonsen and Churchill (1997) describing a twolocus model with sample size 2. The difference is that we distinguish between some states, e.g., ((1,1), (2,2)) and ((1,2), (2,1)), while Simonsen and Churchill (1997) do not. The transition probabilities are given by expressions similar to expressions in Simonsen and Churhill (1997), and follow from the structure of the twolocus model. For example, the probability of going from state ((1,2), (1,2)) to state ((1,^{*}), (^{*},2), (2,1)) is
To obtain the expectations E[T_{s}T_{1} = T_{2}] and E[T_{j}T_{1} = T_{2}], first note that for states S′ and S^{2} we have
6. Similar to Griffiths and Marjoram (1997, Theorem 3.1), we find
Acknowledgments
We thank Mikkel Nygaard Hansen for help with implementation of the simulation program. Bernt Guldbrandtsen is thanked for reading and commenting on the manuscript. J.H. was supported by Danish Research Council grant SNF 9401631.
Footnotes

Communicating editor: R. R. Hudson
 Received February 10, 1998.
 Accepted November 30, 1998.
 Copyright © 1999 by the Genetics Society of America