Estimating Effective Population Size and Migration Rates From Genetic Samples Over Space and Time

Corresponding author: Jinliang Wang, Zoological Society of London, Regent's Park, London NW1 4RY, United Kingdom., jinliang.wang{at}ioz.ac.uk (E-mail)

Communicating editor: J. B. WALSH

	ABSTRACT

TOP ABSTRACT DEFINITIONS AND ASSUMPTIONS IMMIGRATION FROM AN... TWO FINITE DEMES SIMULATION RESULTS APPLICATION TO A REAL... DISCUSSION APPENDIX LITERATURE CITED

In the past, moment and likelihood methods have been developed to estimate the effective population size (N_e) on the basis of the observed changes of marker allele frequencies over time, and these have been applied to a large variety of species and populations. Such methods invariably make the critical assumption of a single isolated population receiving no immigrants over the study interval. For most populations in the real world, however, migration is not negligible and can substantially bias estimates of N_e if it is not accounted for. Here we extend previous moment and maximum-likelihood methods to allow the joint estimation of N_e and migration rate (m) using genetic samples over space and time. It is shown that, compared to genetic drift acting alone, migration results in changes in allele frequency that are greater in the short term and smaller in the long term, leading to under- and overestimation of N_e, respectively, if it is ignored. Extensive simulations are run to evaluate the newly developed moment and likelihood methods, which yield generally satisfactory estimates of both N_e and m for populations with widely different effective sizes and migration rates and patterns, given a reasonably large sample size and number of markers.

ONE of the major applications of population genetics theory has been in the estimation of genetic and demographic properties of populations. Quantities such as the mutation rate, the migration rate, and the recombination rate are all in principle estimable from population genetic data, at least as products with the effective population size.

The variance effective population size (N_e) is an important quantity in evolutionary biology, as a description of the rate at which genetic variance changes due to genetic drift. A variety of methods have been developed to estimate N_e (SCHWARTZ et al. 1999 ), including predictive equations based on life history data (reviewed in CABALLERO 1994 ; WANG and CABALLERO 1999 ), the lethal allelism approach (DOBZHANSKY and WRIGHT 1941 ), and a linkage disequilibrium approach (HILL 1981 ). Several authors have developed methods to estimate N_e for a population from temporal genetic data, that is, data on the frequencies of marker alleles taken at two or more time points from the same populations. In principle, if other population genetic processes such as mutation, migration, and selection are negligible in their effects on allele frequency change over the time frame of the study, then the observed changes in allele frequency can be assumed to be due to the effects of genetic drift. Using what is known about the effects of genetic drift, then, we could estimate N_e from such temporal data. This so-called "temporal approach" was pioneered by KRIMBAS and TSAKAS 1971 , with subsequent substantial statistical refinements by others (e.g., PAMILO and VARVIO-AHO 1980 ; NEI and TAJIMA 1981 ; POLLAK 1983 ; TAJIMA and NEI 1984 ; WAPLES 1989 ; WILLIAMSON and SLATKIN 1999 ; ANDERSON et al. 2000 ; WANG 2001A ). These methods have been applied to a large variety of species and populations (e.g., FUNK et al. 1999 ; LABATE et al. 1999 ; FIUMERA et al. 1999 , FIUMERA et al. 2000 ; KANTANEN et al. 1999 ; TURNER et al. 1999 ; see also many earlier references cited in WILLIAMSON and SLATKIN 1999 ).

A key assumption made by all these approaches, however, is that the only factor involved in allele frequency change is genetic drift. All systematic forces (selection, mutation, and migration) are assumed to be absent. In our view, it is often legitimate to ignore the effects of mutation on the timescale involved in most empirical studies. It also may be reasonable to neglect the effects of selection because direct selection on most markers may be unlikely to be strong enough to cause substantial changes in their frequencies. (The effects of indirect selection acting on other linked loci on the allele frequency change of the marker locus may arguably be included as a "drift" effect and legitimately included in a N_e term.) However, for most populations the effects of migration are not negligible. A major weakness of these temporal methods is that they ignore the effects of migration, which can substantially bias estimates of N_e, either upward or downward (see below). The purpose of this article is to include the effects of migration in the estimation of demographic parameters from temporal genetic data.

Migration can affect the estimation of N_e in two ways. In the short term, migration can change allele frequencies quite quickly. Given that the temporal methods estimate N_e by measuring the pace at which allele frequencies change, migration causes a population to behave as if strong drift was changing allele frequencies very rapidly, resulting in an estimate of N_e that is too low. In the longer term, the problem is reversed. Constant migration and drift would cause the populations to approach an equilibrium level of genetic differentiation and the rate of change of allele frequencies in a deme to slow down to approach that predicted by the effective size of the whole metapopulation. Thus the effective population size of the deme would be substantially overestimated in the long term. These effects are discussed later in this article in more mathematical detail.

This problem can be avoided by simultaneously estimating N_e and the migration rate, which is the approach taken in this article. It turns out that while both migration and drift affect allele frequency changes in a population, they do so differently in terms of the direction and magnitude of these changes. Under migration, the allele frequency of a focal population would tend to become increasingly similar to that of the source population. Under drift, however, the allele frequency of a focal population changes purely at random, irrespective of that of the source population. Furthermore, the magnitude of allele frequency change of a focal population depends on the allele frequency difference between the focal and source populations under migration. Migration has no effect on alleles that happen to be at the same frequencies both in the deme in question and in the populations providing those migrants, but migration is expected to change allele frequencies substantially in cases when the source and recipient demes have very different allele frequencies. Genetic drift, on the other hand, is not affected by differences in allele frequencies among demes, but depends only on the local allele frequency. These differences in the patterns of changes of allele frequencies due to migration and drift allow us to estimate them separately, so long as enough data are from multiple loci. Fortunately, these kinds of data are often available with modern biochemical methods.

The estimation of migration rate and migrant number has a long history. The oldest and most widely used genetic method assumes that populations are at equilibrium in a system approximated by WRIGHT's (1931) island model and then uses estimates of Wright's F_ST to calculate 1/(4N_e m + 1), where m is the local immigration rate (see SLATKIN 1987 ). This method has been roundly criticized (WHITLOCK and MCCAULEY 1999 ), in particular because of the unrealistic assumptions it requires about the demographic structure of the species. Other approaches include genetic stock identification (MILNER et al. 1985 ), assignment methods (RANNALA and MOUNTAIN 1997 ), two-locus descent measures (VITALIS and COUVET 2001 ), coalescent methods (NIELSEN and WAKELEY 2001 ), isolation-by-distance methods (see review by ROUSSET 2001 ), and migration matrix-likelihood models (BEERLI and FELSENSTEIN 2001 ). The first two approaches require detailed knowledge about the possible sources of migrants and do not allow for genetic drift. The latter four approaches assume that the populations are at an equilibrium between the effects of genetic drift and migration. There is a need for a method that does not make restrictive assumptions about genetic and demographic equilibria, yet allows simultaneous estimation of the effects of drift and migration.

Here we extend previous temporal methods for estimating the N_e of a single isolated population so that they can be used to estimate N_e and m jointly for a metapopulation. We first develop an expectation based on a moment estimation approach, which allows an approximate description of the behavior of the estimator and the system. We follow this by developing a likelihood method based on that of WILLIAMSON and SLATKIN 1999 , ANDERSON et al. 2000 , and WANG 2001A . Extensive simulations are run to check the performance of the moment and likelihood estimators in terms of precision and accuracy for the estimation of N_e and m. Two metapopulation models (a focal small population receiving immigrants from an infinitely large source population and two small populations connected by migration) are considered, and the applications of the methods to more complicated metapopulation models (e.g., a small number of demes in island or stepping-stone migration models) were checked by simulations. As an example for the application of our methods, a spatial and temporal genetic data set on Triturus newts (JEHLE et al. 2001 ) was reanalyzed by our methods for the simultaneous estimation of migration and drift.

	DEFINITIONS AND ASSUMPTIONS

TOP ABSTRACT DEFINITIONS AND ASSUMPTIONS IMMIGRATION FROM AN... TWO FINITE DEMES SIMULATION RESULTS APPLICATION TO A REAL... DISCUSSION APPENDIX LITERATURE CITED

We use temporal data to estimate the effective size (N_e) and immigration rate (m) of a partially isolated population. For brevity, we refer to the population for which parameters are being estimated as the focal population and that from which the migrants come as the source population. Data and parameters for the source population are given the subscript A, and those for the focal population, B. In the first section, we assume that the source population is infinitely large and therefore not changing genetically through the time frame of the study. In the section after that, we assume that both populations are finite and jointly estimate the properties of both. In the latter case, source and focal populations are obviously reciprocal. Throughout we assume diploidy. For data from haploid species or loci, the N_e estimated by these models would be double the actual effective size.

For the infinite-source population model, a minimum of three samples is required to estimate N_e and m of the focal population: two temporally spaced ones from the focal population and one taken at any time from the source. For the two small population models, a minimum of four samples is necessary, two temporally spaced ones taken concurrently from each population. Let us denote the size of these samples S, with subscripts A or B to describe whether the sample is from the source or focal population, respectively, and a subscript t to describe the generation in which the sample is taken. Thus S_B,0 is the number of individuals sampled from the focal population in the initial generation. We use x to refer to the relative frequency of a particular allele in a given sample. Therefore x is an estimator of the true frequency (p) of that allele at that time in that population. For convenience, we define q = 1 - p. These quantities, x, p, and q, are subscripted in the same way as S.

Several summary statistics of these quantities turn out to be useful. We define p_AB p_A,0 - p_B,0 to be the initial difference between the allele frequencies of the two populations and p_B,t p_B,t - p_B,0 to be the change in allele frequency by generation t in deme B. In the derivations of the moment-based estimators, we also use _B,i to refer to the change in generation i due to genetic drift of the allele frequency in population B. Each of the quantities can refer to population A if the subscript B is replaced by A.

We make several assumptions throughout. These are that the alleles under study are not under direct selection, that mutation is a negligible source of allele frequency change, that potential sources of migrants to the focal populations can be identified a priori, that all samples are randomly taken from their populations, and that the sampling event does not change the available pool of reproductive individuals. The latter can be realized by sampling after reproduction, sampling with replacement, or sampling from a population much larger than the sample. This sampling scheme therefore corresponds to WAPLES's (1989) plan 1.

Table 1 summarizes most of the terms used in this article.

	IMMIGRATION FROM AN INFINITE-SOURCE POPULATION

TOP ABSTRACT DEFINITIONS AND ASSUMPTIONS IMMIGRATION FROM AN... TWO FINITE DEMES SIMULATION RESULTS APPLICATION TO A REAL... DISCUSSION APPENDIX LITERATURE CITED

Theoretical background:
Imagine that the focal population is receiving migrants at rate m per generation from an infinite, unchanging source population. For effectively neutral alleles, changes in allele frequency over a short period of time are a function of the local effective population size and migration rate. Mutations are negligible compared to drift and migration because of their rare occurrence and the short time span. Under these conditions, and assuming that migration precedes population regulation, the allele frequency at generation i is given by

(1)

Note then that

(2)

Expanding across generations, we can find

(3)

so

(4)

Under the assumptions of Wright-Fisher sampling, the terms are independently binomially distributed. The expected value of all of the terms is zero, because drift does not change on average the allele frequency. The variance of each of these terms is given by p_B,i(1 - p_B,i)/N_e for haploids or p_B,i(1 - p_B,i)/2N_e for diploids. In reality, we do not know p_B,i for any generation between samples, and we can only estimate p_B,i for generations in which there was a sample. This uncertainty is the root of error in the estimation of statistics derived from these equations. In this article, we first estimate this error by using the approximate moments of the distribution of allele frequencies, and then we describe a likelihood approach to the estimation of N_e and m.

Moment approximations:
N_e and m affect the first and second moments of the change in allele frequency over time. Under certain circumstances, we can use these moments to estimate m and N_e. We do this in a two-step process, first estimating m from the mean change in allele frequencies at all loci and then using this value in the estimation of N_e from the second moment of allele frequency change.

Estimating the migration rate: The change in allele frequency is given by Equation 4 above. Taking the expectation of p_B,t we find

(5)

because the expected values of all the terms are zero. Assuming that t is known without error, we can therefore make a one-locus estimate of (1 - (1 - m)^t) from p_B,t/p_AB.

To predict the optimal sampling strategy and to approximate standard errors, we estimate the error variance of an estimate of m. The expected error in estimating a ratio is an approximate function of the error variance of estimating the numerator and the denominator and the error covariance between the two (see LYNCH and WALSH 1998 , Appendix 1). The error in estimating (1 - (1 - m)^t) can be approximated from the error in estimating p_B,t/p_AB. The error variances (V_err[ ]) and error covariance (COV_err) are given by

(6)

Note that the 2's in the denominators drop out when these equations are applied to haploids. The approximation in the first equation here comes from using the initial value of pq to predict the amount of variance in final allele frequency due to genetic drift.

Using these last equations and Equation A1.19b from LYNCH and WALSH 1998 , we can find the expected error variance for the estimate of (1 - (1 - m))^t:

(7)

If the three sample sizes are roughly equivalent, then the error in estimating m is influenced more by the allele frequency in population B than by that in the source population, especially if the product mt is small. The terms involving only the sample sizes, m, and t are presumably constant across loci, so the error variance is largely proportional to 1/(p_AB)². The minimum variance possible when combining L loci can be obtained by weighting each locus by the inverse of its expected error variance. The best estimator of m is given approximately by

(8)

where the weight for locus l is w_l = (x_AB)², F_l = x_B,t/x_AB, and W is the sum of w_l across loci.

Estimating N_e: Under the assumption that the total change in allele frequency per generation is not large (i.e., N_e is not very small and m is not very large), each of the terms in the equation describing allele frequency change is roughly independent of the others. The variance in allele frequency change at generation i due to genetic drift within the focal population is V[_B,i] = p_B,iq_B,i/2N_e,B. If the allele frequency is not changing rapidly, then we can assume that the value of pq is approximately the same each generation. A reasonable estimator of the mean value of pq is (x_B,0(1 - x_B,0) + x_B,t (1 - x_B,t))/2. This is implicitly assuming that the allele frequency on average takes a linear path from the initial to the final value, but in any case is an ad hoc approximation. Then using Equation 4, we can find

(9)

where M = (1 - (1 - m)^2t)/((2 - m)m). Again, the coefficient of N_e,B of 2 drops out when this formula is used for haploids. From Equation 9 and using the estimate of m, we obtain an estimator of 1/(2N_e,B),

(10)

where = (1 - (1 - )^2t)/((2 - )).

In general we do not have exact information about the change in population allele frequencies, but only an estimate based on sample data. Because of sampling variance, E[x_B,t²] is biased upward as an estimator of p_B,t²:

(11)

Thus we can estimate

(12)

The error variance per locus for estimating 1/(2N_e) is approximately proportional to 1/ (calculations not shown), so a minimum-variance estimator of 1/(2N_e) can be obtained by weighting alleles by approximately, estimated by (x_B,0(1 - x_B,0) + x_B,t(1 - x_B,t))/2.

Choosing the right value of t: With a large separation in time between the first and second samples, estimates of m become problematic. The estimate of (1 - (1 - m)^t) becomes more and more similar for different values of m as t gets larger (lim_t[(1 - (1 - m)^t)/m] = 0), yet the error variance in estimating this quantity increases linearly with t. Furthermore, as t increases, the allele frequency change approaches a distribution determined largely by the relative sizes of the effective size and m, rather than by the absolute values of either. When t is large enough that the populations reach an equilibrium between genetic drift and migration, then 0, irrespective of the actual true migration rate. This is clear by inspecting the single-locus estimator of m, = 1 - ^t. As t , x_B,t no longer increases and x_B,t/x_AB is expected to reach a constant smaller than one with time, therefore resulting in 0.

If the time interval is too short and m too small, then there is less power to measure m because too few migration events may have occurred. Larger samples with more informative markers are therefore necessary to obtain a good estimate of m when (1 - m)^t is close to one or zero.

Similarly, the estimation of N_e is also biased if the sampling interval is too long. As t , we have 0 shown above and Equation 10 reduces to Est(1/2N_e,B) = E[p²_B,t]/(t). If the source population is infinite in size, then E[p²_B,t] tends to a constant as t , and therefore _e,B . If the source population has a finite size, then E[p²_B,t] continues to increase as t (so long as the alleles are not fixed in the population), but at a slower rate corresponding to the total size of the focal and source populations. In both cases, the estimate of N_e,B approaches the effective size of the whole species rather than that of the focal population. Again, there is less power for estimating N_e,B if t and S/(N_e,B) are too small. With a small quantity of tS/(N_e,B), the change in allele frequency due to drift would be overwhelmed by that from sampling.

As a numerical example, Fig 1A shows the average estimates of 1/(2N_e,B) and m over 10,000 replicates plotted against sampling interval t, and Fig 1B shows the coefficients of variation (CV) of the estimates when t = 1–12 where the estimation is unbiased. (Methods for these simulations are given in the Appendix.) The true values being estimated are 1/(2N_e,B) = 0.005 (N_e,B = 100) and m = 0.2 for the focal population, which is assumed to be in a metapopulation containing 10 other populations of the same size (N_e = 100) with an island migration model. The effective size of the source population is therefore 1000, sufficiently large although not infinite. As can be seen, estimates of both 1/(2N_e,B) and m are close to the true values when t is small, but approach the true values for the whole species, 0.00045 (N_e = 1112, calculated using the known parameters) and 0, with increasing t. The maximum number of generations allowed for valid joint estimates of 1/(2N_e,B) and m and the number of generations required to reach the equilibrium so that the species N_e is estimated depend on m, the effective sizes of the focal and source populations, and the initial genetic differentiation when the first sample was taken (see below).

View larger version (24K): In this window In a new window Download PPT slide   Figure 1. Means (A) and coefficients of variation (B) of estimates of 1/(2Ne) and m as a function of sampling interval t (in generations). The population contains 11 subpopulations of an equal Ne (100), migrating in the island model with rate 0.2. After the population reaches quasi-equilibrium, two samples separated by t generations are taken from each of the 11 subpopulations. The size of each sample from the focal subpopulation is 200 individuals, and the sample from any source subpopulation is 10 individuals. The 20 samples from the source subpopulations are combined into a single sample. For each sample, 20 loci, each having four alleles, are genotyped and used in the estimation. In the top graph, the averages over 10,000 replicate estimates of m and 1/(2Ne) from the moment estimator are denoted by thin and thick solid lines, respectively. The true value for m is denoted by the horizontal thin dashed line, and the true values of 1/(2Ne) for the focal subpopulation (0.005) and the whole population (1/2224) are denoted by the top and bottom horizontal thick dashed lines. Both axes of this graph are in log scale. In the bottom graph, the coefficients of variation over 10,000 estimates of m and 1/(2Ne) are denoted by thin and thick lines, respectively, from the moment estimator.

For this specific numerical example, the estimation is unbiased when t < 12. The precision of the estimators changes considerably over this short sampling period (Fig 1B). As expected, the best estimates are obtained at an intermediate value of t (t = 2–4 for the example). The optimal sampling interval depends, again, on m, the effective size of the focal population and the initial genetic differentiation when the first sample was taken. With decreasing m and/or increasing initial differentiation, for example, the optimal t is expected to increase.

In experimental design, it is difficult to determine beforehand the optimal interval between sampling for maximum power of the temporal method, because the power depends partially on the parameter being estimated (see e.g., Equation 7). Without any prior knowledge of the demographic history of the populations in question, it is reasonable to consider a short rather than a long sampling interval. While a too long sampling interval results in biased estimates of N_e and m, a too short sampling interval only decreases the precision of the estimates. With powerful statistical methods (e.g., maximum likelihood shown below) and large sample size and marker information, the estimation precision can be improved greatly even if t is as short as one generation.

The consequences of assuming that m = 0: Other methods using temporal data for estimating the effective size have assumed that the immigration rate is zero. We can use the equations above to investigate the effects of making this assumption when it is in fact false.

When m 0, the expected change in allele frequency is zero, M = t, and Equation 12 reduces to

(13)

which is similar to NEI and TAJIMA's (1981) estimator using F_c.

When m > 0, the variance in allele frequency change over the interval between samples is a function of both genetic drift and migration. This can be either greater or smaller than that expected by drift alone, depending on, among other parameters, t. In the short term, the change in allele frequency due to the joint effects of migration and drift will be larger than that expected by drift alone. If t = 1, for example, we have E[(p_B,t - p_B,0)²] = p_AB²m² + p_B,0(1 - p_B,0)/2N_e,B. In this case, therefore, the squared change in allele frequency is increased relative to that due to drift only, which is p_B,0(1 - p_B,0)/2N_e,B. Falsely assuming that m = 0 would on average cause an overestimate of the amount of drift and an underestimate of N_e,B. In contrast, as the time between samples gets large, the change in allele frequency is constrained by migration to a term only somewhat larger than the squared initial difference in allele frequency between the source and focal populations: lim_tE[(p_B,t - p_B,0)²] = p_AB² + p_B,0(1 - p_B,0)/(N_e,Bm(2 - m)). With an increasing sampling interval (t), therefore, a false assumption of m = 0 would cause a downward bias in the estimation of the amount of genetic drift and an upward bias in estimating N_e. As the number of generations between samples continues to get larger, however, the estimate would approach the effective size of the whole species, rather than just the effective size of the focal population, no matter whether the assumption of m = 0 is made or not. These biases can be indefinitely large.

Fig 2 shows the biases of N_e estimation when a false assumption of m = 0 is made. The true parameter values are N_e = 100 and m = 0.1, and the populations are assumed to be either completely differentiated (no previous migration, case A) or at drift-migration equilibrium (case B, where F_st 0.02) when the initial sample is taken. Ignoring migration results in an overestimation of 1/(2N_e) when sampling interval (t) is small and underestimation when t is large. Compared with case A, the initial overestimation of 1/(2N_e) in case B from ignoring migration is much smaller, because the populations are more similar genetically (p_AB small) and therefore migrants have less effect on allele frequency changes. With a larger m and smaller N_e, the initial overestimation of 1/(2N_e) can be dramatic even when the populations are initially at equilibrium (data not shown). For both cases, the estimation is reasonably good for a large range of sampling interval when immigration is accounted for. However, the sampling interval cannot be too long. Otherwise, 1/(2N_e) is increasingly underestimated with t, approaching the true value for the whole species.

View larger version (15K): In this window In a new window Download PPT slide   Figure 2. Bias of 1/(2Ne) estimates as a function of sampling interval (in generations). The true value being estimated is 1/(2Ne) = 0.005 (Ne = 100) denoted by the horizontal dashed line, and the true migration rate is 0.1. The populations are assumed to have no previous migration (case A) or be at drift-migration equilibrium (case B) when the first sample from the focal population is taken. For cases A and B, the averages over 1000 replicate estimates of 1/(2Ne) are denoted by thin and thick solid lines, respectively, from the moment estimator developed herein that takes immigration into account, by thin and thick dashed lines, respectively, from the moment estimator ignoring immigration (Equation 13). In the estimation, we assumed 20 loci, each having four alleles, and two samples from the focal and one sample from the source population, each sample having 200 individuals. Note that both axes are in log scale.

Likelihood model:
Let us first consider a single locus with two alleles, and we focus on a particular allele. For brevity in this subsection, subscript B is dropped for the focal population. To estimate m and N_e of the focal population jointly, at least two temporally spaced samples from the focal and a single sample from the source population (A) are necessary. Because the three sampling events are independent, the probability of getting a set of samples with allele frequencies denoted by x is

(14)

Each of these samples can be assumed to follow a binomial distribution (see WAPLES 1989 for a discussion of when this is appropriate); therefore

where S is the number of diploid individuals sampled. For haploid species, the coefficient 2 of S should be dropped.

The likelihood (L) of a particular set of values of N_e and m is the probability of observing the data given those values. We have, therefore,

(15)

As we can see, the key problem for the likelihood function is to calculate Pr[p₀, p_t, p_A|N_e, m]. By the definition of conditional probability, we have

(16)

Determining the second term in Equation 16, Pr[p₀, p_A|N_e, m], is potentially quite tricky. If we assume, however, that we do not know anything about how close the focal population is to a migration-drift equilibrium nor about the history of these populations, then this can be given reasonably as a product of the uniform distributions of p₀ and p_A, Pr[p₀, p_A|N_e, m] = Pr[p₀|N_e] x Pr[p_A]. As WILLIAMSON and SLATKIN 1999 point out, the range of the possible values for the uniform prior, Pr[p₀|N_e], will depend on whether only alleles segregating in the initial population are used. There are 2N_e + 1 possible configurations for a biallelic locus in a population with size N_e, of which 2N_e - 1 are polymorphic. With the uniform prior, therefore, any possible configuration for the starting frequency p₀ has an equal probability of either 1/(2N_e + 1) or 1/(2N_e - 1) if we consider all or only polymorphic configurations, respectively. For the source population, we use a discrete approximation, Pr[p_A,j = j/J] = 1/(J + 1) for j = 0, 1, 2, ... , J. The estimation does not change essentially once J is reasonably large (say, 1000), as expected.

The first term in Equation 16, Pr[p_t| p₀, p_A, N_e, m], can be calculated by using the transition matrix T (EWENS 1979 ). The probability of a transition from the jth configuration of the parental population at time t - 1 to the ith configuration of the offspring population at time t is

where p^*_j = (j/2N_e)(1 - m) + mp_A and (i, j = 0, ... , 2N_e). This accounts for the deterministic change in allele frequency due to migration as well as the genetic sampling due to drift. Given the initial distribution of allele frequency of the focal population P₀ [a column vector of 2N_e + 1 elements with values 1/(2N_e + 1) or 1/(2N_e - 1) as assumed above] and values of p_A, N_e, and m, the distribution after t generations, P_t, is calculated by the recurrence equation

(17)

These are all the parts necessary to calculate the likelihood, above. For multiple independent, biallelic loci, the joint likelihood of N_e and m is calculated as the product of the likelihoods for each locus.

For loci with more than two alleles, the calculation of the likelihood function is problematic. The number of possible configurations, C = (2N_e + k - 1)!/(2N_e! (k - 1)!) for a diploid locus with k alleles (FELLER 1950 ; WILLIAMSON and SLATKIN 1999 ), increases rapidly with k, prohibiting the evaluation of the exact likelihood due to the limitations of computer memory and processing speed. For a moderate value of N_e, even k = 3 poses a serious computational problem for the likelihood method. ANDERSON et al. 2000 proposed a Monte Carlo method to evaluate the likelihood function with data on multiallelic markers. Their method does not need much storage but is highly demanding computationally and may not converge when k is large (say, k = 10) or some alleles are not observed in all samples. To simplify the likelihood computation for multiallelic loci, we follow our previous approach (WANG 2001A ), which transforms a k-allele locus into k biallelic "loci," each having one of the k alleles with all the other alleles pooled. The overall log-likelihood is approximated by the sum of the log-likelihoods across loci multiplied by the factor of (k - 1)/k to account for the dependence among such converted loci.

The likelihood method can also use more than two samples from the focal population simultaneously to obtain a mean N_e and an average m during the whole sampling period. If samples are taken from the focal population at generations t0, t1, t2, ... , tg, then the joint likelihood of N_e and m is

(18)

Similarly, a continuous growth model, such as exponential growth, could be fitted to the temporal data and joint estimates of initial effective size, growth rate, and immigration rate can be obtained (WANG 2001A ). A similar model to account for the changes of immigration rate over the sampling period could be developed using multiple samples.

	TWO FINITE DEMES

TOP ABSTRACT DEFINITIONS AND ASSUMPTIONS IMMIGRATION FROM AN... TWO FINITE DEMES SIMULATION RESULTS APPLICATION TO A REAL... DISCUSSION APPENDIX LITERATURE CITED

Up to this point, we have considered the demographic properties of a single deme that receives immigrants from an effectively infinite source, that is, a source so large that over the time course of the study it is not changing appreciably in allele frequency. While this is likely to be a reasonable approximation in many cases, there are examples in nature where only a few relatively small demes serve as the only source of migrants to each other. This case is more difficult, both analytically and computationally, but here we consider the case of two small demes as a beginning to such an analysis.

The allele frequency at a neutral locus will change over a single generation in the two populations according to

(19)

Moment estimator:
We can find a moment estimator for the two migration rates for this two-finite-deme case, if we consider the change in the difference in allele frequency between the demes over time. We consider the simple case in which the samples are separated by only a single generation, and so we drop the subscript notation for time. Taking the expectations of Equation 19 above and remembering that the expected values of _A and _B are zero, we find

(20)

The error variance in estimating migration rates would be proportional to 1/p²_AB, so we can find a low variance estimator of the migration rates by weighting information from each allele by p²_AB.

Finding an estimator of the effective size of each deme is also possible. The variance of the change in allele frequency over time due to drift for population A is

(21)

where p^'_A is the expected allele frequency after migration, which can be estimated from Equation 19 by setting _A = 0 and using the estimated migration rates and initial allele frequencies. Thus the effective size can be estimated in this case in a manner similar to the Waples technique, except instead of using the initial allele frequency, one uses the expected initial frequency after migration. Similar derivations for N_e,B can also be obtained.

For the case of t > 1, moment estimators for both m and N_e become much more complicated (especially the latter) and are not shown here.

Likelihood estimator:
The likelihood of the effective sizes (N_e,A, N_e,B) and immigration rates (m_A, m_B) of demes A and B given the temporal data (x_A,0, x_A,t, x_B,0, x_B,t) is

(22)

In Equation 22, each of the four samples follows a binomial distribution for Pr[x|p], the initial allele frequency of each population is assumed to be in a uniform distribution, and the allele frequency distribution of each population at generation t can be calculated by using the transition matrix.

In principle, the exact likelihood could be worked out for any set of the four parameter values, but in practice this would be computationally intensive for anything but small population sizes and small sampling interval (t). The total number of configuration combinations that need to be considered in Equation 22 is ((2N_e,A + 1) (2N_e,B + 1))^t for a biallelic locus, which increases exponentially with sampling interval t. Therefore, we use the transition matrix method for the exact calculation of the likelihood for samples separated by only a single generation (t = 1) and use a Monte Carlo approach to approximate the likelihood for samples with t > 1.

Equation 22 could easily be extended to more than two samples from each deme with little increase in computational demands. It could also be extended to more than two demes, but the computational intensity increases prohibitively with the number of demes, because of the exponential increases in configuration combinations and the number of parameters for joint estimation.

	SIMULATION RESULTS

TOP ABSTRACT DEFINITIONS AND ASSUMPTIONS IMMIGRATION FROM AN... TWO FINITE DEMES SIMULATION RESULTS APPLICATION TO A REAL... DISCUSSION APPENDIX LITERATURE CITED

Infinite source:
Methods for these simulations are summarized in the Appendix. For focal populations with wide ranges of effective size and immigration rates, the maximum-likelihood and moment estimates of 1/(2N_e) and m from simulations are summarized in Table 2. The sample sizes used in the simulations are large, but not unrealistic. Indeed, a large sample size is important for the temporal methods; otherwise, the apparent changes in allele frequency would be dominated by sampling error rather than by genetic drift. For most natural populations, the actual census population sizes could be orders of magnitude larger than the effective sizes (FRANKHAM 1995 ), making both large samples and significant genetic drift possible.

The maximum-likelihood estimator gives unbiased estimates of 1/(2N_e) in all cases considered. The coefficient of variation of the estimates is roughly one-quarter, with 95% of the estimates being in the range of 50 and 150% of the true values. The moment estimator underestimates 1/(2N_e) slightly when m is large and sample size is small relative to N_e. Increasing sample size can decrease the bias dramatically. For the case of N_e = 100 and m = 0.5 in Table 2, for example, the mean ± SD of the moment estimates of 1/(2N_e) improves to 0.9 ± 0.3 (relative to the true value) when the sample size is 500 individuals. The estimates from the moment estimator are generally more dispersed than those from the maximum-likelihood estimator.

The maximum-likelihood method overestimates m slightly only when m is small in very small populations (N_e = 10). This is expected because maximum-likelihood estimates have a lower bound at zero. Increasing sample size or marker information (number of loci and alleles) in a sample decreases the bias and increases precision. In contrast to 1/(2N_e), both the moment and the maximum-likelihood methods give similar estimates of m. However, some moment estimates of m are negative when m is small.

In general, both maximum-likelihood and moment methods give reasonable estimates of N_e and m. The accuracy and precision of the two estimators are similar for m, but the likelihood method is better for N_e. This is partially because moment methods must use the estimate of m in estimating N_e. Such a two-step estimation procedure is expected to be inferior to the joint estimation of m and N_e from the likelihood approach, because the estimation error of m is not accounted for in estimating N_e. Also, when m is large and therefore the changes of allele frequencies are rapid, the approximation to the mean value of pq in deriving the moment estimator of N_e may be no longer good enough.

The distribution of the estimates of N_e [or 1/(2N_e)] generally has a long tail, with a significant proportion of the estimates being much larger (smaller) than the true value. This is similar to the observation made on the case of a single isolated population (NEI and TAJIMA 1981 ; WAPLES 1989 ; WANG 2001A ). For the particular parameter set of N_e = 100, m = 0.1, the distributions of N_e and m estimates over 10,000 replicates are shown in Fig 3. Again, there is very little difference in the distribution of estimates of m between the two techniques, the mean ± SD being 0.114 ± 0.048 for maximum likelihood and 0.115 ± 0.049 for the moment estimator. For N_e estimation, however, maximum likelihood is much more accurate and precise than the moment approach. The mean ± SD of estimated 1/(2N_e) is 0.0049 ± 0.0020 for maximum likelihood and 0.0041 ± 0.0025 for the moment estimator. The latter gives a significant proportion of extremely large N_e estimates.

View larger version (19K): In this window In a new window Download PPT slide   Figure 3. The frequency distributions of Ne (A) and m (B) estimates from maximum-likelihood (thick lines) and moment (thin lines) estimators over 10,000 replicates. Two samples separated by four generations are taken from the focal population, and one sample is taken from the infinite-source population. Each sample contains 100 diploid individuals, which are genotyped for 10 loci, each having four alleles. The true parameter values are Ne = 100, m = 0.1, indicated by the dashed vertical lines.

Typically, the likelihood surface is smooth and has a single maximum. This means that a single arbitrary starting point is sufficient to obtain maximum-likelihood estimates using Powell's quadratic convergence method (PRESS et al. 1996 ). Of course the starting point affects the time spent in the search process, and the values provided by the moment estimates usually act as a good starting point close to the point with the maximum likelihood. For a particular replicate with N_e = 100 and m = 0.1 in Table 2, the likelihood curve is shown in Fig 4. The maximum-likelihood estimates of N_e and m are 102 and 0.107, respectively, with 95% confidence intervals being 65–172 and 0.038–0.191, respectively.

View larger version (68K): In this window In a new window Download PPT slide   Figure 4. Log likelihood as a function of estimates of Ne and m for a particular simulated data set with true parameter values Ne = 100 and m = 0.1. A single sample from the large source and two samples separated by two generations from the focal population were taken with replacement, the sizes of the three samples being 200 diploid individuals. Twenty loci, each having four alleles, were used in the estimation.

Multiple samples over time can increase the estimation precision for the average 1/(2N_e) and m during the entire sampling period. Table 3 summarizes the estimates of 1/(2N_e) and m, relative to the true parameter values of 1/(2N_e) = 0.01 and m = 0.2, from maximum-likelihood and moment estimators using two to six temporal samples from the focal population. As expected, increasing the number of samples has no influence on the accuracy, but increases the precision dramatically, as indicated by the decreasing SD and the width between the 2.5 and 97.5 percentiles. The likelihood method gives unbiased estimation of 1/(2N_e) and m, while the moment estimator results in underestimation of 1/(2N_e) and overestimation of m. The poor accuracy of the moment estimator is partially due to the highly polymorphic loci (10 alleles per locus) and small sample sizes, which lead to a high proportion of rare alleles in the samples (WANG 2001A ).

Quasi-infinite source, the island model and the stepping-stone model:
In reality, the source population is not infinite in size, although it could be very large. Moreover, there can be different numbers of source populations for a focal population with migration in different patterns. As long as we can identify the source populations a priori, however, the parameters for the focal population can be estimated approximately using the methods developed for the infinite-source model by treating the source populations collectively as the single "infinite" source. Here we explore numerical examples to show that the infinite-source model is very robust and works well for small source populations with various models of migration.

Consider a metapopulation consisting of 21 populations of equal effective size N_e. Each population receives, at each generation, a proportion of m migrants taken equally from each of the remaining 20 populations in the island model or from each of the two neighboring populations in the circular stepping-stone model. For the estimation of m and N_e of a focal population, we take two samples separated by two generations from the focal and each of four source populations (the island model) or each of the two neighboring populations (the stepping-stone model). Each sample had 50 individuals and two alleles per locus for N_e = 10 and 100 individuals and eight alleles per locus for N_e = 100. The estimates of 1/(2N_e) and m relative to their true values are summarized in Table 4 for maximum-likelihood and moment estimators. In all cases considered, maximum likelihood provided reasonable estimates of both 1/(2N_e) and m, while the moment estimator gives satisfactory estimation of 1/(2N_e) and m only when drift is not very weak relative to migration. Otherwise, 1/(2N_e) is underestimated, similar to the situation of an infinite large source population shown in Table 2.

We also tested the infinite-source model with the extreme case of two small populations of equal N_e and m, shown in Table 5. Although the source population (A) is small, the estimation of both 1/(2N_e) and m of the focal population (B) is still not bad, at least within a factor of two. The moment estimator gives similar results (not shown). Like the results shown above, the worst case occurs with an extremely high migration rate, where populations can quickly reach the drift-migration equilibrium. In this case, increasing sample size improves the estimation considerably. For the situation of N_e = 100 and m = 0.5 in Table 5, for example, the mean (relative) estimate ± SD is 0.969 ± 0.204 for 1/(2N_e) and 0.809 ± 0.204 for m when the sample size is increased to 400 individuals.

Two-deme models:
As a numerical example for the case of two demes with sampling interval t = 1, Fig 5 shows the frequency distributions of the estimates of N_e,A, m_A, N_e,B, and m_B. As can be seen, both estimators give reasonably good estimates of effective size and migration rate, the precision of the likelihood estimator being always higher than that of the moment estimator.

View larger version (21K): In this window In a new window Download PPT slide   Figure 5. The frequency distributions of Ne,A (A), mA (B), Ne,B (C), and mB (D) estimates from maximum-likelihood (thick lines) and moment (thin lines) estimators over 10,000 replicates, using the two-finite-deme model. Two samples separated by one generation are taken from each of the two populations, each sample containing 200 diploid individuals. Ten loci, each having four alleles, are used in the estimation of the true parameters Ne,A = 10, mA = 0.05, Ne,B = 100, and mB = 0.2, indicated by dashed vertical lines. The long and short vertical lines below the horizontal axis indicate the mean and 2.5 and 97.5 percentiles for the maximum-likelihood (thick lines) and moment (thin lines) estimators. Note that the 97.5 percentiles for Ne,B are >1000 and not plotted in the graph.

For the case of sampling interval t > 1, the moment estimator for N_e becomes quite complicated and the likelihood estimator becomes computationally intensive. In principle, however, both methods are not limited by the sampling interval, except for the case of very large t where populations are approaching an equilibrium between drift and migration (Fig 1 and Fig 2). At intermediate values of t, both methods should actually be more powerful than at t = 1 (Fig 1B). Because of the computational demands of the Monte Carlo approach, we have run only a few replicates for t > 1 for the likelihood method, yielding reasonable estimates of both effective size and migration rate, as expected (data not shown).

	APPLICATION TO A REAL DATA SET

TOP ABSTRACT DEFINITIONS AND ASSUMPTIONS IMMIGRATION FROM AN... TWO FINITE DEMES SIMULATION RESULTS APPLICATION TO A REAL... DISCUSSION APPENDIX LITERATURE CITED

JEHLE et al. 2001 used microsatellite markers to investigate the spatial and temporal genetic differentiation of syntopic newts (Triturus cristatus, T. marmoratus). Individuals from three ponds in western France, situated 4–10 kilometers apart and inhabited by both species, were sampled and genotyped for eight microsatellite loci. Here we focus on T. marmoratus and consider the subpopulations in ponds 1 and 2 only for which temporal samples are available. The subpopulation in pond 3 is ignored because it is much farther away from the other two and has only a single sample. The first sample was collected as sacrificed larvae in 1989 (pond 1) or 1986 (pond 2), and the second was taken nondestructively from both larvae and adults in both ponds in 1998. The allele counts for five polymorphic loci were converted from the original data on allele frequencies and sample sizes listed in the Appendixof JEHLE et al. 2001 , and the data for adults and larvae in 1998 were combined for the likelihood analysis. The generation interval was estimated to be 6 years, and the sampling interval is therefore 1.5 and 2 generations for ponds 1 and 2, respectively (JEHLE et al. 2001 ). Because of the unsynchronized sampling, it is not possible to apply our method for the two-small-population model. However, we can use our method for the infinite-source model as an approximation to estimate the effective size and immigration rate for one pond, with the other pond acting as the source population and samples from it pooled. Since only integer numbers for sampling intervals are allowed in likelihood estimation, the analysis for pond 1 was carried out by using a sampling interval of 2, and the estimates of N_e and m (N^'_e, m') were converted by _e = (1.5/2)N^'_e and = 1 - e^{(2/1.5)log(1-m')}, respectively. Similar results were obtained by using a sampling interval of 1.

The relative profile log-likelihood curves for the effective size and migration rate of subpopulations in ponds 1 and 2 are shown in Fig 6. The maximum-likelihood estimates are N_e = 16.3 and m = 0.39 for pond 1 and N_e = 18.0 and m = 0.19 for pond 2. If isolation between ponds is assumed, the likelihood N_e estimates are 17.1 and 21.9 for ponds 1 and 2, respectively, both being somewhat larger than those allowing migration.

View larger version (19K): In this window In a new window Download PPT slide   Figure 6. Relative profile log-likelihood curves for the effective size (top) and immigration rate (bottom) for subpopulations 1 (Ne1, m1) and 2 (Ne2, m2) of T. marmoratus. The data are published in JEHLE et al. 2001 .

The estimated migration rates between ponds are high. Although newts are found to migrate as juveniles and adults during breeding seasons (HAYWARD et al. 2000 ), there is no direct estimate of migration rate between populations in the literature to support or reject the estimates obtained herein. However, we argue that the migration rate could be high for this population. The average observed number of distinctive alleles per microsatellite locus is seven, which seems to be too large to be maintained in isolated populations with N_e at 10–20. Interspecies introgression may also be partially responsible for the high migration rate estimates. On average 4% of adult newts from these ponds are F₁-hybrids between T. cristatus and T. marmoratus (ARNTZEN and WALLIS 1991 ). Although hybrids were identified morphologically for adults and genetically with diagnostic microsatellites for larvae (JEHLE et al. 2000 ) and excluded from the analysis (JEHLE et al. 2001 ), it is inevitable that some hybrids may be undetected and included in the samples. Coincidentally, pond 2 contains few T. cristatus (JEHLE et al. 2001 ), and the immigration rate for T. marmoratus in this pond is estimated to be much lower than that in pond 1, which contains 75 T. cristatus individuals.

	DISCUSSION

TOP ABSTRACT DEFINITIONS AND ASSUMPTIONS IMMIGRATION FROM AN... TWO FINITE DEMES SIMULATION RESULTS APPLICATION TO A REAL... DISCUSSION APPENDIX LITERATURE CITED

The estimation of demographic properties of a species from its genetic characteristics is now an old topic. Many methods have been generated by statistical geneticists to estimate population sizes, migration rates, population growth rates, and many other interesting and useful ecological parameters. Some of these methods work quite well most of the time, and some are less successful. All share one property: they make simplifying assumptions about the possible range of demographies that can be considered, and the methods proposed here are no exception. It should be our goal, though