What is the migration of individuals between populations?

• Journal List
• HHS Author Manuscripts
• PMC1482677

Evolution. Author manuscript; available in PMC 2007 Jan 1.

Published in final edited form as:

Evolution. 2006 Jan; 60(1): 1–12.

PMCID: PMC1482677

NIHMSID: NIHMS8481

PMID: 16568626

Seongho Song,1,2 Dipak K. Dey,1 and Kent E. Holsinger3,4

Seongho Song

1Department of Statistics, U-4120, University of Connecticut, Storrs, CT 06269-4120

Find articles by Seongho Song

Dipak K. Dey

1Department of Statistics, U-4120, University of Connecticut, Storrs, CT 06269-4120

Find articles by Dipak K. Dey

Kent E. Holsinger

3Department of Ecology & Evolutionary Biology, U-3043, University of Connecticut, Storrs, CT 06269-3043

Find articles by Kent E. Holsinger

Disclaimer

1Department of Statistics, U-4120, University of Connecticut, Storrs, CT 06269-4120

3Department of Ecology & Evolutionary Biology, U-3043, University of Connecticut, Storrs, CT 06269-3043

2current address: Department of Mathematical Sciences, ML:0025, University of Cincinnati, Cincinnati, OH 45221-0025

4Address for correspondence: E-mail: ude.nnocu.bee.niwrad@tnek

Copyright notice

The publisher's final edited version of this article is available free at Evolution

Abstract

Populations may become differentiated from one another as a result of genetic drift. The amounts and patterns of differentiation at neutral loci are determined by local population sizes, migration rates among populations, and mutation rates. We provide exact analytical expressions for the mean, variance, and covariance of a stochastic model for hierarchically structured populations subject to migration, mutation, and drift. In addition to the expected correlation in allele frequencies among populations in the same geographical region, we demonstrate that there is a substantial correlation in allele frequencies among regions at the top level of the hierarchy. We propose a hierarchical Bayesian model for inference of Wright's F-statistics in a two-level hierarchy in which we estimate the among-region correlation in allele frequencies by substituting replication across loci for replication across time. We illustrate the approach through an analysis of human microsatellite data, we show that approaches ignoring the among region correlation in allele frequencies underestimate the amount of genetic differentiation among major geographical population groups by approximately 30%, and we discuss the implications of these results for the use and interpretation of F-statistics in evolutionary studies.

Keywords: F-statistics, genetic drift, migration, mutation, population structure

Sewall Wright pointed out more than 70 years ago that isolated populations tend to diverge from one another as a result of genetic drift and that the amount and pattern of divergence among populations reflects the extent of migration between them (). In particular, individuals belonging to the same population are expected to be more similar to one another than those belonging to different populations, and individuals from different populations within the same geographical region are expected to be more similar to one another than are individuals from different geographical regions. As a result, evolutionary biologists have commonly used statistics that describe patterns of population differentiation as an indication of the extent to which populations are genetically isolated from one another. Wright's F-statistics, and statistics related to them like Nei's G-statistics (), have been the most widely used descriptors of population differentiation for more than fifty years (; ).

While F-statistics are intended to measure the amount of differentiation among populations, they are implicitly based on Wright's island model (), a pseudo-spatial model in which all populations are equally likely to exchange migrants. Much work in theoretical population genetics has been directed towards understanding the patterns of genetic differentiation that arise when more realistic models of migration are incorporated into the models (e.g., ; Kimura and Weiss 1964; Nei and Feldman 1972; Felsenstein 1975; Nagylaki 1976; ; ; ; ; ). Many investigators now calculate pairwise F-statistics to provide an indicator of the amount of differentiation between those populations (e.g., ).

Another aspect of Wright's model has been less widely appreciated. His model implies that allele frequencies among populations at any given point in time are uncorrelated when the number of populations exchaning genes is very large. As a result, the variance in allele frequency among contemporaneous populations is equivalent to the variance in allele frequency within any one population over time. Surprisingly, this result no longer holds when a small to moderately large number of populations are exchanging genes.

show that substantial among-population correlations in allele frequency are expected even for mutation rates as large as 10−3 unless more than 100 populations are exchanging genes. Their results re-emphasize an obervation made by : When a finite number of populations is involved in gene exchange, the entire set of populations “drifts” together. The mean allele frequency changes from generation to generation, and the resulting correlation in allele frequencies among populations may be substantial unless the number of populations exchanging genes is very large. As a result, the variance in allele frequency among contemporaneous populations may be substantially smaller than the variance in allele frequency within any one population over time.

We explore two questions raised by these results: (1) How is the magnitude of among population correlation affected by hierarchical structure in the migration process? How does the magnitude of the within region correlation compare with the magnitude of the among region correlation? By hierarchical structure we mean that there is a higher rate of migration among populations within a region than among populations in different regions. (2) Given that the number of populations actually exchanging genes is larger than the number of populations from which samples are available and is generally unknown, can an approximate estimator that does not depend on the number of populations provide a reasonable estimate of the expected amount of differentiation among populations? After exploring these questions, we illustrate the importance of accounting for among population correlations in estimates of population genetic structure through analysis of a human population data set () with a hierarchical Bayesian model, and we discuss the implications of our results for genetic analyses of population structure.

PROCESS MODEL RESULTS

Analysis of population genetic data must incorporate two aspects of sampling: the usual statistical sampling associated with the sampling of alleles or genotypes within populations and the genetic sampling associated with sampling population allele or genotype frequencies from an underlying stochastic evolutionary process (see , pp. 15–17 for additional discussion). We use the modeling framework introduced in for the stochastic evolutionary process.

Briefly, we focus on a single locus with A allele types, b 1, b 2, ⋯, b A, and assume that there are P populations indexed by i (refer to Table 1 for a summary of notation used in the process model). Let VA×A, be a general mutation matrix which has a element, v rs, the probability of mutation from allele type b r to allele type b s. Let MP×P be a general (backward) migration matrix, i.e., Mij=m↼ij is the probability that the allele in population i came from population j. Backward migration rates are generally used for analysis of migration in population genetic models (see, for example, Nagylaki 1982; , ).

Table 1

Notation used in the single locus process model.

ParameterDefinitionANumber of allelesSNumber of populations in each geographical regionkNumber of geographical regionsPTotal number of populations. P = kSVAn A × A matrix of mutation rates. v rs is the mutation rate from the rth allele to the sth alleleMA P × P matrix of backward migration rates. m ij is the probability that an allele in population i came from population j in the preceding generationBThe Kronecker product of M and VN ijNumber of individuals in population j of geographical region ipij(t)

Open in a separate window

A multi-level hierarchy arises naturally when migration occurs more frequently among populations within the same geographical region than among populations in different geographical regions. Consider a 2-level hierarchy in which there are S populations nested within each of k geographical regions (kS = P). Assume that population j in geographical region i is of size N ij, and let pij(t) be the A × 1 vector of allele frequencies in that population at generation t. We concatenate the population allele frequency vectors within geographical region i to an SA × 1 vector pi(t) and concatenate the k resulting cluster vectors to a kSA × 1 vector p(t). We define

p∗(t) = (M ⊗ V′)p(t) , 

(1)

where M ⊕ V′ denotes the Kronecker product of M and V′. Given N ij, we assume that the population is diploid (so that the number of allele copies is 2N ij) and that the pij(t+1) are conditionally independent. Thus,

2Nijpij(t+1)∼M(2Nij,pij∗(t)),

(2)

where M denotes a multinomial distribution.. Through (1) and (2), we pass from p (t) → p*(t) → p (t+1).

Stationary equations for means and covariances

The Markov Chain defined through (1) and (2) has a finite state space. If all entries of V′ are nonzero for some t ≥ 1, then this chain has no absorbing states. In fact, it is aperiodic and irreducible and thus has a unique stationary distribution (see for an introduction to stationary distributions in population genetic models). The stationary mean vector is identical in all populations and is given by the unique left eigenvector of V corresponding to an eigenvalue of 1. The stationary variance-covariance matrix is given by the solution of the following set of equations:

Σij,ij=(1−12Nij)(BΣB′)ij,ij+12Nij[Diag(uij)−uiju′ij]Σij,i′j′=(BΣB′)ij,ij,

(3)

where B = M ⊕ V′ (see for details of the derivation, which follows ).

Results for a Finite Island Model

To make further analytical progress we consider a special case of the general structured migration model presented in the preceding section. The finite island model studied in imagined that a single backward migration rate applied to all populations. A natural generalization is to consider a hierarchical model in which migration occurs among populations within a single geographical region at the same rate but in which migration occurs among populations in different geographical regions at a different (and smaller) rate. For notational simplicity, we consider a 2-level hierarchy in which there are S populations within each of k geographical regions. We specify the migration matrix, M as follows:

M=(M1M2…M2M2M1…M2⋮⋮⋱⋮M2M2…M1),

(4)

where

M1=(m0m1S−1…m1S−1m1S−1m0…m1S−1⋮⋮⋱⋮m1S−1m1S−1…m0)

(5)

and

M2=(m2S(k−1)m2S(k−1)…m2S(k−1)m2S(k−1)m2S(k−1)…m2S(k−1)⋮⋮⋱⋮m2S(k−1)m2S(k−1)…m2S(k−1)),

(6)

with m 0 = 1 − m 1 − m 2, and m 1 > m 2. The remainder of the formulation is unchanged.

Lengthy calculations (outlined in ) provide exact expressions for the stationary mean vector and for the stationary variances and covariances. To illustrate the properties of this specialization we consider an example with two alleles and asymmetric mutation.

Example

With two alleles at a locus, the mutation matrix is given by

V=(v11v12v21v22),

where v 11 = 1 − v 12 and v22 = 1 − v 21. This mutation matrix leads to the stationary mean vector u = (μp, 1 − μp), where μp = ν21/(ν12 + ν21). In addition,

σ2=12Nμp(1−μp)1−(1−12N)δv2[(1−r1,S−r2,k+r3,k,s)+(r1,S−r3,k,S)ρ1+r2,kρ2],

where δv = v 11 − v 21 = 1 − (ν12 + v 21), σ2 is the variance in allele frequency within any one population over time, ρ1 is the allele frequency correlation for populations in the same geographical region, and ρ2 is the allele frequency correlation for populations in different geographical regions. ρ1 and ρ2 are given by

ρ1=δv2(r1,S−r3,k,S)[1−δv2(1−r2,kk−1)]c(S−1)+δv4r2,k2cS(k−1)ρ2=r2,kδv2[1−δv2(1−r1,S−r3,k,SS−1−r2,k)]cS(k−1)+(r1,S−r3,k,S)r2,kδv4cS(k−1)

where

c=[1−δv2(1−r2,kk−1)][1−δv2(1−r1,S−r3,k,SS−1−r2,k)]−(S−1)δv4S(k−1)r2,k2r1,S=2m1−SS−1m12r2,k=2m2−kk−1m22r3,k,S=2m1m2−S−1S(k−1)m22.

In the case where each geographical region contains only a single population, i.e., S = 1, these results reduce to those presented in for the finite island model.

Figure 1 illustrates the magnitude of ρ1 and ρ2 for several parameter combinations. Notice that with mutation rates that may be typical of protein-coding loci (v = 5 × 10−7), both the correlation of allele frequencies among populations within a geographical region, ρ1, and the correlation among geographical regions, ρ2, are very high unless either the number of populations within each region, S, or the number of regions K is on the order of 200 or more. Moreover, the correlation decays more slowly as a function of increasing numbers of populations and regions when migration is common (m 1 = 0.04, m 2 = 0.01) than when it is rare (m 1 = 0.008, m 2 = 0.002). Notice also that the within-region correlation is relatively insensitive to the number of regions. It depends largely on the number of populations within each region. Similarly, the among-region correlation depends primarily on the number of regions and less on the number of populations within regions. Finally, with mutation rates that may be more typical of microsatellites (v = 0.005) the within-region correlation remains substantial with as many as 40 populations per region and the among-region correlation is small only when migration is rare (m 1 = 0.008, m 2 = 0.002).

Open in a separate window

Figure 1

Plots of the correlation among populations within the same geographical region, ρ1, and among geographical regions, ρ2, as a function of the number of geographical regions, k, and the number of populations within each region, S, for several choices of migration rate with two mutation rates in the two-allele model.

INTERPRETATION FOR F-STATISICS

In general, the expressions for the stationary variance depend on the stationary mean. We can, however, calculate a scaled variance,

FST=σp(t)2μp(t)(1−μp(t)),

(7)

that removes this dependence. In this case F ST can be regarded as an intraclass correlation coefficient (; ). In fact, , p.294) defined F ST as “the correlation between random gametes within populations, relative to gametes of the total set of populations.” Since it was first introducted by and , F ST has been the most widely used statistic for describing hierarchical structure in genetic data.

For a finite set of populations evolving according to (1) and (2) there are two quantities that might correspond with (7). (Refer to Table 2 for a summary of notation used in relation to F-statistics.) One of these quantities, which we will denote θ, takes σp(t)2 to be the allele frequency variance in a single population across time (or equivalently as the allele frequency variance in a single population associated with different realizations of the stochastic evolutionary process). The second takes σp(t)2 to be the allele frequency variance across populations at a single point in time, i.e.,

θ(p11(t),…,pkS(t))=∑i=1k∑j=1S(pij(t)−μp(t))2∕kSμp(t)(1−μp(t)),

(8)

where μp(t)=∑i=1k∑j=1Spij(t)∕kS.

Table 2

Notation used in the two-allele, multiple locus hierarchical beta model.

ParameterDefinitionθ(I)An estimate of θθ(II), θ(III)Alternative estimates of θ(p11(t),…,pkS(t))INumber of loci (two alleles per locus)KNumber of geographical regionsS kNumber of populations in geographical region kp iksFrequency of allele A 1 at locus i in population s of geographical region kπikMean allele frequency at locus i in geographical region kπiMean allele frequency at locus iπMean allele frequency across lociπik(1 − πik)θxVariance in allele frequency at locus i within geographical regionsπi(1 − πi)θyVariance in allele frequency at locus i among geographical regionsπ(1 − π)θvVariance in allele frequency across lociρ1Allele frequency correlation among populations within geographical regionsρ2Allele frequency correlation among geographical regions

Open in a separate window

We use θ in our notation to emphasize that both quantities treat populations as random effects, as in Cockerham's random-effects model for θ-statistics (; ; ; ). In particular, θ in our formulation has precisely the same interpretation as θ in the Cockerham random-effects model (Weir, personal communication; ). Moreover, θ is the parameter most directly related to features of the allele frequency distribution associated with the process model described above, as well to allele frequency distributions arising from the pure isolation model underlying evolutionary interpretations of Weir and Cockerham's θ-statistics (; ; ). Thus, we regard θ as the quantity that corresponds to Wright's F ST. Notice that θ(p11(t),…,pkS(t)) corresponds to a random effects interpretation of G ST (1973).

Estimation of θ is straightforward. Given a sample of genetic data from some set of populations we estimate μp(t) and σp(t)2 from the data, and the estimate corresponding to θ is θ(I)=σ^p(t)2μ^p(t)(1−μ^p(t)). Notice that if the number of populations and population clusters is very large, then the allele frequency correlations illustrated in Figure 1 are negligible, and θ≈θ(p11(t),…,pkS(t)) (see also ). Recalling that k and S in (8) refer to the number of population clusters and populations among which migration is occurring, not the number of populations included in the sample, an estimate of θ(I) might also be regarded as an estimate of θ(p11(t),…,pkS(t)) when ρ1 and ρ2 are small.

Estimation of θ(p11(t),…,pkS(t)) is less straightforward. Let Num denote the numerator on the right hand side of (8), and Denom the corresponding denominator. θ(III) = E(Num/Denom) would correspond precisely with the definition of θ(p11(t),…,pkS(t)). To estimate θ(III), however, would require knowledge of k and S, and they are generally unknown. An alternative is to estimate θ(p11(t),…,pkS(t)) with θ(II) = E(Num)/E(Denom), which is is relatively insensitive to k and S for moderate values. In fact, we can show that θ(II) → θ(III) in probability as k and S tend to infinity. Furthermore, by using θ(II) to estimate θ(p11(t),…,pkS(t)), we can provide an approximate expression for θ(II) that does not depend on k or S (see ).

Tables ​Tables33 and ​and44 provide comparisons among θ(I), θ(II) and θ(III) for 2N = {100, 1000}, v = v 12 = v 21 = {5 × 10−3, 5 × 10−6}, and several choices of m 1 and m 2 with m 1 + m 2 = 0.05. We calculate the expectations for θ(III) from simulations using the process model in (1) and (2) conditional on 0<(1∕kS)Σijpij(t)<1; efficiency refers to the fraction of simulated samples satisfying this condition. We calculate θ(I) and θ(II) from the corresponding analytical expressions. Notice that θ(I) can be two orders of magnitude greater than θ(II) and θ(III) under model (1) and (2). Notice also that θ(II) and θ(III) are quite similar, except when the population size, N is small. These patterns also hold with larger local population sizes (results not shown). For large mutation rates (5 × 10−3) the differences among θ(I), θ(II), and θ(III) are fairly small. For smaller mutation rates (5 × 10−6) the differences between θ(II) and θ(III) are negligible, except when local populations are small (N = 50), but the differences between either of them and θ(III) can be substantial. With N = 500, k = S = 5, m 1 = 0.04, and m 2 = 0.01, for example, θ(II) ≈ θ(III) = 0.01, but θ(I) = 0.67.

Table 3

Comparison of statistics related to F ST under models (1) and (2) with migration rate = 5 × 10−2 and mutation rate = 5 × 10−3. θ(I) measures the variance in allele frequency over time. θ(III) measures the variance in allele frequency among populations at one time. θ(II) is an approximation of θ(III) that does not require the number of populations and geographical regions involved in gene exchange to be known.

2NkSm1m2θ(I)θ(II)θ(III)SDefficiency100550.040.010.1048570.0882870.088165(0.0261710)1.00000.0450.0050.1149850.0987910.098595(0.0319790)1.0000250.0921990.0751520.075144(0.0196840)1.0000

Open in a separate window

Table 4

Comparison of θ(I), θ(II), and θ(III) under models (1) and (2) with migration rate = 5 × 10−2 and mutaion rate = 5 × 10−6.

2NkSm1m2θ(I)θ(II)θ(III)( SD )efficiency100550.040.010.9530720.1184880.035089( 0.054724 )0.47550.0450.0050.9532080.1649430.033850( 0.065785 )0.5904250.9530470.0872270.062394( 0.041973 )0.1781

Open in a separate window

INFERENCE UNDER A HIERARCHICAL BETA MODEL

Results presented in the preceding sections refer to variation at a single locus observed across time, but as evolutionary biologists we more commonly have samples from several or many different genetic loci collected from populations at the same time. If we regard the genetic sample as having been drawn from populations at stationarity with respect to the same underlying stochastic evolutionary process, then the samples at each locus can be regarded as (approximately independent) samples from the stationary distribution of this process. In fact, even if the populations are not at stationarity it may be reasonable to regard samples from different loci as (approximately independent) samples from a common underlying process, provided that the relevant evolutionary parameters (mutation, migration, and population size) are shared across the loci.

Assuming a common N and m across loci seems reasonable, because these parameters refer to demographic features of the populations from which the samples are drawn. It may seem problematic to assume common mutational parameters, but notice that in the example studied above, the scaled variance and correlations depend only on the total mutation rate, ν11 − ν21 = 1 − (ν12 + ν21). Similarly, in a multiallele model, the scaled variance and correlations depend only on the total mutation rate, provided that all alleles are equally likely to mutate (see equations (12) and (13) in ). Thus, we need not assume a common mutational bias across loci, only that the total mutation rate at each locus can be regarded as an independent sample from the same underlying distribution of mutation rates. In short, we propose to substitute replication across loci for replication across time. , pp. 174–175) justified bootstrap resampling across loci to obtain the sampling distribution for a method-of-moments estimate for θ using a similar argument (see also ).

We assume that our sample consists of data from K population clusters, each of which has S k populations, k = 1, … , K, and that the genotype of each individual is determined at I loci, each with two alleles A 1 and A 2. Let pIKS×1 denote the allele frequencies of A 1, i.e., p iks is the allele frequency of A 1 in s th population of the k th population cluster at locus i, i = 1, … , I; k = 1, … , K; s = 1, … , S k. Notice that this structure corresponds to the stochastic evolutionary process illustrated in our example. Classical results (summarized in and in ) suggest that the stationary distribution of allele frequencies among populations is well-approximated by a Beta distribution when only two alleles are segregating. ; see also ) extends these results to show that the distribution of allele frequencies is well-approximated by a Dirichlet distribution when mutation is equally frequency among all alleles at a multi-allelic locus. Moreover, has shown that a Dirichlet distribution arises naturally in population genetics whenever alleles are regarded as exchangeable, regardles of whether the underlying stochastic evolutionary process has reached stationarity. Several authors have used Beta or Dirichlet distributions to construct likelihood or Bayesian methods for analysis of genetic structure (e.g., ; ). Thus, we assume a hierarchical beta model for the allele frequency distribution:

piks∣πik,θx∼Beta(1−θxθxπik,1−θxθx(1−πik)),πik∣πi,θy∼Beta(1−θyθyπi,1−θyθy(1−πi))

(9)

and

πi∣π,θv∼Beta(1−θvθvπ,1−θvθv(1−π)).

It follows that

E(piks)=π,σ2=Var(piks)=π(1−π){θv+θy(1−θv)+θx(1−θy)(1−θv)},σ2ρ1=Cov(piks,piks′)=π(1−π){θv+θy(1−θv)}

for s ≠ s′,

σ2ρ2 = Cov(piks, pik′s′) = π(1 − π)θv

for k ≠ k′ and s ≠ s′ and

Cov(piks, pi′k′s′) = 0

for all k, k′, s and s′ but i ≠ i′. Extension of this modeling framework to multiple alleles per locus is straightfoward, but tedious (see for an outline of the approach). describe an alternative approach with the same variance-covariance properties based on a mixture of beta random variables. We assume the prior distribution for the parameters is given by independent uniform densities. Posterior inference is done through Markov Chain Monte Carlo (MCMC) simulation with the Metropolis-Hastings algorithm ().

FST in the Hierarchical Beta Model

Natural estimates of θ and θ(p111(t),…,pIKS(t)) arise from this model. Specifically,

θ(I)=σp(t)2μp(t)(1−μp(t))=σ2π(1−π)=θv+θy(1−θv)+θx(1−θy)(1−θv),

corresponds directly with the definition of θ as a function of the temporal mean and variance in allele frequencies. Notice also that

1 − θ(I) = (1 − θv)(1 − θy)(1 − θx) .

Thus, θ(I) exhibits the conventional mathematical properties of F-statistics in hierarchically structured populations (, p. 295; ). θv corresponds to the scaled allele frequency variance across loci, θy to the scaled allele frequency variance across geographic regions within a locus, and θx to the scaled allele freqency variance across populations within a geographic region.

An estimate of θ(p111(t),…,pIKS(t)) corresponding to θ(II) arises from noticing that as the number of populations and geographical regions tends to infinity, for i = 1, … , I,

E[∑k=1K∑s=1Sk(piks−p¯...)2∕KS¯]E[p¯...(1−p¯...)]

converges in probability to

θ(II)=θ(I)[1−(w1ρ1+w2ρ2)]1−θ(I)(w1ρ1+w2ρ2),

with w1=∑k=1KSk2∕K2S¯2, w 2 = 1 − w 1 and S¯=∑k=1KSk∕K where the expectation is across the stochastic evolutionary process. Notice that ∑k=1KSk2=∑k=1K(Sk−S¯)2+KS¯2. Thus, w 1 tends to 0 as K tends to infinity, provided that there exists ε in (0, ∞) such that ∑k=1K(Sk−S¯)2<ε for all K. As a result,

θ(II∗)=θ(I)(1−ρ2)1−θ(I)ρ2=θy+θx(1−θy)

(10)

provides a reasonable estimate of E(Num)/E(Denom) when the number of populations and geographical regions is unknown but large (compare ). Finally,

θ(III)=E[∑i=1I∑k=1K∑s=1Sk(piks−p¯...)2∕(IKS¯)p¯...(1−p¯...)]

where the p iks are replaced with samples from the posterior distribution and the expectation is calculated across the joint posterior corresponds directly with the definition of θ(p111(t),…,pIKS(t)), if K and S k are known. For completeness, we note that

ρ1=θv+θy(1−θv)θv+θy(1−θv)+θx(1−θy)(1−θv)

and

ρ2=θvθv+θy(1−θv)+θx(1−θy)(1−θv).

Decomposition of θ(p111(t),…,pIKS(t))

In the formulation just presented θ(p111(t),…,pIKS(t)) summarizes the population structure while ignoring the hierarchy. As is clear from the discussion of θ(I), however, it is natural to consider the partitioning of variation at each level of the hierarchy. There we pointed out that θx corresponds to the scaled allele frequency variance among populations within geographical regions and θy to the scaled allele frequency variance among geographical regions. In this section we consider an approach to decomposing θ(p111(t),…,pIKS(t)). Specifically, let

θK(II∗)=E[∑k=1K∑s=1SkSk(p¯ik.−p¯i..)2∕KS¯]E[p¯...(1−p¯...)]=θ(I)(ρ1−ρ2)1−θ(I)ρ2=θy

and

θE(II∗)=E[∑k=1K∑s=1Sk(piks−p¯ik.)2]∕KS¯E[p¯...(1−p¯...)]=θ(I)(1−ρ1)1−θ(I)ρ2=θx(1−θy),

where the expectations are taken across the stochastic evolutionary process. Then θK(II∗) and θE(II∗) correspond to the scaled allele frequency variance among contemporeneous geographical regions and among contemporaneous populations within geographical regions, respectively. Notice that

θ(II∗)=θK(II∗)+θE(II∗).

Thus, θK(II∗) and θE(II∗) provide a readily interpretable additive decomposition of θ(II*). Specifically, θK(II∗)=θy is the proportion of contemporaneous genetic variation attributable to allele frequency differences among regions, and θE(II∗) is the proportion of contemporaneous genetic variation attributable to allele frequency differences among populations within regions. θx and θy retain the multiplicative relationship with θ(II*) characteristic of Wright's F-statistics:

1 − θ(II∗) = (1 − θx)(1 − θy).

APPLICATION

Human Population Data

describe the HGDP-CEPH Human Genome Diversity Cell Line Panel. We analyze the genetic diversity at 377 autosomal microsatellite loci in 1056 individuals included in the panel. These individuals represent 52 worldwide populations, and the genotypes are derived from those used in a related analysis reported in . Using an individual-based Bayesian classification algorithm (), were able to gather the 52 populations into 5 geographically coherent groups. We treat these groups as geographical regions (K = 5) in our analysis: Eurasia (21 populations, i.e., S 1 = 21), Africa (6 populations, i.e., S 2 = 6), East Asia (18 populations, i.e., S 3 = 18), America (5 populations, i.e., S 4 = 5) and Oceania (2 populations, i.e., S 5 = 2). To illustrate analysis within the framework developed above, we designate the most frequent allele type at each locus as A 1 and group all other allele types into a pseudo-allele type A 2. Analyses in which we designated the second most frequent allele type at each locus A 1 and grouped the remainder into a pseudo-allele class produced results concordant with those described below (results not shown).

Results

Results from analysis of these data are presented in Table 5. To illustrate the effect of ignoring the region-level allele frequency correlation Table 5 includes estimates from a model in which we assume independent uniform (0,1) prior distributions for πi. Specifically, the hierarchical model developed above uses

πi∣π,θv∼Beta(1−θvθvπ,1−θvθv(1−π)),

while the reduced model uses

πi ∼ Unif(0, 1) .

Because the actual number of population clusters and populations is unknown, we present values for θ(III) in which we assume that the actual numbers of population clusters and populations are equal to those in our sample.

Table 5

Parameter estimates for the human microsatellite data. See the text for an explanation of the difference between the two models.

Hierarchical beta modelReduced modelParameterMean(2.5%, 97.5%)Mean(2.5%, 97.5%)π0.3401(0.3301, 0.3506)n.a.1θx0.0294(0.0283, 0.0306)0.0293(0.0283, 0.0305)θy0.0784(0.0729, 0.0845)0.0760(0.0704, 0.0818)θv0.0298(0.0237, 0.0370)n.a.ρ10.8010(0.7869, 0.8149)0.7371(0.7191, 0.7552)ρ20.2253(0.1847, 0.2708)n.a.

Open in a separate window

1Parameters not estimated in the reduced model are reported as not applicable.

2θ(III) based on the observed number of population clusters and populations.

Notice that in the hierarchical model θ(II*) = 0.1055 is quite similar to θ(III) = 0.1046, suggesting that our estimate of θ(p11(t),…,pkS(t)) is little affected either by uncertainty about the actual number of population clusters and populations or by approximating E(Num/Denom) by E(Num)/E(Denom). Similarly, notice that in the reduced model, θ(I) = θ(II*) = 0.1031 is quite similar to θ(III) = 0.1046 in the hierarchical model. More importantly, notice that θ(I) = θ(II*), because the allele frequency correlation among geographical regions is neglected. As a result, the estimate of θ(I) = 0.1031 from the reduced model is about 30% smaller than the estimate of θ(I) = 0.1322 in the hierarchical model. Method of moments estimates for the full, multi-allele data set (reported in ) give an estimate quite similar to those reported here for the variance component among populations within regions: 0.025 versus θx = 0.0294. The smaller value reported for the variance component among regions, 0.043 versus θy = 0.0784, probably reflects the effect of grouping less frequent alleles into a single pseudo-allele class.

The discrepancy between estimates of F ST in the hierarchical and reduced models (θ(I) = 0.1322 versus θ(I) = 0.1031, respectively) illustrates that analyses of F-statistics must either incorporate among region correlations or interpret the results in terms of differentiation relative to an unknown set of geographcial covariances (see for similar results in an analysis of variance context). In these data, the price for ignoring the population covariances is a 30% underestimate of F ST . In a similar analysis of data from Platanthera leucophaea, an orchid, showed that neglecting the among population correlations would lead to an estimate of F ST less than half of the best estimate (0.2122 versus 0.4837).

These data also show that there are substantial correlations among both populations (ρ1 = 0.8010) and geographical regions (ρ2 = 0.2253). Moreover, the additive partitioning of θ(II*) indicates that 75% of the total allele frequency variation is a result of differences among geographic regions while 25% is a result of differences among populations within geographical regions (θK(II∗)=0.0784,θE(II∗)=0.0271). Notice both models lead to the same additive partitioning, because it is done relative to the set of geographical covariances.

Alternatively, we can use the traditional multiplicative partitioning of F-statistics and note that there is more than twice as much allele frequency differentiation among population clusterss than among populations within population clusters (θy = 0.0784 vs. θx = 0.0298). For these data, the two partitionings produce similar numerical results. For data sets in which θy is large, however, the differences may be substantial. Which partitioning is to be preferred will depend on whether interest centers on describing the contribution of each level of the hierarchy to contemporaneous variation in allele frequency (θK(II∗) and θE(II∗)) or on describing the magnitude of allele frequency variance at each level across replicates of the stochastic evolutionary process (θy and θx).

DISCUSSION

Wright's F-statistics have been the most widely used tool for describing the genetic structure of populations since they were first introduced more than fifty years ago (; ). Unfortunately, they paid little attention to the distinction between population and sample frequencies. clearly understood the distinction between population and sample frequencies and the additional role of genetic sampling, but gene diversity statistics (1973) were the first widely used method that recognized some of the sampling properties. Unfortunately, Nei's statistics accounted only for the statistical sampling associated with sampling from actual populations and did not consider the genetic sampling associated with sampling the populations from an underlying stochastic evolutionary process (). To address these shortcomings, introduced what is now the most widely used approach for characterizing population structure, a method-of-moments approach that accounts for both statistical and genetic sampling.

While the approach is satisfactory in many ways, it suffers from several important limitations. First, interval estimates of the parameters either depend on asymptotic convergence of mean squares calculated from Bernoulli random variables to a chi-squared distribution () or on bootstrapping of samples across loci. In contrast, interval estimates are easily derived from the MCMC sample produced during analysis of our hierarchical model. Second, the chi-squared approximation depends on the assumption that the correlation among populations is zero, while the results in section 2 clearly show that the among-population correlation can be very large for reasonable mutation and migration rates, unless the number of populations and population clusters is very large. Third, as point out, the estimate of θ provided by Weir and Cockerham's method estimates F ST only when the correlation among populations is zero. When the correlation is non-zero, then the estimate provided is

β=θ−θA1−θA,

where θ is equivalent to F ST and π(1 − π)θA is the covariance in allele frequencies among populations. If samples at each locus are independent samples from a common underlying stochastic evolutionary process, however, then θA is equivalent to our θv and β is equivalent to our θ(II*). Specifically,

θ(II∗)=θ(I)−θv1−θv,

where θ(I) is our estimate of F ST. In short, the parameters of the hierarchical statistical model we develop have direct, natural interpretations in terms of parameters related to the underlying stochastic evolutionary process, θ and θ(p11(t),…,pkS(t)), even in the face of substantial among population correlations. Ours is not the only approach that can account for these correlations, but existing alternatives depend on approximating the among population allele frequency distribution with a normal distribution (; ) or use a mixture approach whose interpretation is less natural (). Of course, a normal approximation will work well only when allele frequency distributions are unimodal, roughly symmetric, and have little mass near 0 or 1, conditions satisfied only when min [((1 − θ) ∕ θ)p, ((1 − θ) ∕ θ)(1 − p)] ≫ 0.

Applying our method to human microsatellite data illustrates that the among population correlation in allele frequencies can be substantial, even in species like humans with large, highly dispersive populations. The estimated allele frequency correlation among geographical regions (Eurasia, Africa, East Asia, America, and Oceania), ρ2, is greater than 20%. The estimated allele frequency correlation among populations, ρ1, is even greater, about 80%. As a result of the correlation among major population groupings our best estimate of F ST, θ(I), is about 30% greater than estimate of F ST that ignore the region-level correlation..

While the numerical difference betewen θ(I) and θ(II*)) (0.13 vs. 0.10) may not seem great, one reason for interest in estimates of θ is their possible relationship to evolutionary process parameters. For example, if we were to treat the data as coming from a 1-level hierarchy, assume that mutation rates are much smaller than migration rates, and neglect the among population correlation in allele frequencies, then F ST = θ ≈ 1/(4Nm + 1) (). Thus, we might estimate Nm as Nm^=(1∕f−1)∕4. Using θ(I) as our estimate of F ST we would obtain Nm^=2.25, while using θ(II*) we would obtain Nm^=1.67.

As the results in Tables ​Tables33 and ​and44 show, however, we might vastly overestimate the extent to which gene flow is responsible for retarding population differentiation when the actual number of population clusters and populations is small or the migration rates are large or both. In the case of Platanthera leucophaea, for example, ignoring the correlation among populations would lead to an inferred estimate of Nm of about 0.98 with a 95% credible interval of 0.73-1.20, while a better estimate would be that Nm is only 0.27 with a 95% credible interval of 0.19-0.37 (). In short, because existing methods for analysis of population structure typically neglect among population correlations in allele frequency and because such correlations may be substantial, past inferences about migration derived from F-statistics may substantially overstate the importance of gene exchange in preventing differentiation of partially isolated populations.

ACKNOWLEDGEMENTS

Comments from Jody Hey and two anonymous reviewers led to substantial improvements in this paper. We are indebted to them for the time and energy they devoted to reading and commenting on an earlier draft. This research was supportted in part by a grant from the U.S. National Institutes of Health, 1 R01 GM068449-01A1.

APPENDIX A

At stationarity, the distribution for p(t) is the same for all t. Thus, let ukSA×kSA denote the stationary mean vactor, i.e., E[p(t)|M, V, N] where N is the vector of population sizes.

Since E[pij(t+1)∣pij(t),M,V,N]=pij∗(t),

E[p(t+1)∣M,V,N]=EE[p(t+1)∣p(t),M,V,N]=E[p∗(t)∣M,V,N]=(M⊗V′)E[p(t)∣M,V,N]=B(M,V)E[p(t)∣M,V,N],

(A.1)

where B(M, V) = (M ⊗ V′). Under stationarity, u = B(M, V)u. In fact, u does not depend on M. It is the unique left eigenvector of V corresponding to an eigenvalue of 1 (see for details).

Turning to stationary covariances, let Σ be the kSA dimensional stationary variance-covariance matrix of p. We define Δ=diag(Δ11, …, ΔkS) where

Δij=[Diag(pij∗(t))−pij∗(t)(pij∗(t))′]∕2Nij,

with Diag(pij∗(t))=Diag{pij1∗(t),…,pijA∗(t)}. For convenience, we denote E(Δ) = diag(E(Δ11), …, E(ΔkS)). Then

ΣkSA×kSA=Var(p(t+1)∣M,V,N)=BVar(p(t)∣M,V,N])B′+E(Δ∣M,V,N)

(A.2)

since P*(t) = Bp(t) and pij∗(t)=∑l=1k∑m=1SBij,lmplm(t) where B = B(M, V) and Bij,lm denote the (ij, lm)th A × A block matrix of B.

We calculate E[Diag(pij∗(t))∣M,V,N] and E[pij∗(t)(pij∗(t))′∣M,V,N] to obtain E(Δij|M, V, N). First, we observe that

E(Δij∣M,V,N)=Diag(∑l=1k∑m=1SBij,lmE[plm(t)∣M,V,N]).

Next,

E[pij∗(t)(pij∗(t))′∣M,V,N]=Var[pij∗(t)∣M,V,N]+E[pij∗(t)∣M,V,N]{E[pij∗(t)∣M,V,N]}′=∑l=1k∑m=1S∑l′=1k∑m′=1SBij,lmCov[plm(t),pl′m′(t)∣M,V,N]B′ij,l′m′+Bij,lmE[plm(t)∣M,V,N]{E[pl′m′(t)∣M,V,N]}′B′ij,l′m′.

Therefore,

E[Δij∣M,V,N]=12Nij{Diag(∑l=1k∑m=1SBij,lmE(plm(t)∣M,V,N))−Bij,lmE[plm(t)∣M,V,N]{E[pl′m′(t)∣M,V,N]}′B′ij,l′m′−∑l=1k∑m=1S∑l′=1k∑m′=1SBij,lmCov[plm(t),pl′m′(t)∣M,V,N]B′ij,l′m′}

(A.3)

Under stationarity, E(Δij|M, V, N) becomes

12Nij[Diag(uij)−uiju′ij−(BΣB′)ij,ij],

(A.4)

since u = Bu. Thus, for i, i′ = 1, …, k and j, j′ = 1, …, S we have the A × A diagonal matrix of Σ,

Σij,ij=(1−12Nij)(BΣB′)ij,ij+12Nij[Diag(uij)−uiju′ij]

(A.5)

and the A × A diagonal matrix,

Σij,i′j′ = (BΣB′)ij,ij, 

for i ≠ i′ or j ≠ j′.

APPENDIX B

Recalling that V is the mutation matrix, we let

B(M,V)kSA×kSA=M⊗V′=(M1⊗V′M2⊗V′…M2⊗V′M2⊗V′M1⊗V′…M2⊗V′⋮⋮⋱⋮M2⊗V′M2⊗V′…M1⊗V′).

Similarly,

B′ = (M⊗V′)′ = M ⊗ V , 

since M is symmetric matrix. Thus ∑kSA×kSA reduces to

ΣkSA×kSA=(Σ11Σ12…Σ12Σ12Σ11…Σ12⋮⋮⋱⋮Σ12Σ12…Σ11)

For simplifying the , we need to get (B∑B′) which consists of the following matrices, (M i ⊗ V′) ∑i′ j′ (Mj ⊗ V′) where i, i′, j, j′ = 1, 2:

(M1 ⊗ V′)∑11(M1 ⊗ V′) consists of the following A × A submatrices: diagonal matrix is

V′[(m02+m12S−1)Σ111+(2m0m1+(S−2)m12S−1)Σ211]V

and off-diagonal matrix is

V′[(2m0m1S−1+(S−2)m12(S−1)2)Σ111+{m12S−1+(m0+(S−2)m12(S−1)2)2}Σ211]V.

Further

(M1⊗V′)Σ11(M2⊗V′)=V′[m2(m0+m1)S(k−1)Σ111+(S−1)m2(m0+m1)S(k−1)Σ211]V=(M2⊗V′)Σ11(M1⊗V′),(M2⊗V′)Σ11(M2⊗V′)=V′[m22S(k−1)2Σ111+(S−1)m22S(k−1)2Σ211]V,(M1⊗V′)Σ12(M1⊗V′)=V′[(m0+m1)2Σ312]V,(M1⊗V′)Σ12(M2⊗V′)=V′[(m0+m1)m2k−1Σ312]V=(M2⊗V′)Σ12(M1⊗V′)

and

(M2⊗V′)Σ12(M2⊗V′)=V′[m22(k−1)2Σ312]V.

After a straightforward but tedious calculation, it can be shown that

(BΣB′)111=V′[(1−r1,S−r2,k+r3,k,S)Σ111+(r1,S+r3,k,S)Σ211+r2,kΣ312]V,(BΣB′)211=V′[r1,S−r3,k,SS−1Σ111+(1−r2,k−r1,S−r3,k,SS−1)Σ211+r2,kΣ312]V

and

(BΣB′)312=V′[r2,kS(k−1)Σ111+S−1S(k−1)r2,kΣ211+(1−r2,kk−1)Σ312]V.

Finally, we obtain

Σ111=(1−12N)V′[(1−r1,S−r2,k+r3,k,S)Σ111+(r1,S−r3,k,S)Σ211+r2,kΣ312]V+12N(Diag(u11)−u11u11′),Σ211=V′[r1,S−r3,k,SS−1Σ111+(1−r2,k−r1,S−r3,k,SS−1)Σ211+r2,kΣ312]V

What is the movement of individuals between populations?

What is it called when people move into a new population?

What is the effect of migration between two populations?

When migration of a section of population to another place is called?