Abstract
We developed an empirical Bayes procedure to estimate genetic distances between populations using allele frequencies. This procedure makes it possible to describe the skewness of the genetic distance while taking full account of the uncertainty of the sample allele frequencies. Dirichlet priors of the allele frequencies are specified, and the posterior distributions of the various composite parameters are obtained by Monte Carlo simulation. To avoid overdependence on subjective priors, we adopt a hierarchical model and estimate hyperparameters by maximizing the joint marginallikelihood function. Taking advantage of the empirical Bayesian procedure, we extend the method to estimate the effective population size using temporal changes in allele frequencies. The method is applied to data sets on red sea bream, herring, northern pike, and ayu broodstock. It is shown that overdispersion overestimates the genetic distance and underestimates the effective population size, if it is not taken into account during the analysis. The joint marginallikelihood function also estimates the rate of gene flow into island populations.
AS a stock management tool to counteract decreased or depleted fishery resources, stock enhancement programs have been undertaken in many countries for salmonid (Hilborn and Winton 1993; Ritter 1997; Kaeriyama 1999; Knapp 1999) and for other marine species (Bartley 1999; Kitada 1999). Concerns about the genetic effects of hatchery releases on wild populations have increased and aroused discussion (Walters 1988; Waples 1991; Hilborn 1992; Utter 1998; Waples 1999). Campton (1995) reviewed the genetic effects of hatchery releases on natural stocks of salmon and brown trout and concluded that the empirical data supporting those arguments are absent or largely circumstantial. This is a complex topic that needs further research (Waples 1999). A 10point approach for a responsible stock enhancement program has been proposed, which includes the need to use genetic resource management to avoid deleterious genetic effects (Blankenship and Leber 1995). Using wild individuals as broodstock may possibly reduce genetic risks (Bartleyet al. 1995; Haradaet al. 1998).
The genetic identity between produced progenies and the wild stock will be required before one can release the progenies. To examine the genetic identity, statistically significant differences are required. The homogeneity χ ^{2} test of allele frequencies is commonly used for testing genetic differences and the Roff test (Roff and Bentzen 1989) is used when minor alleles exist. If the null hypothesis is not rejected, the statistical power is required to be reported from a conservation viewpoint (Peterman 1990; Dizonet al. 1995). However, when the genetic difference is small, the corresponding statistical power may also be small with small sample sizes, making it difficult to conclude that there is no genetic difference. The statistical power is the probability of detecting the alternative hypothesis when it is correct. A considerably large sample size is required if one wants to obtain satisfactory large statistical power to reject the null hypothesis and detect small genetic differences. The hypothesis testing framework implies rejecting the null hypothesis, so it does not work well for detecting the genetic identity, and calculating the power is meaningless. Problems of the null hypothesis testing framework are discussed in Cohen (1994) and Hagen (1997).
An effective method of determining genetic identity is to examine the genetic distances between populations. If an estimated confidence interval of the genetic distance between two populations includes 0, we could conclude that the populations are genetically identical or not statistically significantly different. There are several measures for the genetic distance (Nei 1987). However, the sample distributions of these genetic distances are unknown. It is then inappropriate to estimate the confidence intervals of the genetic distance using asymptotic variances of the estimator.
In this article, we develop a Bayesian estimating procedure to measure genetic distances between populations from allele frequencies. We can directly evaluate the probability distribution of the genetic distance from its posterior distribution. The general method developed here is extended to estimate the effective population size termed N_{e} in a population based on temporal changes in allele frequencies. The joint marginallikelihood function derived here coincides with the likelihood function to estimate the rate of gene flow into island populations using the sample allele frequency from a number of islands (Rannala and Hartigan 1996).
METHODS
Let the frequencies of k alleles of two populations to be compared be p_{11}, · · ·, p_{1}_{k} and p_{21}, · · ·, p_{2}_{k}. Sanghvi (1953) proposed that the genetic distance between two populations be determined by
Prior and posterior distribution of D: It is not easy to describe a reasonable prior distribution of D, especially when we compare more than two populations. Alternatively, we set a prior for allele frequencies. Let the allele frequency of a population be p = (p_{1}, · · ·, p_{k})′ and the sample count be n = (n_{1}, · · ·, n_{k})′, where Rp_{i} = 1 and Rn_{i} = 2n (n individuals). When the sample is collected by a simple random sampling procedure with replacement, n follows a multinomial distribution. A β distribution is known as a conjugate prior of the binomial parameter p. A Dirichlet distribution is a conjugate prior of multinomial proportions, which is an extension of a β distribution (Johnson and Kotz 1969; Lee 1989):
The posterior distribution is obtained by multiplying the likelihood function, which is multinomial distribution in this case, by the prior. The posterior distribution of p is then given by
Empirical Bayes procedure: The primary disadvantage of using a Bayesian analysis for allele frequency estimation is that there is no obvious way of selecting a reasonable prior (Lange 1995). The Dirichlet distribution with α_{1} = · · · = α_{k} = ½ is a noninformative prior (Box and Tiao 1992). Here we adopt an empirical Bayes procedure to avoid dependence to priors. This procedure estimates the hyperparameters α by maximizing the marginallikelihood function (Maritz and Lwin 1989),
Lange (1995) estimated the hyperparameters from singlelocus data using Newton’s method. Here we estimate them from multilocus data by maximizing Equation 3 using a simplex minimization for the negative logarithm of Equation 3. Assuming that α is the same for H populations (samples) to be compared and Rα_{i} is also the same for J loci, the joint marginal likelihood is then given by
Parameter σ^{2} is the dispersion parameter that defines the magnitude of overdispersion; i.e., the variance of the response Y exceeds the nominal variance (McCullagh and Nelder 1983). For example, the expectation of the binomial random variables of a sample size m is E[Y] = mp and the variance is V[Y] = mp(1  p). If there is overdispersion, the variance is V[Y] σ^{2}mp(1  p) though the expectation remains the same, where Y has a density of a βbinomial distribution. For a multinomial event with overdispersion, the variancecovariance matrix of Y is σ^{2} times larger than that of the multinomial distribution, where Y has a density of a Dirichletmultinomial distribution.
Johnson and Kotz (1969) showed that the variancecovariance matrix of a Dirichletmultinomial distribution is
The binomial and multinomial counts are assumed to be taken by a simple random sampling, so the dispersion parameter σ^{2} is considered to indicate the magnitude of overshooting from a simple random sample. McCullagh and Nelder (1983) stated that “The simplest and perhaps the most common mechanism of overdispersion is clustering in the populations.” Kitada et al. (1994) estimated the dispersion parameter for fish tag recovery data and showed that the variances of the estimated mortality rates were
Standardized genetic distance: When allele frequencies at J loci are obtained from genetically identical populations,
Here we standardize D and propose a general index for the genetic distance. Performing a square root transformation to make the variance independent of the mean (Snedecor and Cochran 1967) and subtracting the expected value of I (the derivation is given in the appendix), we obtain a standardized genetic distance as
The χ^{2} distribution of
Effective population size: The effective population size is estimated from the temporal variation of allele frequencies in a population. Since the observed variance of the allele frequencies includes the sampling variance in addition to the genetic drift, we subtract the sampling variance when estimating the effective population size. Let us assume that we have two samples with sizes n_{0} and n_{t} from the population at generations 0 and t, respectively. The empirical Bayes procedure developed here can be extended to obtain the posterior distributions of the effective population size N_{e} by using the posterior distribution of Fstatistics calculated from the posterior distribution of allele frequencies.
The standardized variance of allele frequency change measured by Fstatistics has been used to estimate N_{e} (Krimbas and Tsakas 1971; Nei and Tajima 1981; Pollak 1983; Waples 1989). Among Fstatistics, F_{k} proposed by Pollak (1983) is similar in form to D and is given by
Equation 9 assumes that 2n_{0} and 2n_{t} genes are taken by a binomial sampling from the population. If there is overdispersion, F_{k} is overestimated, which leads to underestimation of N_{e}. Since the effective sample size is obtained by discounting the apparent sample size by dispersion parameters, Equation 9 is modified as
CASE STUDIES
Red sea bream: To evaluate genetic distances, we first analyzed the data of four populations of red sea bream (Pagrus major) from Tabata and Mizuta (1997; Table 1). From the fragment pattern of mtDNA Dloop regions with six restriction enzymes, 48 haplotypes were obtained for four wild populations. We decomposed the haplotype frequency to six allelic frequencies for each restriction enzyme and eliminated HaeIII and Msp, which showed little or no polymorphism, from the analysis.
The estimate of the total hyperparameters was 106.553 for each locus (Table 1), and the dispersion parameter was estimated at 1.80 by maximizing Equation 3. Here 2n. = (72 + 95 + 93 + 90)/4 = 87.5 because mtDNA is a haploid. With a prior distribution specified by these parameters, we obtained the posterior distribution of D by dividing RD_{j} over four loci by the number of loci (Table 2). As an example, the histogram and estimated density function of the posterior distribution of D at HinfI between Tanabe Bay and Tomogashima Channel is shown in Figure 1. D_{12} in Figure 2 was obtained as the mean of such four posterior genetic distances at the four loci. It should be noted that the posterior distances were overestimated including the overdispersion and the posterior distributions in Figure 2 were then overestimated. The means and SDs of D_{12}, D_{13}, and D_{14} were about two times larger than D_{23}, D_{24}, and D_{34}; however, they might include the effect of the smaller sample size of population 1 (Tanabe Bay).
The posterior distribution of the standardized genetic distance took the overdispersion and sample size difference into account. The means of I_{12}, I_{13}, and I_{14} ranged from 0.2706 to 0.4671, whereas those of I_{23}, I_{24}, and I_{34} took negative values. The SDs ranged from 0.53 to 0.58 and took almost the same values. The genetic differences with population 1(I_{12}, I_{13}, and I_{14}) looked larger than the others (Table 3). However, the posterior distributions of I overlapped well with the theoretical distribution of no genetic difference (Figure 3).
We estimated the 95% confidence interval of the dispersion parameter to be from 1.72 to 1.88 by the likelihoodratio test. The lower limit of the dispersion parameter corresponds to the upper limit of the genetic distance, from which we evaluate the difference. The means of the posterior distributions for the lower limit of the dispersion parameter were increased from 8 to 27% and SDs remained the same (Table 3), but the posterior distributions of I were almost the same as those for the point estimate of the overdispersion and still overlapped well with the theoretical distribution (Figure 3).
The value of 95% upper limit of the credibility region of the theoretical normal distribution of I with mean of 0 and variance of 0.5 is 1.16. All posterior means were <1.16, and the credibility regions included 0; hence we concluded that there was no genetic difference between the four populations of red sea bream. This finding agreed with the result of the original authors, who reported that the Roff test did not reject the homogeneity of the haplotype frequencies (Tabata and Mizuta 1997, p = 0.219). Nevertheless, they rejected the hypothesis by the pairwise comparison. From the results of our test, however, we argue that it was inappropriate analysis.
Herring: Stock enhancement of herring (Clupea pallasii) has been performed in Akkeshi Bay, Hokkaido (Japan). Because the matured herrings migrate to Akkeshi Bay to spawn, they are considered to have originated from Lake Akkeshi and Akkeshi Bay. Although wild adult fish that migrated to the bay are used for artificial spawning to produce juveniles every year, it still may be important to monitor the genetic change and estimate the effective population size to maintain the wild stock.
Temporal changes in allozyme allele frequencies were obtained by combining two studies on the same loci by Ando and Ohkubo (1997) and Hotta et al. (1999; Table 4). Independent samples were taken in March and April 1993. In 1996, males and females were taken separately from the sample, hence they were not independent. For the purposes of our case study, we treated the data as if they were taken independently.
The estimate of the total hyperparameters was 130.956 for each locus (Table 4), and the dispersion parameter was estimated at 2.39, with 2n. = (206 + 214 + 168 + 148)/4 = 184. The posterior distribution of F_{k} estimate was calculated from each of two sets of posterior distributions of allele frequencies for 1993 and 1996 by
F_{k} ranged from 0.0014 to 0.0814 with the mean and SD of 0.0226 and 0.0105, respectively. The posterior distribution of F_{k} is shown in the left side of Figure 4.
Most of the matured herring migrating to Akkeshi Bay to spawn are in their second year of life; the remainder are in their third year. The age composition of the spawners was surveyed and estimated at 0.9 and 0.1 for each age class by the Japan SeaFarming Association. The expected number of generations can be used for t in the estimating equations of N_{e} because the expectation of the Fstatistics was approximated to be linear with t as shown in Waples (1989, p. 382). We estimated the expected number of generations at 1.48 (= 3/2.1) divided by the number of years between samples by the mean age of spawners as Miller and Kapuscinski (1997) did. The samples were taken before reproduction, so Equation 10 was used by eliminating the term of 1/N and substituting (206 + 214)/2 = 210 for 2n_{0} and (168 + 148)/2 = 158 for 2n_{t}. The posterior means of N_{e} estimates and 95% credibility region of N_{e} are given in Table 5.
The dispersion parameter was estimated at 2.39 with a 95% confidence interval from 1.00 to 7.51. From a conservation viewpoint, it is better to consider the lower limit of N_{e}. The lower limit of the dispersion parameter evaluates the upper limit of F_{k} corresponding to the lower limit of N_{e}. The lower limit of σ^{2} was 1.00. Corresponding with that, no overdispersion arose and no subpopulation existed. The number of simulations with a negative value of N_{e} estimate was 1221 in 10,000 trials. When F_{k} ≤ [1/(2n_{0})  1/(2n_{t})], the only feasible estimate of N_{e} is infinity (Waples 1989). The mean of the positive N_{e} estimate was 350, and the 95% credibility region was estimated from 20 to infinity (Table 5). The posterior distribution of the positive N_{e} estimate is shown in the right side of Figure 4; it suggests the order of the effective population size of herring, even though the upper limit was not estimated.
Northern pike: We analyzed the data from Miller and Kapuscinski (1997) for comparison. Temporal changes in microsatellite allele frequency at seven loci were typed in the northern pike (Esox lucius) population of Lake Escanaba, Wisconsin. Of the seven loci, five had two alleles and two had three alleles. Following the data processing techniques of Williamson and Slatkin (1999), we also combined the two least common allelic classes for the two loci with three alleles and created a diallelic frequency set. To estimate the hyperparameters, we allocated sample sizes of 86 for 1961, 110 for 1977, and 72 for 1993 according to the diallelic frequencies for each locus to obtain the number of individuals corresponding to the frequencies. Because the data had one sample for each year, it was not possible to estimate the hyperparameters for each locus. Therefore, we assumed that the seven loci were independent of each other and had common hyperparameters. The two common hyperparameters and the dispersion parameter were estimated at 16.232, 3.089, and 9.74, with 2n. = (172 + 220 + 144)/3 = 178.7.
F_{k} between the year of 1977 and 1993 ranged from 0.0087 to 0.1264 with the mean and SD being 0.0479 and 0.0150, respectively. The posterior distribution of F_{k} is given in the left side of Figure 5. The posterior distribution of the N_{e} estimate was obtained by substituting the posterior distribution of F_{k} into Equation 10, eliminating the term 1/N, with t = 4, as given in Miller and Kapuscinski (1997). The mean of the positive N_{e} estimate was 73 and a 95% credibility region neglecting the overdispersion was estimated from 29 to 190 (Table 6). The number of negative N_{e} estimates was 6 in 10,000 trials. Miller and Kapuscinski’s estimates for 1977 and 1993 data were 0.038 for F_{k} and 72 for N_{e} with a 95% confidence interval from 17 to 258. Our posterior mean of F_{k} was 1.26 times larger and that of the N_{e} estimate coincided with a 67% narrower confidence interval (Table 6).
The 95% confidence interval of the dispersion parameter was estimated from 5.51 to 18.80. For the 95% lower limit of the dispersion parameter, the number of simulations with a negative value of N_{e} estimate was 8480 in 10,000 trials. The mean of the positive N_{e} estimates was 1065 and the 95% credibility region was estimated from 123 to infinity (Table 6). The posterior distribution of the positive N_{e} estimate, neglecting the overdispersion (σ^{2} = 1), and for the lower limit of
Ayu broodstock: Ayu (Plecoglossus altivelis) is the most popular target species of recreational anglers in rivers and streams in Japan. A total of 300 million juveniles are released every year, of which hatcheryproduced fish comprise ∼30%. The life span of ayu is 1 year. They spawn in a river from September to November and die after spawning. Hatched larvae go down to the sea and winter there. The upstream run of wild ayu juveniles begins from the coast in late March to early April and is over by early July. Soon after, they mature, spawn from September to November, and then die after spawning.
Hatcheries have commonly cultured broodstocks over generations. At the Gunma Prefecture Fisheries Experimental Station, adult ayu have been cultured over 27 generations. About 30004000 fish have been reared every year as broodstock, some of which are used for artificial fertilization. In 1996, ∼850 females and 650 males were used. The temporal changes in allozyme allele frequencies of the ayu broodstock were reported by Yoshizawa (1997; Table 7). These samples in Table 7 were taken after artificial fertilization from two rearing tanks in which males and females were kept separately.
The total of the hyperparameters was estimated at 28.576 for each locus (Table 7), and the dispersion parameter was estimated at 5.86, with 2n. = (190 + 210 + 128 + 130 + 100 + 100)/6 = 143. The 95% confidence interval of the dispersion parameter was estimated to range from 3.04 to 13.49. We calculated four F_{k}’s on the basis of temporal changes in allele frequencies observed in the first (19961997; F_{1}) and second (19971998) time intervals (F_{2}), over the entire interval (19961998; F_{3}), and for the entire interval based on the pooled F for the first two intervals, as Miller and Kapuscinski (1997) did. For the last case, Miller and Kapuscinski (1997) used the sum of F_{1} and F_{2} for the entire interval, but we used the mean for the two intervals (F_{mean}), which gives the same form of N_{e} estimate by Equation 10, substituting n_{0} by the harmonic mean of the sample size of the first and second year, and n_{t} by that of the second and third year. When t is equal for the two intervals, the estimate of N_{e} derived from pooled F_{k} is the harmonic mean for both sampling periods (Nei and Tajima 1981; Pollak 1983; Waples 1989; Miller and Kapuscinski 1997; Williamson and Slatkin 1999).
The values of F_{k} calculated by Equation 11 were varied using four combinations for two samples in each sampling year, reflecting the variation in allele frequencies for each sampling period. The posterior mean and SD of F_{3} were the largest, and the SD of F_{mean} was the smallest but almost the same as that of F_{1} (Table 8). The posterior distributions of F_{k} are illustrated in the left side of Figure 6.
We fixed the generation time at t = 1.0 for the first and second time intervals and had t = 2.0 for the entire interval because the life span of ayu is 1 year. The samples were taken after reproduction, so we used Equation 10, substituting 2n_{0} = (190 + 210)/2 = 200, 2n_{t} = (100 + 100)/2 = 100, and N = 1500, which was the total number of individuals used for artificial fertilization.
We made four estimates of N_{e} on the basis of F_{1}, F_{2}, F_{3}, and F_{mean}. Estimates for the 95% lower limit of the dispersion parameter (= 3.04) are given in Table 8. The posterior mean obtained using F_{3} was the largest with the largest SD, and that for F_{mean} was second with a smaller SD. The N_{e} estimate that took ∞ in 10,000 simulations was 8677 when using F_{mean}, and that for F_{3} was 4927. This is because the sampling correction term in the denominator of Equation 10 took the very similar values of 0.0912 for F_{3} and 0.0928 for F_{mean}, despite the posterior mean of F_{mean} being smaller than that of F_{3}, as shown in the left side of Figure 6. The posterior distributions of the positive N_{e} estimate obtained by using F_{1}, F_{2}, F_{3}, and F_{mean} for the lower limit of
We failed to estimate the upper limit of the credibility regions because of large sampling variances with overdispersion. However, when the numbers of breeding males N_{m} and females N_{f} are given, which is difficult to know in a wild population but possible in hatcheries, the effective population size is obtained by N_{e} = 4N_{m}N_{f}/(N_{m} + N_{f}) (Wright 1931). We obtained a value for N_{e} of 1473 by using the equation (4 × 850 × 650/(850 + 650)), which referred to the effective population size where a random mating was performed by artificial fertilization. In a spawning season of ayu, males and females eligible for spawning were selected every day from the broodstock and used for artificial fertilization. The number of females used in an artificial fertilization ranged from ∼10 to 20, and the ratio of males to females was ∼0.8. The eggs were squeezed from the females and stocked in a stainless bowl and then fertilized by squeezing sperm from individual males. This method of fertilization might not guarantee a random mating of the males and females used; hence, 1473 should be used as the upper limit of the credibility regions instead of ∞ (Table 8). If we neglect overdispersion, the 90% credibility region could be obtained at [13589] with the posterior mean of 136, which was underestimated.
DISCUSSION
The empirical Bayes procedure developed here makes it possible to describe the skewness of the genetic distance and evaluate genetic differentiation between populations while taking full account of the uncertainty of the sample allele frequencies. When we compare populations in which the genetic differentiation is small, the hypothesistesting framework cannot accept the null hypothesis of no genetic differentiation in almost all cases, because of the poor statistical power with relatively small sample sizes. The empirical Bayes procedure is effective even in such cases. So we believe it could play an important role in the field of conservation.
This general method can easily be extended to any parameter that is a function of multinomial frequencies. When the parameter of interest is a function of allele frequencies, the posterior distribution of that parameter can be obtained through the function by using the posterior distribution of the allele frequencies, instead of assuming a prior distribution for the parameter.
Overdispersion and empirical Bayes: Until now, models based on a simple random sampling from the gamete pool have been assumed when evaluating allele frequencies. However, as shown in the four case studies treated in this article, a simple random sampling is not necessarily guaranteed. If there are subpopulations divided spatially in a survey area, or a cluster of a genotype is taken, a sample allele frequency might be overdispersed. If overdispersion arises, a sampling variance becomes σ^{2} times larger than that for a simple random sampling. This can seriously affect the precision of the estimate of genetic distance and the effective population size. As a result, the genetic distance and Fstatistics can be overestimated, and the effective population size can be underestimated, if overdispersion is not taken into account in the analysis. Therefore, it is quite important to take overdispersion into account when estimating genetic distance and effective population size.
Sample sizes: If we use the noninformative prior of the Dirichlet distribution (Box and Tiao 1992), the dispersion parameter might be overestimated for most cases. For example, the dispersion parameter of herring was estimated at 92.5 from the noninformative prior, which was estimated at 2.39 from the empirical Bayes procedure, illustrating the importance of estimating the hyperparameters from the data.
We examined the relationship between sample size and the estimates of the hyperparameters using the data of red sea bream given in Table 1. We estimated the hyperparameters with multipliers of 0.5, 2, 3, and 4 to test each population with the same sample allele frequencies. The estimates of the hyperparameters were stable and not dependent upon sample sizes. This confirmed the robustness of the empirical Bayes procedure (Table 9). However, the dispersion parameter became larger as sample size increased. This is to be expected from the relationship between the total of the hyperparameters, the sample sizes, and the dispersion parameter given by Equation 5.
Suppose there are several subpopulations and the population allele frequencies are largely varied spatially. If the sample sizes are small, one might consider that the large variation in the sample allele frequencies is a function of the small sample size. On the other hand, if the same allele frequencies are obtained for larger sample sizes, one could consider that the large variation comes from the subpopulation structure with confidence. The more samples one draws, the more precisely one can estimate the dispersion parameter. In addition, increasing the number of polymorphic loci to be surveyed may also increase the information available for estimating the dispersion parameter, e.g., the precision of the dispersion parameter estimate in the case of red sea bream, in which the narrowest confidence interval was obtained among our four case studies.
It is also quite important to consider sampling strategies to minimize overdispersion caused by sampling procedures. For example, sampling from different sites and times may be useful to avoid sampling clusters of individuals having the same genotype. Such multiple samples contribute simultaneously to a more precise estimation of the dispersion parameter.
Standardized genetic distance: As I follows the normal distribution, it takes values between ∞ and +∞. For simplicity, let
Overdispersion and gene flow: Weir (1996) stated that for populations that have reached an equilibrium under the joint effects of drift and mutation or migration, Wright (1945) found that allele frequencies for loci with two alleles had a β distribution, and for multiallele loci the distribution was Dirichlet (Wright 1951). We assumed that the hyperparameters for the β or Dirichlet distributions were common for every sample and locus that was in an overdispersed population. Our assumption corresponds with Wright’s (1945, 1951) theories. So if the random sampling is performed, the estimated hyperparameters and dispersion parameter both describe a kind of genetic differentiation between populations that have reached an equilibrium. If all populations mate randomly, the total variance of allele frequency p with two alleles of 2n genes is given by Weir (1996)
Rannala and Hartigan (1996) proposed the pseudomaximumlikelihood method (PMLE) for estimating the rate of gene flow into island populations using the distribution of alleles in samples from a number of islands. We confirmed that their likelihood function for multiple loci (p. 149 Equation 10) coincides with Equation 3 by using the relationship of Γ(n) = (n  1)!. In PMLE, α_{i} is treated by θp_{i}. Here, p_{i} is a nuisance parameter estimated from the data as pˆ_{i} = n_{i}./n.., and then pˆ_{i} is substituted for p_{i} in the loglikelihood function, and the only unknown parameter θ is estimated by using the Newton method. By contrast, we directly maximize the negative loglikelihood function and estimate
Wright (1969) proposed the estimator of θ for a discretegeneration island model of a population at equilibrium, based on F_{ST} as
The essential idea for estimating overdispersion is to compare the variation of sample allele frequencies obtained from the different locations to the multinomial variance. In addition, the effective population size is based on the changes in allele frequencies between generations. Conversely, overdispersion provides insight into the spatial variation of allele frequencies. By evaluating the spatial variation, it might become possible to discriminate the overdispersion resulting from the variation between generations. Hence, the procedure needs to evaluate overdispersion as a function of the spatial variation and then measure the variation between generations taking overdispersion into account.
In the three case studies we looked at for estimating the effective population size, direct information on the spatial variation was scarce. Therefore, the precision of the dispersion parameter was marginal. When subpopulations exist, overdispersion arises and affects the estimation of the effective population size. It is then important to collect data on the spatial variation. At the same time, when many isolated subpopulations exist, the effective population size is considered to be close to the size of a subpopulation. When this occurs, it seems dangerous to dismiss the variation between generations as overdispersion. It needs further consideration.
Practical considerations on estimating N_{e}: From the approximate variance formula of N_{e} estimate (Pollak 1983, Equations 28 and 29; Waples 1989, Equation 17), it is clear that increasing the sample size, the number of loci, and the number of generations t simultaneously ensures greater precision for the estimate of N_{e} (Waples 1989). Miller and Kapuscinski (1997) stated that if N_{e} is expected to be moderately large, the sample size, the number of loci, and the number of generations should all be as large as possible. To improve the precision of the estimate of N_{e}, it is essential to reduce the sampling variance and increase information on genetic drift.
Sample size: The idea of the temporal method is to estimate N_{e} from the genetic change over time described by Fstatistics estimated from the sample allele frequencies. Fstatistics, then, consist of the genetic drift and the sampling variance. To evaluate the genetic drift, we have to subtract the sampling variances from the Fstatistics. The second and third terms in the denominator of Equation 10 are the sampling variances at generations 0 and t. If N_{e} is large, the genetic drift may be small, so the denominator of Equation 10 would take a negative value, which leads to an infinite N_{e} for small sample sizes n_{0} and n_{t}. If overdispersion arises, the effect of subtracting the sampling variances becomes σ^{2} times larger, which is why we failed to estimate the upper limit of the credibility region of N_{e}. As pointed out by Waples (1989), the temporal method should be useful for cases of small N_{e}, where larger genetic drift is expected. Even in the case of a small N_{e}, the problem of an infinite N_{e} estimate can occur due to large sampling variance, as shown in the ayu studies, because of the small sample sizes. When one uses the temporal method, reducing the sampling variance is indispensable. The sample size should be kept as large as possible. A larger sample size also provides greater information on the genetic drift.
Number of loci: Williamson and Slatkin (1999) developed a maximumlikelihood temporal method to estimate N_{e} and compared estimates with those derived with the Fstatistic method. The simulation result in their Table 1 showed that increasing the number of loci reduced the variance and bias in both estimators, although when the number of loci was >50, the corresponding reduction of variance and bias was not large, and the total information on allele frequency changes did not increase much. The results of Williamson and Slatkin (1999) suggest that information on genetic drift was not improved much even if the number of loci was >100, because their simulation was based on diallelic alleles. Fstatistics measure a magnitude of changes in allele frequencies per allele, which can be regarded as a sample mean. So, the estimation precision can be improved if the number of loci is increased. This suggests that increasing the number of alleles is essential, which can be attained by increasing the number of loci.
Number of years between samples: The number of years between samples is correlated with the number of generations, and it then affects the precision of the estimate of N_{e}. A large number of generations between samples can improve the precision of the estimate of N_{e} (Waples 1991), because information on genetic drift increases as the number of generations increases. Williamson and Slatkin (1999, Table 1) showed through simulated populations sampled at generations 04 and 08 that the variance and bias in both estimators were reduced when the number of years between samples was doubled, although the effect of reducing the bias was not clearly observed with the Fstatistic method.
For the case study of ayu, the posterior mean of N_{e} was 350 based on F_{1} and 796 based on F_{3}, and SDs for the two estimates were 4024 and 18,538, respectively, showing the reduction of precision despite the fact that the number of years between samples was doubled (Table 8). This is because the doubled number of generations increased the variance of allele frequency changes. The numbers of infinite N_{e} estimates in 10,000 simulations were 8453 on F_{1} and 4927 on F_{3}, and the smaller F values increased the estimated value of N_{e}. The result was similar for northern pike. The point estimate of N_{e} based on F_{3} (= 125) was larger than those based on F_{1} (= 35) and F_{2} (= 72), and the confidence interval for F_{1} was the largest (Miller and Kapuscinski 1997).
Miller and Kapuscinski (1997) discussed the question of which estimate obtained from F_{3} and F_{mean} to use for the entire time interval. If N_{e} changes between the two sampling intervals, it should be evaluated by using F_{1} and F_{2} for the respective intervals. The numbers of adult ayu used for the artificial fertilization in 1997 were 600 females and 480 males, which are lower numbers than those used in 1996. The posterior mean of F_{2} was larger than that of F_{1}, which resulted in a smaller N_{e} estimate based on F_{2}, supporting the fact that N_{e} actually changed (Table 8 and Figure 6). One can use the estimate of N_{e} based on F_{mean} as the harmonic mean of the effective population size in the respective intervals. The precision was improved by using F_{mean}, in which a larger quantity of information was included. The greater precision that occurred with using F_{mean} and the lower precision that occurred with using F_{3} for the case study of ayu are shown clearly in Figure 6.
When N_{e} does not change in the entire interval, F_{3} is expected to have more information on genetic drift than F_{1} and F_{2}. Williamson and Slatkin (1999, Table 1) also showed by their simulation that the estimates of N_{e} based on F_{mean} had smaller variances and biases than those based on F_{1} and F_{3}. This suggests that the decision about which estimate obtained from F_{3} and F_{mean} to use should be made on the basis of the relative effect of improving precision by using F_{mean} and a doubled number of years. Which estimate has more information on the genetic drift must be determined on a casebycase basis.
Overlapping generations: The basic theory of the temporal method assumes generations to be discrete. The expected number of generations used in Equation 10 directly affects the estimate of N_{e}. We take time to be measured in years. The expected number of generations between samples can be estimated by dividing the number of years between samples by the mean generation time, which corresponds to the mean age of maturity. In the case of ayu, since the life span is 1 year, 1 year coincides with one generation, which makes it possible to estimate E[t] by the abovementioned method.
In the case of herring and salmon, however, where there are overlapping year classes of spawners, the estimation of E[t] is complex. When generations overlap, the agespecific birth rate may essentially affect the estimate of E[t]. Hill (1979) showed that estimates of N_{e} are robust with overlapping generations if demographic parameters are stable. If demographic parameters change over time, Fˆ may be biased upward, leading to an estimate of N_{e} that is too small (Pollak 1983; Waples 1989). Jorde and Ryman (1995) proposed a correction method for the bias and showed by using simulations for short time intervals that the bias was larger for a case where agespecific birth rates were extremely different compared with a case with equal agespecific birth rates. Waples (1990) developed a statistical method for this situation that can be applied to Pacific salmon populations that have an unusual life history of semelparity with overlapping year classes. Tajima (1992) developed a formula to estimate the expected number of generations without computer simulations and showed the estimates agreed well with estimates obtained by the method of Waples (1990), which requires computer simulations.
In the case of herring, agespecific survival and birth rates were unknown, so it was not possible to apply the method of Jorde and Ryman (1995), which requires these demographic parameters. If an individual continues to spawn every year after the first spawning, like herring, E[t] may be estimated by dividing the years between samples by the mean age of spawners, leading to the estimate of E[t] at 1.48 (= 3/2.1). If a distribution of agespecific birth rates is concentrated to a specific age, E[t] may be close to the estimate obtained by the methods of Waples (1990) or Tajima (1992). We estimated the number of steps at 12 and the expected number of generations at 5.71 (= 12/2.1) for a time interval of 3 yr between samples by using the computer program given in Tajima (1992), where we substituted f(2) = 0.9, f(3) = 0.1, and f(4) = 0. The downward bias of Nˆ_{e} when demographic parameters change over time with overlapping generations should be corrected upward. From a conservation viewpoint, the estimate of N_{e} without the correction must be conservative for the overlapping generations.
Rate of inbreeding: As another evaluation of breeding population size, the inbreeding coefficient may be useful, especially for cases where the population size is estimable, as it is in the field of fishery science. Crow (1954) used the inbreeding effective size, making a distinction between that and the variance effective size, which was defined by an inverse of the inbreeding coefficient. However, it is known that there is no great difference between the two effective sizes (Nei 1987); hence we calculated the posterior distribution of the rate of inbreeding defined as 1/(2N_{e}) (Falconer and Mackay 1996), as an infinite N_{e} corresponds to an inbreeding coefficient of 0, obtained from the 95% credibility region. The means, SDs, and 95% credibility regions of the posterior distributions of the rate of inbreeding were 0.0083, 0.0070, and [0, 0.0251] for herring; 0.0004, 0.0012, and [0, 0.0041] for northern pike; and 0.0011, 0.0041, and [0, 0.0140] for ayu.
The posterior mean of herring was 7.5 times larger than that of ayu with a righttailed credibility region shown in Figure 7. Overdispersion for herring was underestimated, because the samples for males and females taken in 1996 were analyzed as independent samples even though they were the same sample, leading to a smaller value of N_{e}, which caused the rate of inbreeding to be overestimated. The mean for northern pike was smallest with the narrowest credibility region. However, these values may be underestimated because of an overestimated dispersion parameter of northern pike, which was the largest among our four case studies. There was only one sample for one sampling year, and we assumed that the seven loci had common hyperparameters, so the estimated dispersion parameter may include the change of the allele frequencies.
Multistage sampling in hatcheries: All existing methods assume that N_{e} is drawn from a gamete pool by a simple random sampling. This is an appropriate assumption for the reproduction of a wild population. However, for broodstocks cultured over generations in hatcheries, candidates of the next broodstock are sampled from the progenies produced by the broodstock. Therefore, N_{e} is drawn from the progenies by a twostage sampling. If artificial fertilization using a part of the candidates is performed, as in the case study of ayu, N_{e} is drawn from the progenies by a threestage sampling and the sample is drawn from the candidates to estimate the allele frequencies, which is therefore a twostage sampling of the progenies. The multistage sampling must lead to the different form of V(x  y) given in Waples (1989). This is a problem that needs further research, but it should be noted that the variances corresponding to the twostage and threestage sampling become small when the sample sizes are large. In the case of ayu, a total of 30004000 candidates were sampled from the progenies and cultured in rearing tanks, and 1500 adult fish from the candidates were used for artificial fertilization. Hence the sample allele frequencies of ayu were expected to represent those of the progenies produced by the broodstock. However, if the sample sizes are small, V(x  y) is seriously affected.
APPENDIX
Expectation of
Acknowledgments
We thank ZhaoBang Zeng and two anonymous referees for their comments on an earlier version of this article. We also thank Ray Timm for critical review of the manuscript, Fumio Tajima for important suggestions made during our research, Kazutomo Yoshizawa for biological information on ayu broodstocks including unpublished data, and Masashi Yokota for helpful discussions.
Footnotes

Communicating editor: ZB. Zeng
 Received February 17, 2000.
 Accepted July 26, 2000.
 Copyright © 2000 by the Genetics Society of America