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 marginal-likelihood 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 marginal-likelihood 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 10-point 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 Ne in a population based on temporal changes in allele frequencies. The joint marginal-likelihood 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 p11, · · ·, p1k and p21, · · ·, p2k. 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 = (p1, · · ·, pk)′ and the sample count be n = (n1, · · ·, nk)′, where Rpi = 1 and Rni = 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 marginal-likelihood function (Maritz and Lwin 1989),
Lange (1995) estimated the hyperparameters from single-locus 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] σ2mp(1 - p) though the expectation remains the same, where Y has a density of a β-binomial distribution. For a multinomial event with overdispersion, the variance-covariance matrix of Y is σ2 times larger than that of the multinomial distribution, where Y has a density of a Dirichlet-multinomial distribution.
Johnson and Kotz (1969) showed that the variance-covariance matrix of a Dirichlet-multinomial 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 n0 and nt 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 Ne by using the posterior distribution of F-statistics calculated from the posterior distribution of allele frequencies.
The standardized variance of allele frequency change measured by F-statistics has been used to estimate Ne (Krimbas and Tsakas 1971; Nei and Tajima 1981; Pollak 1983; Waples 1989). Among F-statistics, Fk proposed by Pollak (1983) is similar in form to D and is given by
Equation 9 assumes that 2n0 and 2nt genes are taken by a binomial sampling from the population. If there is overdispersion, Fk is overestimated, which leads to underestimation of Ne. 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 D-loop 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.
Allele frequencies of the mtDNA D-loop region from Tabata and Mizuta (1997) and estimated hyperparameters for four populations of red sea bream from eastern Japan
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 RDj 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. D12 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 D12, D13, and D14 were about two times larger than D23, D24, and D34; 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 I12, I13, and I14 ranged from 0.2706 to 0.4671, whereas those of I23, I24, and I34 took negative values. The SDs ranged from 0.53 to 0.58 and took almost the same values. The genetic differences with population 1(I12, I13, and I14) 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 likelihood-ratio 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 Fk estimate was calculated from each of two sets of posterior distributions of allele frequencies for 1993 and 1996 by
Fk ranged from 0.0014 to 0.0814 with the mean and SD of 0.0226 and 0.0105, respectively. The posterior distribution of Fk 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 Sea-Farming Association. The expected number of generations can be used for t in the estimating equations of Ne because the expectation of the F-statistics 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 2n0 and (168 + 148)/2 = 158 for 2nt. The posterior means of Ne estimates and 95% credibility region of Ne are given in Table 5.
Means and SDs of the posterior distribution of the genetic distance for the red sea bream populations
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 Ne. The lower limit of the dispersion parameter evaluates the upper limit of Fk corresponding to the lower limit of Ne. 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 Ne estimate was 1221 in 10,000 trials. When Fk ≤ [1/(2n0) - 1/(2nt)], the only feasible estimate of Ne is infinity (Waples 1989). The mean of the positive Ne estimate was 350, and the 95% credibility region was estimated from 20 to infinity (Table 5). The posterior distribution of the positive Ne 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.
—A histogram and the estimated density function of the posterior distribution of the genetic distance between red sea bream populations of Tanabe Bay and Tomogashima Channel for HinfI using data given in Table 1. D12 in Figure 2 was obtained as the mean of the four posterior genetic distances for the four loci.
Fk 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 Fk is given in the left side of Figure 5. The posterior distribution of the Ne estimate was obtained by substituting the posterior distribution of Fk into Equation 10, eliminating the term 1/N, with t = 4, as given in Miller and Kapuscinski (1997). The mean of the positive Ne estimate was 73 and a 95% credibility region neglecting the overdispersion was estimated from 29 to 190 (Table 6). The number of negative Ne estimates was 6 in 10,000 trials. Miller and Kapuscinski’s estimates for 1977 and 1993 data were 0.038 for Fk and 72 for Ne with a 95% confidence interval from 17 to 258. Our posterior mean of Fk was 1.26 times larger and that of the Ne estimate coincided with a 67% narrower confidence interval (Table 6).
—Posterior distributions of the genetic distance (D) between four populations of red sea bream using data given in Table 1.
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 Ne estimate was 8480 in 10,000 trials. The mean of the positive Ne estimates was 1065 and the 95% credibility region was estimated from 123 to infinity (Table 6). The posterior distribution of the positive Ne 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 hatchery-produced 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.
Means and 95% credibility regions of the posterior distribution of the standardized genetic distance for the red sea bream populations
Hatcheries have commonly cultured broodstocks over generations. At the Gunma Prefecture Fisheries Experimental Station, adult ayu have been cultured over 27 generations. About 3000-4000 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.
—Posterior distributions of the standardized genetic distance (I) for the red sea bream populations taking the point estimate of the dispersion parameter 1.80 (left) and the 95% lower limit 1.72 (right) into account.
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 Fk’s on the basis of temporal changes in allele frequencies observed in the first (1996-1997; F1) and second (1997-1998) time intervals (F2), over the entire interval (1996-1998; F3), 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 F1 and F2 for the entire interval, but we used the mean for the two intervals (Fmean), which gives the same form of Ne estimate by Equation 10, substituting n0 by the harmonic mean of the sample size of the first and second year, and nt by that of the second and third year. When t is equal for the two intervals, the estimate of Ne derived from pooled Fk is the harmonic mean for both sampling periods (Nei and Tajima 1981; Pollak 1983; Waples 1989; Miller and Kapuscinski 1997; Williamson and Slatkin 1999).
Temporal changes in allozyme allele frequencies and estimated hyperparameters of herring in Akkeshi Bay
—Posterior distributions of the standardized variance of allele frequency change Fk and effective population size Ne of herring for the 95% lower limit of the dispersion parameter (σ2 = 1) using data in Table 4.
The values of Fk 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 F3 were the largest, and the SD of Fmean was the smallest but almost the same as that of F1 (Table 8). The posterior distributions of Fk 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 2n0 = (190 + 210)/2 = 200, 2nt = (100 + 100)/2 = 100, and N = 1500, which was the total number of individuals used for artificial fertilization.
Means and 95% credibility regions of the posterior distribution of the effective population size of herring in Akkeshi Bay obtained using data in Table 4
We made four estimates of Ne on the basis of F1, F2, F3, and Fmean. Estimates for the 95% lower limit of the dispersion parameter (= 3.04) are given in Table 8. The posterior mean obtained using F3 was the largest with the largest SD, and that for Fmean was second with a smaller SD. The Ne estimate that took ∞ in 10,000 simulations was 8677 when using Fmean, and that for F3 was 4927. This is because the sampling correction term in the denominator of Equation 10 took the very similar values of 0.0912 for F3 and 0.0928 for Fmean, despite the posterior mean of Fmean being smaller than that of F3, as shown in the left side of Figure 6. The posterior distributions of the positive Ne estimate obtained by using F1, F2, F3, and Fmean for the lower limit of
—Posterior distributions of the standardized variance of allele frequency change Fk and effective population size Ne of northern pike for the 95% lower limit of the dispersion parameter (σ2 = 5.51) and that with no overdispersion (σ2 = 1.00) using 1977-1993 data from Miller and Kapuscinski (1997).
Means and 95% credibility regions of the posterior distribution of the effective population size of northern pike obtained using 1977-1993 data from Miller and Kapuscinski (1997)
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 Nm and females Nf are given, which is difficult to know in a wild population but possible in hatcheries, the effective population size is obtained by Ne = 4NmNf/(Nm + Nf) (Wright 1931). We obtained a value for Ne 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 [13-589] with the posterior mean of 136, which was underestimated.
Temporal changes in allozyme allele frequencies from Yoshizawa (1997) and estimated hyperparameters of ayu broodstock cultured over 27 generations in Gunma Prefecture
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 hypothesis-testing 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.
Means and 95% credibility regions of the posterior distribution of the effective population size of ayu for 95% lower limit of the dispersion parameter (3.04) obtained using data in Table 7
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 F-statistics 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.
—Posterior distributions for the four estimators of the standardized variance of allele frequency change Fk and effective population size Ne of ayu broodstock using 1996-1998 data in Table 7.
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
Estimated hyperparameters and the dispersion parameter for four populations of red sea bream for sample sizes of 0.5, 2, 3, and 4 times larger than the original one with the same sample allele frequencies
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 pseudomaximum-likelihood 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 θpi. Here, pi is a nuisance parameter estimated from the data as pˆi = ni./n.., and then pˆi is substituted for pi in the log-likelihood function, and the only unknown parameter θ is estimated by using the Newton method. By contrast, we directly maximize the negative log-likelihood function and estimate
Estimated hyperparameters and the dispersion parameter from the mtDNA haplotype distribution among islands for Channel Island foxes (Table 2 of Rannala and Hartigan 1996), using the full-likelihood function (Equation 3)
Wright (1969) proposed the estimator of θ for a discrete-generation island model of a population at equilibrium, based on FST 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 Ne: From the approximate variance formula of Ne 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 Ne (Waples 1989). Miller and Kapuscinski (1997) stated that if Ne 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 Ne, 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 Ne from the genetic change over time described by F-statistics estimated from the sample allele frequencies. F-statistics, then, consist of the genetic drift and the sampling variance. To evaluate the genetic drift, we have to subtract the sampling variances from the F-statistics. The second and third terms in the denominator of Equation 10 are the sampling variances at generations 0 and t. If Ne is large, the genetic drift may be small, so the denominator of Equation 10 would take a negative value, which leads to an infinite Ne for small sample sizes n0 and nt. 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 Ne. As pointed out by Waples (1989), the temporal method should be useful for cases of small Ne, where larger genetic drift is expected. Even in the case of a small Ne, the problem of an infinite Ne 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 maximum-likelihood temporal method to estimate Ne and compared estimates with those derived with the F-statistic 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. F-statistics 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 Ne. A large number of generations between samples can improve the precision of the estimate of Ne (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 0-4 and 0-8 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 F-statistic method.
For the case study of ayu, the posterior mean of Ne was 350 based on F1 and 796 based on F3, 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 Ne estimates in 10,000 simulations were 8453 on F1 and 4927 on F3, and the smaller F values increased the estimated value of Ne. The result was similar for northern pike. The point estimate of Ne based on F3 (= 125) was larger than those based on F1 (= 35) and F2 (= 72), and the confidence interval for F1 was the largest (Miller and Kapuscinski 1997).
Miller and Kapuscinski (1997) discussed the question of which estimate obtained from F3 and Fmean to use for the entire time interval. If Ne changes between the two sampling intervals, it should be evaluated by using F1 and F2 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 F2 was larger than that of F1, which resulted in a smaller Ne estimate based on F2, supporting the fact that Ne actually changed (Table 8 and Figure 6). One can use the estimate of Ne based on Fmean as the harmonic mean of the effective population size in the respective intervals. The precision was improved by using Fmean, in which a larger quantity of information was included. The greater precision that occurred with using Fmean and the lower precision that occurred with using F3 for the case study of ayu are shown clearly in Figure 6.
When Ne does not change in the entire interval, F3 is expected to have more information on genetic drift than F1 and F2. Williamson and Slatkin (1999, Table 1) also showed by their simulation that the estimates of Ne based on Fmean had smaller variances and biases than those based on F1 and F3. This suggests that the decision about which estimate obtained from F3 and Fmean to use should be made on the basis of the relative effect of improving precision by using Fmean and a doubled number of years. Which estimate has more information on the genetic drift must be determined on a case-by-case 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 Ne. 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 above-mentioned 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 age-specific birth rate may essentially affect the estimate of E[t]. Hill (1979) showed that estimates of Ne 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 Ne 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 age-specific birth rates were extremely different compared with a case with equal age-specific 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, age-specific 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 age-specific 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 Ne without the correction must be conservative for the overlapping generations.
—Posterior distributions of the rate of inbreeding of herring, ayu broodstock, and northern pike for the 95% lower limit of the dispersion parameter.
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/(2Ne) (Falconer and Mackay 1996), as an infinite Ne 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 right-tailed 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 Ne, 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 Ne 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, Ne is drawn from the progenies by a two-stage sampling. If artificial fertilization using a part of the candidates is performed, as in the case study of ayu, Ne is drawn from the progenies by a three-stage sampling and the sample is drawn from the candidates to estimate the allele frequencies, which is therefore a two-stage 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 two-stage and three-stage sampling become small when the sample sizes are large. In the case of ayu, a total of 3000-4000 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 Zhao-Bang 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: Z-B. Zeng
- Received February 17, 2000.
- Accepted July 26, 2000.
- Copyright © 2000 by the Genetics Society of America
