THE early reports of mitochondrial DNA (mtDNA) haplotype diversity among humans indicated a recent large increase in population size starting from a small number of founders (Cannet al. 1987; Vigilantet al. 1991). This pattern, which is marked by an increase in rare variants, has not subsequently been observed for nuclear loci. Instead, nuclear loci show an excess of polymorphic sites segregating at intermediate frequencies (Hey 1997). This is illustrated by Tajima's (1989) statistic, D, which is positive, although not significantly so, for four extensive nuclear DNA sequence data sets (Hardinget al. 1997; Clarket al. 1998; Zietkiewiczet al. 1998; Harris and Hey 1999). In addition, single nucleotide polymorphism (SNP) data do not show evidence of population growth (Nielsen 2000). Here, I consider a subdivided population model, which allows for a single change in the effective size of the population, and apply the results to restriction fragment length polymorphism (RFLP) data from a worldwide sample of humans (Bowcocket al. 1987; Matulloet al. 1994; Poloniet al. 1995). The data provide strong evidence for a change in demography: a decrease in effective size. This is most likely due to a shift from a more ancient subdivided population to one with less structure today.

Numerous models of subdivision with migration have been proposed to explain patterns of genetic variation among local populations. The most important distinction among these is that some models, such as Kimura and Weiss' (1964) stepping-stone model, generate isolation by distance (Wright 1943), whereas others, such as Wright's (1931) island model, do not. The model considered here is of the latter type. Wakeley (1998) studied the ancestral genealogical process for samples from a D-deme version of the finite island model, in which the number of demes is assumed to be large. In the finite island model (Maruyama 1970; Latter 1973), all demes exchange migrants regardless of their geographic location. The single-generation probability of migration is assumed to be the same between each pair of demes. Specifically, m is the probability that a gene came from one of the other D − 1 demes to its present one in the previous generation. This model is realistic for species in which individuals choose either to stay put or to move to a location picked uniformly from the entire species range. In certain situations, it may also apply when individuals choose either to migrate a great distance or not at all.

When the number of demes is large, the genealogical history of a sample taken from such a population can be divided into two parts, which I will call the scattering phase and the collecting phase; see Figure 1. The scattering phase is the very recent history of the sample, during which coalescent and migration events bring it rapidly to the start of the second, and much longer, collecting phase of the history. The collecting phase begins when each ancestral lineage is in a separate deme and is characterized by migration events that move lineages from deme to deme, the great majority of which are unoccupied. Eventually, a migration event will place a pair of lineages into the same deme, at which time they can coalesce. The ability to break the genealogy into these two parts, and to consider them separately, depends only on the number of demes being large. When this is true, the time spent in the scattering phase can be ignored because the collecting phase dominates the history (Wakeley 1998). A similar separation of time scales has been found in the study of genealogies under partial selfing (Nordborg and Donnelly 1997).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

A hypothetical genealogy of a sample from d = 3 demes under the model. In this case there was just one migration event during the scattering phase, and this was in the sample from deme 2. Thus , , , and n′ = d + 1 = 4.

Consider first the equilibrium version of this island model. We assume that each of D demes is of constant diploid size N, but any results will also apply to haploids if N is placed by N/2. This model is characterized by two parameters: θ = 4NDu, in which u is the neutral mutation rate, and M = 2NDm/(D − 1). When D is large, M becomes simply 2Nm. A sample of n items, taken from d different demes, which are arbitrarily labeled 1 through d, is denoted n = (n₁, n₂, …, n_d), where n_i is the sample size from the ith deme. Of course n=Σi=1dni . This notation for the sample differs from that used in Wakeley (1998). We assume that θ is finite even though D and N are large. Thus, we have D ⪢ n in addition to the usual coalescent condition, N ⪢ n. Simulations show that D does not have to be terribly big for the infinite-demes approximation to hold (Wakeley 1998).

The scattering phase of the genealogy takes the singledeme sample, n_i, through a series of migration and coalescent events to a new configuration where each remaining lineage is in a separate deme. The number of these recent ancestral lineages, n′_i (1 ≤ n′_i ≤ n_i), depends on the sample size and the migration rate. It is important to note that each of these n′_j lineages is in a separate deme at the end of the scattering phase. The distribution of n′_i is given by P[ni′∣ni]=Sni(ni′)(2M)ni′(2M)(ni) (1)(Wakeley 1998; Equation 32), where x_(r) = x(x + 1) … (x + r − 1), and S⁽ⁱ⁾_j is an unsigned Stirling number of the first kind. Abramowitz and Stegun (1964) define Stirling numbers and list many of their properties. The notation |s(i, j)| was used instead of Sj(i) in Wakeley (1998). Slatkin's (1989) Equation 2 for the probability that n alleles are descended from nonimmigrants is a special case of Equation 1 above. (Note that there is a typographical error in Slatkin's equation; there should be a 2Nm in the numerator.)

Thus, the distribution of the number of ancestral lineages of each deme's sample at the start of the collecting phase is the same as the distribution of the number of alleles in the Ewens sampling formula (Ewens 1972; Karlin and McGregor 1972) but with mutation replaced by migration. There is one major difference: the n′_i items at the start of the collecting phase represent exactly that number of ancestral lineages, whereas the alleles counted by the Ewens sampling formula do not. Let n′=(n1′,n2′,…,nd′) and n′=Σi=1dni′ . Because events in different demes are independent, the joint distribution of all the n′_i is given by P[n′∣n]=∏i=1dP[ni′∣ni]. (2)The collecting phase of the history then begins with n′ (d ≤ n′ ≤ n) lineages, and each of these is in a separate deme. Equation 2 above follows from Equation 36 in Wakeley (1998).

The genealogical history of these n′ ancestral lineages during the collecting phase is a coalescent (Kingman 1982a, b) when time is measured in units of twice the effective population size Ne=ND[1+12M] (3)(see appendix a). Nei and Takahata (1993) give a similar formula for the special case of two sampled sequences under the finite island model, which they noted was implicit in earlier work, and which quickly converges on Equation 3 as the number of demes increases. Thus, the duration of the collecting phase is at least of order ND generations, whereas the duration of the scattering phase is of order N generations. This is why the time of the scattering phase can be ignored when D is large.

Because the collecting phase is a coalescent, it is relatively simple to allow for changes in effective population size. Here I consider the possibility of a single abrupt change in demography at generation t in the past. Thus, the collecting phase itself is divided into two parts. Before time t, the population is assumed to have been of effective size N_eA = N_AD_A[1 + 1/(2M_A)], which can be compared to Equation 3. Thus, the change in demography might have been caused by a change in deme size, in the number of demes, or in the migration rate. Part of the flexibility of the model is that the change might have been any combination of these. If the ancestral population was panmictic we have the demographic model that Takahata (1995) proposed for human history, but without extinction and recolonization.

This nonequilibrium model necessitates two parameters in addition to θ and M. The first is the time of the event, T = t/2ND), measured in units of 2ND generations. The second is R = N_eA/(ND), the ratio of the ancestral effective size to the current total population size. Thus, values of R/[1 + 1/(2M)] greater than one indicate that the ancestral effective size was larger than the current effective size. Given an estimate of R and assuming that ND = kN_AD_A, where k is some constant, we can retrieve M^A=1∕[2(kR^−1)] . This is of considerable interest below in the application of this model to human history. Further, given estimates of both R and M and assuming M_A = kM, we could estimate (N_AD_A)/(ND) as R^∕[1+1∕(2kM^)] , which in conjunction withû could be used to estimate θ_A.

Neutral mutations are modeled as a Poisson process along the branches of the genealogy (Hudson 1983, 1990). When a mutation happens it is assumed to occur at a previously unmutated site. This assumption is shared by Watterson's (1975) infinite sites model without recombination and Kimura's (1969) model, which assumes free recombination between sites. We apply these models to genetic data when the mutation rate per site is so small that multiple changes at any one site are unlikely. Watterson (1975) and others have studied the distribution of segregating sites under the neutral infinite sites model, and Sawyer and Hartl (1992) have done so for both selection and neutrality under Kimura's (1969) model. Under neutrality, both models make the same prediction about the expected number of segregating sites in a sample, either taken as a whole or broken into categories by frequency (Hudson 1983; Tajima 1989; Sawyer and Hartl 1992; Fu and Li 1993). In other words, expectations derived under particular assumptions about the recombination rate are valid for all rates of recombination. Here, I assume Watterson's (1975) model. The effect of linkage and recombination is on the variances, and these are largest under when the recombination rate is lower.

Under these mutations-at-unique-sites models, each segregating site will divide the sample into two groups: one that has inherited the mutant base and one that still shows the ancestral base. Let i = (i₁, i₂, …, i_d) be the counts of mutant bases at a single polymorphic site in the sample n = (n₁, n₂, …, n_d). Further, let z_i be the number of sites in the sample that show the pattern i = (i₁, i₂, …, i_d). In this article, I calculate L(M,R,T;i)=Pr(i∣site is variable;M,R,T), (4) the likelihood of a particular mutant site pattern given that a site is polymorphic. I do this by first deriving E(z_i), the expected number of sites that show the pattern i = (i₁, i₂, …, i_d), for a sample of sequences under Watterson's infinite sites model. As discussed above, this expectation will hold for arbitrary rates of recombination. Sawyer and Hartl (1992) showed that with free recombination between sites, z_i is Poisson distributed with parameter E(z_i), so L(M,R,T;i)=E(zi)ΣiE(zi). (5)The denominator in Equation 5 is simply equal to the expected number of segregating sites in the sample. As E(z_i) is linear in θ, the likelihood depends only on the demographic parameters, M, R, and T, and not on the mutation parameter, θ. Unless an outgroup sequence is available, it is impossible to distinguish the mutant from the ancestral base. Thus, except for the special case of i_j = n_j/2 for all j, the likelihood of a particular site pattern is L(M, R, T; i) + L(M, R, T; i^c), where i^c = (n₁ − i₁, n₂ − i₂, …, n_d − i_d). With this modification, Equation 5 is suitable for the analysis of SNP and RFLP data. If different sites can be assumed to be independent of each other, as they can for the human RFLP data analyzed below, the likelihood of the entire data is the product of the likelihoods for individual sites.

APPLICATION TO HUMAN RFLP DATA

The method was applied to some data derived from the extensive RFLP surveys of Bowcock et al. (1987), Matullo et al. (1994), and Poloni et al. (1995). These data overlap with those analyzed by Nielsen et al. (1998). I considered data from four localities, which I call demes: Japan, New Guinea, Senegal, and Italy. Joanna Mountain kindly supplied these data, which consisted of n_j = 20, for j = 1, 2, 3, 4, at 45 independent loci. Due to the computational burden of calculating the likelihood 6, I made a smaller data set of n_j = 7, for all j, by randomly sampling haplotypes within demes (without replacement). One locus was lost because it was not polymorphic in the subsample. This left 44 independent loci for the analysis, and these data are shown in Table 1. I also performed all of the analyses described below on subsamples of size n_j = 3 and n_j = 5, and found nearly identical results. The full data sent to me by Joanna Mountain included two more African sample locations: Central African Republic and Zaire. I excluded these two samples for reasons discussed below. However, I also analyzed just the three African samples and the results were identical in overall pattern to what is reported here for the four widely separated samples. The biggest difference was that the estimate of the current migration rate was almost five times larger.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

A contour plot of the likelihood of the data in Table 1 over a range of R and T, when M = 1.83 (see text). Any waviness in the surface results from the fact that the likelihood was calculated at a finite number (25 × 25 = 625) of points.

For the data in Table 1, the maximum-likelihood estimate of the migration rate was M^=1.80 , with a 95% confidence interval of (1.22, 2.83). Thus, current levels of gene flow among human subpopulations appear to be relatively high. Estimates of M varied only slightly over the entire range of R and T. For example, the equilibrium migration model—i.e., setting R equal to 1 + 1/(2M)—gave M^=1.83 . Therefore, to save time, a single value of M = 1.83 was used to compute the likelihood surface for R and T shown in Figure 2. I confirmed that this had little effect on the results by separately maximizing the likelihood, over M, at nine evenly distributed points on the surface, and comparing these results with those obtained using M = 1.83. Figure 2 is a contour plot of the likelihood of the data in Table 1 over a broad range of both R and T. Each contour in Figure 2 represents a 2-unit change in log likelihood, so that the approximate joint 95% confidence region for R and T is enclosed within the first contour.

A number of conclusions can be drawn from this figure. First, it is clear that small values of R and T (lower left in Figure 2), which could represent recent population expansion, are strongly rejected. Next, the present model reduces to an equilibrium model when R = 1 + 1/(2M), that is when the effective size was the same before and after T. For M = 1.83, this gives R = 1.27 for all T, a value that falls between the third and the fourth contours in Figure 2. These contours define a broad flat area more than 6 log-likelihood units below the maximum, which occupies much of the graph. Thus, equilibrium migration (i.e., no change in demography) is also strongly rejected. As T either increases or decreases for a given finite value of R, the likelihood returns to this broad flat area, which is consistent with an equilibrium model. This makes intuitive sense because the change in demography will have little effect if it is either extremely recent or extremely ancient.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

The likelihood surfaces (a) for R and (b) for T when R = ∞. Only values within two log-likelihood units of the maximum are shown. Thus, a traces the ridge of the surface shown in Figure 2, and b shows a slice of the surface but for R much greater than any in Figure 2.

As Figure 2 implies, there is no upper bound for R. The 95% confidence interval for R is (5.0, ∞), and in fact the likelihood of the data continues to increase as R grows. As R increases, E_n′(z_i′) and E_n′(S)—see Equations 7 and 15—become linear in R as well as θ. For instance, E_n′(S) converges on θR∑n″=2n′gn′n″(τm)∑i″=1n″−11i″ (16)as R goes to infinity. Both θ and R then drop out of the likelihood. When R is large, our power to say anything about its exact value diminishes. Using the large-R limiting expressions gives a maximum log-likelihood value of −309.17535. Figure 3a, which traces the ridge of maximum likelihood on the surface in Figure 2, shows that the likelihood is within 1% of its maximum value by the time R reaches 10³. Thus, while we cannot reject huge values of R, we might argue that the evidence for R being greater than about 10³ is not strong. In sum, there is clear evidence for at least a fivefold difference in recent vs. ancient effective population size in humans. If we assume that N_AD_A = ND, then from the confidence interval for R we find that the ancestral migration rate among human demes, M_A, must have been 0.125 or smaller.

The timing of this change in demography is also difficult to establish precisely. Figure 2 shows that the lower bound of the 95% confidence interval for T does not depend much on R, but that the upper bound is positively correlated with it. If we were to draw a line along the ridge of the maximum-likelihood surface in Figure 2, it would run from T = 0.81 at R = 10 to T = 3.42 at R = 10⁶. The more ancient the event was, the more extreme it has to have been to explain the data. Using the limiting expressions for E_n′(z_i′) and E_n′(S). gives a 95% confidence interval for T as (0.2, ∞). Figure 3b shows that the likelihood increases with increasing T, but that the surface becomes quite flat above T = 2.0. Translating the confidence interval for T into years requires values for the recent effective size and the generation time for humans. Some typical values are N_e = 10⁴ and g = 15 yr for the generation time (Takahataet al. 1995; Mountain 1998). Assuming that N_e is given by Equation 3, that M = 1.8, and using the lower bound of 0.2 for T, we estimate that this demographic event could have occurred as recently as 2N_eTg/[1 + 1/(2M)] ≈ 47,000 years ago, but might have been much more ancient.

• Received May 6, 1999.
• Accepted August 12, 1999.

• Copyright © 1999 by the Genetics Society of America