Abstract
A new method for the estimation of migration rates and effective population sizes is described. It uses a maximum-likelihood framework based on coalescence theory. The parameters are estimated by Metropolis-Hastings importance sampling. In a two-population model this method estimates four parameters: the effective population size and the immigration rate for each population relative to the mutation rate. Summarizing over loci can be done by assuming either that the mutation rate is the same for all loci or that the mutation rates are gamma distributed among loci but the same for all sites of a locus. The estimates are as good as or better than those from an optimized FST-based measure. The program is available on the World Wide Web at http://evolution.genetics.washington.edu/lamarc.html/.
SEVERAL methods for the estimation of migration rates between subpopulations have been proposed. We can subdivide them into two very different approaches: (1) marking individuals and tracking their individual movements and then extrapolating these individual movements to migration rates; or (2) surveying genetic markers in the populations of interest and calculating a migration rate from allele frequencies or sequence differences. Of course, one should be aware that these two approaches do not estimate the same quantity: approach 1 estimates the actual “instantaneous” migration rate of individuals, whereas approach 2 reflects an average over a time period whose length is determined by the rate of mutation per generation of the locus under study or by the time scale of genetic drift. The genetic approach tends to gives a lower estimate than the individual migration rates approach because the method is looking at changes that become established in the subpopulation gene pool.
Current estimation methods for genetic data are methods such as those related to FST (e.g., Slatkin 1991; Slatkin and Hudson 1991), the rare allele approach of Slatkin (1985), a maximum-likelihood method using gene frequency distributions (Rannala and Hartigan 1996; Tuftoet al. 1996), and approaches based on coalescent theory (Kingman 1982a,b), such as the cladistic approach of Slatkin and Maddison (1989), the method outlined in Wakeley (1998), and maximum likelihood using coalescent theory (Nath and Griffiths 1993, 1996). We describe here a new method using a maximum-likelihood- and coalescent theory-based approach.
Our method integrates over all possible genealogies and migration events. This should be superior to pairwise estimators such as those using FST or related statistics (cf. Felsenstein 1992b) and also more powerful than the cladistic approach of Slatkin and Maddison (1989), which needs to know the true genealogy. Our approach estimates similar parameters as the program of Bahlo and Griffiths (http://www.maths.monash.edu.au/~mbahlo/mpg/gtree.html/), which is based on the work of Griffiths and Tavaré (1994) and Nath and Griffiths (1996). It differs in the way we search the genealogy space, and our methods support mutation models for different types of data: the infinite allele model for electrophoretic markers, a one-step model for microsatellites, and a finite-sites model for nucleotide sequences. The sequence model is more useful than the infinite sites model, which forces the researcher to discard data when there are multiple mutations at the same sites.
MODEL
We propose a method to make a maximum-likelihood estimate of population parameters for geographically subdivided populations. The general outline of such estimates involves extending coalescence theory (Kingman 1982a,b) to include migration events. Migration models with coalescents were first developed by Takahata (1988) and Takahata and Slatkin (1990) for two gene copies and discussed more generally for n gene copies in two populations by Hudson (1990), with generalization to multiple populations and different models of migration by Notohara (1990). Nath and Griffiths (1993, 1996) used this migration-coalescent process for maximum-likelihood estimation of one parameter, the effective number of migrants
—Two-population model with four parameters: m1 is the migration rate per generation from population 2 to 1, m2 from population 1 to 2, and 
is the effective population size.
Our migration model consists of two populations. These have effective population sizes
Because we observe only genetic differences and do not know the mutation rate μ, we must absorb μ into the parameters so that we use as our parameters
The likelihood formula for this family of methods (Felsenstein 1988, 1992a; Kuhner et al. 1995, 1998) is based on the product between the genealogy likelihood Prob(D|G) (e.g., Felsenstein 1973; Swoffordet al. 1996) and the prior probability of the coalescent genealogy
Takahata and Slatkin (1990) found that Kingman’s results could not be extended to cases with geographical structure. They were unable to obtain the distribution of times to coalescence in the presence of migration. We have avoided this difficulty by having the genealogy G specify not only the coalescences but also the times and places of migration events. With this information in G, its probability density becomes for two populations a product of terms for the intervals between events in the genealogy
If we extend the parameter estimation from a single locus to many unlinked loci, we face the additional problem that, although μ may be constant per locus, it can vary between loci. Neglecting variation of μ, combining
Note that the gamma distribution described here is the distribution of μ across loci, not the more commonly used gamma distribution of rates across sites within a locus. If needed the latter can be approximated by Hidden Markov model methods in the case of DNA sequences.
IMPLEMENTATION
To calculate the likelihood
The parameter estimation is simple in principle: we have to find the maximum of the likelihood function, which has five or six dimensions, depending on whether μ is allowed to vary between loci. This is done by a modified damped Newton-Raphson method (Dahlquist and Björk 1974) using explicit equations for the first and second derivatives.
The genealogy sampler is of more interest. A large sample of genealogies G1, G2,...,Gm is drawn from a distribution whose density is proportional to Prob(D|G)
In the program Coalesce (Kuhneret al. 1995) a region of internal nodes in the genealogy is erased and rebuilt according to the coalescent model. This scheme is not applicable to migration models. We adopt a new scheme that is much more general and that can be used in any extension to the coalescent model. The first genealogy for each locus is constructed with a UPGMA method (as implemented by M. K. Kuhner and J. Yamato in Felsenstein 1993), and the minimal necessary migration events are inserted using a Fitch parsimony algorithm (Fitch 1971). The times for coalescent events or migration events on this first genealogy are calculated with (3) using a uniform random number for pj and solving for uj, and parameter values, which are either guessed by the researcher or calculated using an FST-based method (see appendix).
The genealogies are sampled as follows (cf. Figure 2): A coalescent node or tip z on this genealogy is chosen at random, the lineage below it is dissolved, and the node is used as the starting point to simulate the ancestry using a migration-coalescent process through a series of time intervals (Figure 2). The other lineages do not change and are taken as given and form the partial genealogy Gp. For the further description of the rearrangement process we use “up” and “top” for being or moving closer to the tips of the genealogy and “down” and “bottom” for being or moving toward the root of the genealogy. Starting at the chosen node z (see Figure 2), we draw a new time interval by solving for uj in
—Transition from an old genealogy Go to a new genealogy Gn. Dotted lines are the times of coalescences or migrations. (A) Genealogy Go with migration. The black bar marks a migration from the white population to the black or vice versa; z is the node to be picked. (B) Partial genealogy Gp after drawing the coalescent node z at random and dissolving the branch to the next coalescent below. (C) Simulation of the coalescent with migration. One possible outcome with three consecutive steps is shown: (1) using Equation 7 a new time interval is drawn and a migration event from white to black is also drawn (Equation 8) and the lineage is extended down to that event; (2) a new time interval is drawn: it extends too far back, so the lineage advances down to the time j; (3) a new time interval is drawn with a coalescent event at its bottom end. The process stops at time k. (D) the final configuration Gn.
In cases where all the lineages have not coalesced by the time of the root of the partial genealogy, simulation on the active lineage and the bottom lineage of Gp continues until the lineages coalesce (Figure 3). When simulating on two lineages we draw new time intervals using (3) instead of (7). There is one exception to this rule at the bottom of the genealogy: if by chance the dissolved lineage is the bottommost lineage of one of the populations on Go, we keep simulating with a single active lineage below the root of Gp until the former root of Go is reached and then start to simulate on both of the two remaining lines. The lineage between the root of Gp and the old root on Go on the partial genealogy has a known history, so there is no need to simulate this history again.
—Simulation on one or two lineages at the bottom of a genealogy. Striped lineages are active. On these lineages the migration-coalescent process is used to find new time intervals and events. (A) old genealogy Go, with 1 indicating the branch to be cut to get genealogy B, and 2 indicating the branch to be cut to get C. (B) The root R is below the cut-point D, and so below the root we need to calculate the migrations and coalescences for two lineages, because there is no previous information present. (C) In this genealogy a single lineage simulation is needed, using the information present between R and D. Below D we need to simulate on both remaining lineages.
A change to the newly found genealogy Gn is accepted with probability
This results in cancellation of the
If the genealogy is accepted, it is the starting point for a new cycle; otherwise the previous one continues to be used. After having sampled a large number of genealogies we estimate the parameters by maximizing (6). In theory a single long Markov chain would be enough to find the genealogies that contribute most to the likelihood, but if we start from a very bad
Simulation of 100 single-locus datasets with 25 sampled individuals for each population and 500 or 1000 bp, respectively, for different known parameter combinations
Influence of number of sites and number of loci on parameter estimates
SIMULATION RESULTS
We have simulated genealogies with known true parameters
For the comparison with an FST-based estimator (see appendix) we used a symmetrical parameter setup for the simulations, Θ1 = Θ2 and γ1 = γ2. However, the parameters that were estimated were not constrained to be symmetric.
Our simulations generally recover the true ΘT very well (Table 1). With 1000 bp the means are better than with 500 bp. The means for γ are not close to the true γT. With low ΘT values γ is grossly overestimated. The upward bias in γ is huge with 500-bp data and shrinks when longer sequences are used (Tables 1 and 2).
Sequences in two populations can be similar to each other due to migration or due to low variation caused by a small Θ. The simulated datasets with low Θ sometimes show by chance no variation at all in short sequences. In these cases one would expect the estimate of Θ to be 0.0 and the estimate of
The number of cases in which the FST-based estimator (Table 3) is usable is often dramatically lower than that with the ML estimator (Table 1; see also appendix). For example, with ΘT = 0.001 and γT = 10 the FST-based estimator could be used in only 45 out of 100 runs. This behavior seems not to improve with longer sequences. Increasing the number of sites and the number of loci helps to reduce the variance of the ML estimates (Table 2 and Figure 4). Estimates of γ from datasets with very long sequences are biased and are smaller than γT. This is the result of low acceptance rates during a run, when the starting genealogy is so well defined by the data that in our simulation runs the chains were too short and too few different migration scenarios were tried. The estimation of nonsymmetrical true parameters works quite well (Table 4). We encounter similar biases as with symmetrical parameters: with low Θ the γ estimates are biased upward and with high Θ they are biased downward.
FST-based estimator with three parameters (Θ ≡ Θ1 ≡ Θ2, M1, M2, see also appendix)
The assumption that mutation rate varies according to a gamma distribution adds more noise to the estimation. All biases are similar to the ones shown in the tables. The estimation of the shape parameter 1/α is certainly not very precise with only a few loci (Figure 5). One of the runs shown in Figure 5 obviously contains almost no information about 1/α and its log-likelihood curve is nearly flat, so that almost any value of 1/α can be accepted for this dataset. The increase of the log-likelihood values with very small 1/α and far from the maximum is due to imprecisions in the calculation using the discrete gamma approximation. In evaluations using a more exact but much slower quadrature scheme (Sikorskiet al. 1984) the peaks of the log-likelihood curves are at similar values, but log-likelihood for very low 1/α values are monotonically decreasing and not, as with the discrete gamma approximation, slightly increasing.
DISCUSSION
The estimation of the migration parameter γ seems to be rather imprecise, and this is true for both migrate and FST (e.g., Tables 1 and 3). The standard deviations for γ are large for single-locus estimates. More sites per locus give better estimates but the improvement of the variance is small, and given the obstacles of sequencing long stretches of DNA for several individuals one might better invest in the investigation of an additional unlinked locus. Each new locus has its own genealogy that is not correlated with that of the other loci and so increases the power of estimating parameters more, as we can see in Table 2.
There is considerable bias in the estimates of γ when datasets either have no or little genetic variation or have very high variation. Datasets with almost no variation inflate
In single runs we can also see a “fatal attraction” to 0.0 if the true parameters are very close to 0.0. The genealogy sampler is using the current parameter values to propose coalescences and migration events. If one or more of these parameters are very close to 0.0, then it can be seen by inspecting (2) that as a result our procedure will either most often propose an immediate coalescent event if one of the Θ is close to 0.0 or rarely propose an immigration event with very low γ. If a chain never proposes any instances of an event, its rate of occurrence will be estimated to 0.0, and this situation may persist indefinitely, a fatal attraction. Even if the parameters are prevented from becoming exactly zero this means that it may take large amounts of simulation to escape these values.
Because long sequences define the genealogy very well, our Markov chain Monte Carlo sampler will often reject newly proposed genealogies and only a few different genealogies are sampled for the parameter estimation. Therefore, the program does not readily explore the whole migration-genealogy space and tends to stick to the good starting genealogy. This generates a bias in
—Example of estimation of population parameters for two populations with 1, 3, and 10 loci. The graphs are cross-sections through the parameter space passing through the peak of the likelihood surface. The axes are on a logarithmic scale. The dashed lines give the true values of the parameter values: Θ1 = 0.05, Θ2 = 0.005, γ1 = 10.0, γ2 = 1.0. The gray area is an approximate 95% confidence set based on a likelihood-ratio test criterion defined as the range of parameter values with log-likelihood values equal to or higher than the maximum -0.5 × 9.4877 
.
Summing over loci assuming that the mutation rate is constant between loci delivers much better estimates of the parameters than single-locus estimation. This model is rather unnatural for real data, especially microsatellite data or electrophoretic data. Summing over loci using a gamma-distributed mutation rate and estimating the shape parameter 1/α increases the variation of the parameters but allows more realistic estimation of the parameters with several loci.
Simulation with unequal known parameters of 100 two-locus datasets with 25 sampled individuals for each population and 500 bp per locus
Conclusion: Our method based on coalescents with migration delivers similar or better estimates than the FST method based on the expectations derived by Maynard Smith (1970), which for certain data also produces fine results. Qualitative comparison with other methods than the FST measure, for example Rannala and Hartigan (1996), or the cladistic approach of Slatkin and Maddison (1989; see also Hudsonet al. 1992, their Table 1) shows that all these methods have a tendency to overestimate γ when Θ is not very high and all methods have quite large variances. The maximum-likelihood estimators have smaller variances than the other approaches. But the biggest drawback of the FST methods compared to ML methods seems to be their inability to estimate population size and migration rate for data where the homozygosity between populations happens to be equal to or less than the homozygosity within a population (see also Hudsonet al. 1992).
Direct comparison of the different methods is currently difficult because each one estimates a different number of parameters or uses a different migration model. We expand our model to accept more populations and different migration schemes. This should then facilitate direct comparison. We believe that taking the history of mutation into account and summing over all possible genealogies gains the most information from the data. We believe that our approach and the approaches of Nath and Griffiths (1996) and Bahlo and Griffiths (http://www.maths.monash.edu.au/~mbahlo/mpg/gtree.html/) will produce better estimates than other approaches.
—Log-likelihood curves of shape parameter 1/α. The true 1/α is 1.0. The curves were evaluated with Θ and γ parameters from the maximum-likelihood estimate. Left, 10 runs with 10 loci each. Right, sum of the log-likelihood curves over 10 runs with 10 loci each.
The methods described here are available in the program Migrate over the World Wide Web at http://evolution.genetics.washington.edu/lamarc.html/. We distribute it as C source code and also as binaries for several platforms including PowerMacintoshes and Windows 95/98 for Intel-compatible processors.
APPENDIX
Derivation of the importance sampling function: We want to calculate Equation 1 but would need to sum over all possible genealogies. This function can be transformed into an importance sampling function by assuming that
Data models: Our migration estimation divides naturally into two parts: calculation of
For sequences we are using the model of change originated by one of us (J.F.) as implemented in PHYLIP 3.2 in 1984, described in Kishino and Hasegawa (1989) and described as F84 by Swofford et al. (1996). It is a variant of the Kimura two-parameter model, which allows for different transition and transversion ratios and variable base frequencies. In migrate the model is the same as in PHYLIP 3.5 and includes also rate variation among sites (Felsenstein and Churchill 1996), but the inclusion of rate variation between sites was not tested.
For microsatellite data, we have implemented a one-step mutation model in which the probability of making a net change of i steps in time interval u is
Our “infinite” allele model is approximated with a (k + 1)-allele model. The observed alleles are
FST calculations: We use a calculation based on
FST-based parameter estimation
Acknowledgments
We thank Mary K. Kuhner and Jon Yamato for many discussions, much help in debugging, and comments on the manuscript. This work was funded by National Science Foundation grant no. BIR 9527687 and National Institutes of Health grant no. GM51929, both to J.F., and in part by a Swiss National Funds fellowship to P.B.
Footnotes
-
Communicating editor: S. Tavaré
- Received January 5, 1998.
- Accepted February 23, 1999.
- Copyright © 1999 by the Genetics Society of America
