Inferring Population History From Microsatellite and Enzyme Data in Serially Introduced Cane Toads, Bufo marinus

Corresponding author: Arnaud Estoup, Centre de Biologie et de Gestion des Populations, Campus International de Baillarguet, CS 30 016, 34980 Montferrier/Lez, France., estoup{at}ensam.inra.fr (E-mail)

Communicating editor: J. HEIN

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Much progress has been made on inferring population history from molecular data. However, complex demographic scenarios have been considered rarely or have proved intractable. The serial introduction of the South-Central American cane toad Bufo marinus in various Caribbean and Pacific islands involves four major phases: a possible genetic admixture during the first introduction, a bottleneck associated with founding, a transitory population boom, and finally, a demographic stabilization. A large amount of historical and demographic information is available for those introductions and can be combined profitably with molecular data. We used a Bayesian approach to combine this information with microsatellite (10 loci) and enzyme (22 loci) data and used a rejection algorithm to simultaneously estimate the demographic parameters describing the four major phases of the introduction history. The general historical trends supported by microsatellites and enzymes were similar. However, there was a stronger support for a larger bottleneck at introductions for microsatellites than enzymes and for a more balanced genetic admixture for enzymes than for microsatellites. Very little information was obtained from either marker about the transitory population boom observed after each introduction. Possible explanations for differences in resolution of demographic events and discrepancies between results obtained with microsatellites and enzymes were explored. Limits of our model and method for the analysis of nonequilibrium populations were discussed.

THE evolutionary history of natural or introduced populations typically involves complex changes in (effective) size and occasionally genetic admixture between more or less differentiated forms. This represents a major challenge to population genetic theory, most of which is based on equilibrium or simple nonequilibrium models. Much emphasis has been put on inferring historical processes from molecular data (reviewed in BEAUMONT 1999 ; STEPHENS and DONNELLY 2000 ). Recent progress in this field has been possible thanks to the development of the coalescent theory (KINGMAN 1982 ; HUDSON 1991 ), as well as the advent of a new category of genetic markers, microsatellites (reviewed in ESTOUP and ANGERS 1998 ). Under the coalescent framework, changes in population size and structure correspond to rescaling time in the coalescent tree (reviewed in DONNELLY and TAVARE 1995 ). Therefore patterns of coalescence provide indirect information about past demography. However, coalescent events inferred from data are usually imprecise. If the mutation rate is high, then any trace of the earliest coalescences is eliminated by subsequent mutations, whereas if the mutation rate is low there will be insufficient mutations to "document" many coalescence events. Hence, analyzing multiple categories of markers with different rates and forms of mutation (e.g., microsatellites and enzymes as in the present study) can be particularly fruitful. Imprecision also arises from the fact that genealogical models such as the coalescent model describe a complex stochastic process evolving through time, whereas most data sets, including those examined in this article, provide information at only a single time point. In this context, genetic data appear to be profitably combined with information from other sources such as paleontology, archeology, demography, and historical records, for instance, using a Bayesian approach (e.g., TAVARE et al. 1997 ; WILSON and BALDING 1998 ; PRITCHARD et al. 1999 ).

In recent years, important advances have been made toward the goal of estimating the probability of obtaining a given gene sample configuration to make "fully likelihood-based" statistical inference from molecular data, rather than drawing inferences based on summary statistics (BEAUMONT 1999 ; BAHLO and GRIFFITH 2000 ; STEPHENS and DONNELLY 2000 ). However, those approaches require an implementation of models and algorithms, which remains extremely challenging for populations with a complex evolutionary history, so that most applications have focused on relatively simple demographic models (e.g., WILSON and BALDING 1998 ; BEAUMONT 1999 ). For a complex evolutionary history, inferential methods that are not fully likelihood based still appear to be the best option available (reviewed in FU and LI 1999 ).

In this context, the biological and demographic model studied here presents four major points of interest: (i) The serial introduction of cane toad (Bufo marinus) in the Caribbean and Pacific islands exemplifies a complex historical process involving genetic admixture and a series of population size fluctuations; (ii) it is a recent process, which contrasts with the larger timescale of the demographic events previously studied, for instance, in human populations (e.g., WILSON and BALDING 1998 ; PRITCHARD et al. 1999 ); (iii) large molecular data sets can be compared for two different categories of genetic markers, microsatellites and enzymes; and (iv) a large amount of historical, geographical, and demographic information is available for those introductions (see MATERIALS AND METHODS). We used a Bayesian approach to combine this information with large molecular data sets and used a rejection algorithm to estimate simultaneously demographic parameters describing the four major phases of the introduction history. Historical inferences were compared for the two categories of studied markers.

	MATERIALS AND METHODS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Demographic model:
The cane toad is native to the American tropics and was deliberately introduced as a biocontrol agent between the middle of the 19th and 20th centuries in numerous islands of the Caribbean and the Pacific. A large amount of historical and demographic information about those introductions is available (reviewed in EASTEAL 1981 ). Here we focus on the information that was useful in formulating our model (Fig 1).

View larger version (22K): In this window In a new window Download PPT slide   Figure 1. Summary of the introduction history of cane toad in Caribbean and Pacific islands from two South American source populations and demographic parameters describing our model (italic letters). This starts from the initial release points in Australia. A similar scheme without the introduction events in Australia would summarize the introductions in the Hawaiian islands of Oahu, Kauai, Molokai, Maui, and Hawaii. Not shown to scale.

The early introductions were made from both French Guyana and British Guyana to Barbados in 1833, from Barbados to Jamaica in 1844, from Jamaica to the northeast of Puerto Rico in 1923, from the northeast of Puerto Rico to the island of Oahu (Hawaiian archipelagos) in 1932, and from Oahu to Gordonvale in Australia (northeast Queensland) in 1935. Forty individuals were introduced in northeast Puerto Rico, 148 in Oahu, and 101 in Gordonvale. The number of introduced individuals is unknown for Barbados and Jamaica. B. marinus is a highly prolific species (7500–20,000 eggs/female; ALFORD et al. 1995 ). This feature at least partly explains the observation that, following introductions, populations undergo a sudden and transitory increase in size preceding demographic stabilization (cf. the "boom-and-bust" phenomenon reviewed by WILLIAMSON 1996 for introduced species).

Information on the origin and dates of introductions allowed us to fix the times for the admixture and population size fluctuation events in the coalescent gene tree and hence to limit the number of parameters in our model (Fig 1). The population size fluctuations were specified by four variable parameters: the number of introduced individuals (N₃), the effective population size after the demographic boom (N₂) and its duration (D_B), and the effective population size after stabilization (N₁). The sudden character of the population explosion documented in all islands made it sensible to assume that the introduced population goes from N₃ to N₁ in a single generation and from N₁ to N₂ in a single generation too. Population reduction to N₁ also appears to be rapid (EASTEAL 1981 ; A. ESTOUP, personal observations) and hence modeled (for simulation convenience too) by a single generation transition from N₂ to N₁. We also assumed that the effective number of individuals introduced was the same in each island so that this part of our model was specified by a single variable parameter (N₃). Since the individuals introduced to Barbados from South America have two geographical origins, we included in our model the possibility of genetic admixture in Barbados. This evolutionary event was specified by two variable parameters, the admixture rate (A_R), defined as the proportion of individuals originating from each source population, and the time since separation between the two source populations (T_S).

We used genetic data obtained from populations introduced to Australia and the Hawaiian islands. In both cases the toads were reared and released initially in one area and the offspring were then distributed to five discrete areas: in Hawaii, the islands of Kauai, Molokai, Maui, Hawaii, and Oahu (from 1933 to 1935) and, in Australia, Bundaberg, MacKay, Ayr, Ingham, and Gordonvale districts (from 1936 to 1937). The number of individuals involved in these distributions was large, of the order of tens of thousands (reviewed in EASTEAL 1985 ). Our model (Fig 1) took these aspects into account by considering that a single (bottlenecked) generation of population size N₃ occurred prior to independent evolution of the five Australian or the five Hawaiian populations without further bottlenecks. Because of the large geographical distance between them or their island nature, the five Australian and Hawaiian introduction sites were considered to evolve as totally isolated demes until the sampling dates.

Population sampling and marker analysis:
The B. marinus populations sampled for this study correspond to the five main introduction points in Australia (Bundaberg, Mackay, Gordonvale, Ayr, and Ingham), the five main introduction points in the Hawaiian islands (Kauai, Molokai, Maui, Hawaii, and Aiea), and the two source populations (Kourou in French Guyana and Georgetown in British Guyana). The exact locations of the 10 introduced populations are given in EASTEAL 1985 . Samples were collected in 1980 and 1999 for the enzyme and microsatellite analyses, respectively.

In making inferences about population history it is desirable to use data from many unlinked loci (e.g., DONNELLY and TAVARE 1995 ). In this study, 10 microsatellite and 22 enzyme markers were analyzed. Previous population studies did not show evidence for statistical association between loci other than that expected under particular demographic scenarios (EASTEAL 1985 ; LEBLOIS et al. 2000 ; TIKEL et al. 2000 ). An ascertainment problem occurs frequently in analysis of molecular markers. Loci are commonly chosen because they were known to be polymorphic (i.e., variable) in the population studied, while monomorphic loci were not genotyped or discarded from analysis. Inferences that can be drawn under these circumstances will differ in general from inferences that would be drawn from loci genotyped "at random." This problem was minimized here as our 10 microsatellites were chosen randomly, and we included 12 monomorphic enzyme loci that were scored initially in the same populations and discarded from the data set published by EASTEAL 1985 .

The data set of 22 enzyme loci genotyped across the 10 points of introduction was produced by EASTEAL 1985 . We analyzed the 5 Australian introduction points for 10 microsatellite markers. DNA extractions were performed on individual pieces of toe stored in 95% ethanol, using a Chelex-based protocol (ESTOUP et al. 1996 ). Ten microsatellite loci (BM118, BM121, BM229, BM235, BM239, BM279, BM322, BM217, BM224, and BM231) were genotyped, using fluorescent PCR and an ABI sequencing machine (Applied Biosystems, Foster City, CA; Perkin-Elmer, Norwalk, CT) as described in TIKEL et al. 2000 .

Estimation procedure for demographic and marker parameters:
Because of the large amount of background information in the present setting, the Bayesian paradigm provides an appropriate framework for statistical inference. In the Bayesian language, the demographic models specify prior distributions for the genealogical tree, and inference then proceeds via the posterior distribution of the tree given the observed data, the coalescent prior for the genealogy, and the priors for the demographic and marker parameters. Posterior distributions tend to give most support to the neighborhood of the true value, although the prior can have an important influence on the posterior. This method allows the data to speak simultaneously for all parameters of interest and to estimate both marginal and joint distributions of those parameters.

The coalescent model (KINGMAN 1982 ; HUDSON 1991 ) provided the theoretical bases required to describe the genealogical processes underlying the patterns of shared ancestry among the genes in the sample(s). In this study, simulations of the coalescent processes were performed according to SIMONSEN et al. 1995 , assuming no gene exchange between the five sampled populations postfounding. Coalescence times were corrected as described in CORNUET and LUIKART 1996 when they overlapped historical periods characterized by different population sizes.

Because of the complexity of our demographic history, we used a nonfully likelihood method based on a rejection algorithm modified from TAVARE et al. 1997 and PRITCHARD et al. 1999 . With this method the allele frequency spectrum of a sample of m chromosomes is summarized by statistics. We chose four summary statistics for microsatellites: a, the mean number of different alleles; H, the mean gene diversity (NEI 1987 ); V, the mean of the variance in repeat number (all three means computed across loci and populations); and F_st as an estimator of differentiation between populations (WEIR and COCKERHAM 1984 ). Enzyme data were summarized by three of these statistics, a, H, and F_st. These statistics were chosen because they exhibit different behaviors under different population histories (e.g., NEI et al. 1975 ; EASTEAL 1985 ; MARUYAMA and FUERST 1985 ; CORNUET and LUIKART 1996 ; PRITCHARD and FELDMAN 1996 ; KIMMEL et al. 1998 ).

Posterior distributions for the parameters of interest were estimated using the following rejection algorithm:

1. Simulate independently values from the distributions of priors (see below) for the six demographic parameters N₁, N₂, N₃, D_B, A_R, and T_S.

2. Simulate, using the coalescent process and the above demographic parameter values, a genealogical tree representing the relationships between the m chromosomes for each of the unlinked loci (10 for microsatellites and 22 for enzymes).

3. Simulate values from the prior distributions (see below) for one or two marker parameters, namely, the mutation rate (µ) and, for microsatellites only, the variance of the geometric distribution for changes in repeat number (²); then simulate genotypes for each of the m/2 individuals at each locus, using those marker parameter values and the genealogical tree.

4. Compute a*, H*, F_st*, and V* (for microsatellites) from the simulated genotypes.

5. If all of |a - a*|/a, |H - H*|/H, |F_st - F_st*|/F_st, and |V - V*|/V (for microsatellites) are less than a small number (), then record N₁, N₂, N₃, D_B, A_R, T_S, µ, and ² (for microsatellites) for the posterior distributions.

6. Return to 1.

This procedure gives a sample from the posterior distributions of N₁, N₂, N₃, D_B, A_R, T_S, µ, and ² (for microsatellites), conditional on a, H, F_st, and V (for microsatellites) being within of the observed values. PRITCHARD et al. 1999 chose values of 0.1 for similar statistics. In this study, simulation showed a relatively low sensitivity of results to values between 0.1 and 0.2 for microsatellites and between 0.05 and 0.1 for enzymes. However, those tests showed that the acceptance rate decreased rapidly with lower values (e.g., for microsatellites, 1/3000 and 1/15,500 for = 0.14 and 0.1, respectively). Therefore, most treatments were done using intermediate values of 0.08 for enzymes (acceptance rate, 1/1800) and 0.14 for microsatellites.

Prior distributions of demographic parameters:
Information available from literature on the cane toad introductions was used to inform prior beliefs about demographic parameters (Table 1). Cane toads reach sexual maturity at 1 year and are then immediately reproductively active (ZUG and ZUG 1979 ). The death rate is high among adults so that the assumption of a generation time of 1 year seems sensible (see also EASTEAL 1985 ; EASTEAL and FLOYD 1986 ).

EASTEAL 1985 used WRIGHT's (1931) classical relationship between the effective population size, the amount of differentiation between isolated populations estimated by F_st, and the number of generations since establishment to estimate the effective population size (roughly corresponding to our parameter N₁) for the five cane toad populations released in Australia and in the Hawaiian islands from 10 polymorphic enzyme loci. This estimation did not consider the large size fluctuations that typically occurred before stabilization of toad populations. EASTEAL 1985 found an effective population size of 350–400 individuals with a "95% confidence interval" between 100 and 800 individuals. This interval was relatively narrow but only considered errors in estimating F_st from genetic samples, but not uncertainty in estimating the effective population size from F_st. Ecological, historical, and demographic data gave neighborhood size estimations of Australian populations between 470 and 46,200 individuals (EASTEAL and FLOYD 1986 ). Although imprecise and referring to a different population model than the one considered by EASTEAL 1985 and by us, those estimations suggest the possibility of large N₁ values. Therefore, we adopted a rather diffuse prior probability distribution on N₁: a lognormal distribution (LN) with parameters (6, 2) and N₁ N₃ (N₃ being the number of founding individuals), which gives support to values from a few to >20,000 individuals.

Less information is available on the effective population size during the boom period (N₂). In several islands, extremely large densities with "hundreds of thousands of individuals" have been mentioned, without much detail on the geographical scale (reviewed in EASTEAL 1981 ). Using mark-recapture methods, FREELAND 1986 showed an increase of density >15 in cane toad populations actively colonizing the northern part of Australia. The low precision of the background information on N₂ was taken into account by choosing a very diffuse prior distribution: a LN(8, 3) distribution with N₂ N₁, which gives support to values of a few tens to several millions of individuals. The boom at population size N₂ is transitory. Its duration (D_B) appears to be 5–30 years (e.g., between 1923 and 1935 in Puerto Rico; EASTEAL 1981 ), with substantial variation among populations (FREELAND 1986 ; A. ESTOUP, personal observations). The prior adopted for D_B was an exponential distribution with a mean of 20 and supported values between 0.5 and 75 years.

The number of individuals introduced in each island varies between 40 (Puerto Rico) and 148 (Hawaii Island) with a mean of 96. We assume the number introduced to Jamaica and Barbados to be in this range. The mean number of effective individuals (N₃) is thus bounded between 2 and 150. N₃ is also likely to be considerably low because of the possibility of an unequal sex ratio, the introduction of immature individuals, and the large reproductive variance between adults. This was taken into account by choosing a LN(3, 1) as prior distribution for N₃ with lower and upper bounds of 2 and 150, respectively. The range of supported values for N₃ was between 3 and 110 (Table 1). We also recorded the founding ratio , the prior on F_R being obtained by combining the prior on N₁ and N₃.

No information was available regarding the level of divergence between the two source populations and the admixture rate in Barbados between those two populations. Genotyping of 10 microsatellite loci from two samples, French Guyana and British Guyana, showed that the source populations were highly differentiated, supporting the hypothesis of a hybrid origin for all introduced island populations (A. ESTOUP and C. MORITZ, unpublished data). The differentiation was measured by computing the mean genetic distance (µ)² (GOLDSTEIN et al. 1995 ), giving (µ)² = 19.140. Assuming a strict stepwise mutation model (SMM; OHTA and KIMURA 1973 ) and a mean mutation rate of 5 x 10^-4 (ESTOUP and ANGERS 1998 ), this value of (µ)² equates to a mean time since separation of the source populations of 19,140 generations (GOLDSTEIN et al. 1995 ), with a 95% confidence interval computed by bootstrapping over loci equal to [830–52,000] generations. Hence, we represented our knowledge of the time since separation of the two source populations (T_S) by an exponential distribution with mean 20,000, truncated at the number of generations difference between the sampling dates in Australia or Hawaiian islands and the date of introduction in Barbados (Fig 1). This distribution supported values between 700 and 73,500 generations (Table 1). In the absence of prior information, we adopted a flat prior for the admixture rate in Barbados (A_R) corresponding to a uniform distribution with lower and upper bounds equal to 0 (no admixture, all individuals originate from a single source population) and 0.5 (50% of individuals originate from each source population), respectively.

Prior distributions of marker parameters:
Interlocus variability in the mutation rate and model (reviewed in ESTOUP and CORNUET 1999 ; and SEIELSTAD et al. 1999 for microsatellites) is another source of noise in an already very noisy system. Information available from literature on the mutation processes of microsatellites and enzymes was used to inform prior beliefs about the parameters describing those processes (Table 1).

Microsatellites: We modeled the mutation rate at microsatellite loci, allowing for variation of mutation rate across loci, while maintaining an appropriate amount of uncertainty in the mean mutation rate. On the basis of data from humans (DIB et al. 1996 ; i.e., 768 mutations observed over 1,238,709 meioses), we used a gamma(768, 1,238,710) prior distribution for the mean mutation rate, with mean 768/1,238,710 6.2 x 10^-4. Hence, in each replicate, we first simulated a mean mutation rate from our gamma(768, 1,238,710) prior, and then, for each locus, the mutation rate was drawn from an exponential distribution with that mean. The mutation rate variance across loci combining these priors is substantial with 2.5 and 97.5% quantiles equal to 1.6 x 10^-5 and 2.3 x 10^-3, respectively. Those values are similar to interval values typically considered for autosomal microsatellites in human (10^-3–10^-5; WEBER and WONG 1993 ; SEIELSTAD et al. 1999 ) and reflect our poor knowledge on the variance of this parameter among loci as well as among species (e.g., SCHUG et al. 1998 ).

On the basis of recent observations (reviewed in ESTOUP and ANGERS 1998 and SEIELSTAD et al. 1999 ; ELLEGREN 2000 ; GONSER et al. 2000 ), we adopted the generalized stepwise mutation model (GSM) in which the change in the number of repeat units forms a geometric random variable (p; KIMMEL and CHAKRABORTY 1996 ). The geometric distribution in our GSM model refers to a change expressed in an (absolute) number of repeat units subsequently added or withdrawn to the mutating allele with equal probability for a symmetric mutation model. The first part of the model is specified by the variance of the geometric distribution ² = ⁺_i=1[i²(1 - p)ⁱ] - = . The large data set of microsatellite mutations for humans of DIB et al. 1996 suggested an estimate of ² near 0.36. Hence, we modeled the variance in mutation size, ², as having an exponential distribution with mean 0.36 as prior.

Enzymes: Information available on the mean and the interlocus variance of the mutation rate is even poorer for enzyme than for microsatellite loci. Studies of the accumulation of allozyme mutants in Drosophila melanogaster lines gave spontaneous mutation rates between 5.14 x 10^-6 (VOELKER et al. 1980 ) and 1.21 x 10^-5 (MUKAI and COCKERHAM 1977 ) or 1.8 x 10^-6 and 1.3 x 10^-6, respectively, if mutations to null alleles are excluded. Similar numbers were obtained for protein loci in humans (NEEL et al. 1986 ). We modeled uncertainty in the mean mutation rate for enzyme by using a gamma(10, 2,000,000) distribution as prior. This distribution has a mean of 5 x 10^-6 and permits most of the above estimates. As for microsatellites, in each replicate, we first simulated a mean mutation rate from the gamma(10, 2,000,000) prior, and then the mutation rate at each locus was drawn from an exponential distribution with that mean. Single locus mutation rates ranged between 1.0 x 10^-7 and 2.0 x 10^-5 (2.5 and 97.5% quantiles, respectively). We chose an infinite allele model (IAM; KIMURA and CROW 1964 ), which does not require any parameters additional to the mutation rate.

Under the standard coalescent, no information about the values of population sizes (N) and µ can be obtained from allelic data except through their product, = 2Nµ. Some information can be inferred for N separately if information about µ and the mutation model can be obtained from other sources, as done in this study. However, if the information on µ and the mutation model is biased or incorrect, it is expected to change our inferences about N as well as those on other demographic parameters such as the admixture rate and time since separation of source populations. To address this issue, simulations were done in which the mutation rates were drawn from a gamma(2, µ/2) prior distribution with µ = 10^-3 or 10^-4 for microsatellites and µ = 10^-5 or 10^-6 for enzymes. Those priors allowed exploration of mutation rate conditions that were substantially different from those used previously while keeping a reasonable variance for µ (e.g., 2.5 and 97.5% quantiles equal 1.2 x 10^-5 and 2.8 x 10^-4, respectively, for µ = 10^-4). In additional test simulations focusing on mutation models, we assumed a SMM for both categories of markers and, for microsatellites only, the occurrence of constraints on allele size or of a mutation bias resulting in the gain of a repeat unit that was more frequent than the loss of a repeat unit (reviewed in ESTOUP and CORNUET 1999 ). Allele size constraints were included in our simulations by imposing reflecting boundaries to the allele size range (e.g., GOLDSTEIN et al. 1995 ; FELDMAN et al. 1997 ). An arbitrary value of k = 20 possible continuous allelic states was chosen as a situation corresponding to strong allele size constraints for microsatellite sequences (e.g., ESTOUP et al. 1999 ). A mutation bias resulting in the gain of a repeat unit that was 2.8 times more frequent than the loss of a repeat unit was suggested by COOPER et al. 1999 .

Tests for a constant demography:
The null hypothesis of constant demography was tested in two different ways. First, a population that has recently experienced a population reduction or a genetic admixture between two differentiated forms will generally develop a heterozygosity excess, while populations that have experienced a population expansion will show a heterozygosity deficit (CORNUET and LUIKART 1996 ). This means that the heterozygosities observed in a sample of the population are larger or smaller than expected assuming mutation-drift equilibrium, a given mutation model, and the same number of alleles in the sample. For microsatellites, an analogous test can be achieved using the allele size variance as test statistic. Expected heterozygosity or allele size variances were estimated assuming a GSM with the variance of the geometric distribution (²) fixed to 0.36. Computations were done using the package Bottleneck (PIRY et al. 1999 ) for the tests on the heterozygosity and a personal program for the test on the allele size variance. The Wilcoxon signed-rank test (SPRENT 1989 ) was used to compute an exact probability over all loci (only polymorphic loci were considered here).

In a second set of computations, we aimed to test whether the enzyme and microsatellite data for the introduced cane toad populations supported our complex model summarized in Fig 1 rather than a model of constant demography. This was done by opposing in a Bayesian framework the two demographic models, the constant demography model being specified by only two (or three) parameters (i.e., N₁, µ, and ² for microsatellites) and the complex model by seven (or eight) parameters (i.e., N₁, N₂, N₃, D_B, A_R, T_S, µ, and ² for microsatellites). As the probabilities of the two models could not be calculated directly, we used acceptance rates for our simulations, which act as an estimate of the posterior probability (PRITCHARD et al. 1999 ). Half of the prior weight was placed on the constant demography model (p_m1 = 0.5) and half on the complex model (p_m2 = 0.5). Then, for a given (0.08 and 0.14 for enzymes and microsatellites, respectively), we estimated the acceptance rates for each model over 2 x 10⁶ and 1 x 10⁶ iterations for microsatellites and enzymes, respectively, and hence an estimate of the relative posterior probabilities of models p_m1 and p_m2. Finally, we applied a correction to take into account the different number of parameters in both models. Because our models are nested, we considered that the statistic G = 2 ln(p_m2/p_m1) followed a ² distribution with the number of degrees of freedom equal to the number of parameters present in one model and absent in the other. This is analogous to a likelihood-ratio test for nested models (HILBORN and MANGEL 1997 ).

	RESULTS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Summary statistics and tests of the constant demography model:
The summary statistics used for the rejection algorithm (number of alleles, gene diversity, allele size variance, and F_st) are given for the Australian and Hawaiian populations and for the microsatellite and enzyme markers in Table 2.

As expected, microsatellites (9 polymorphic and 1 monomorphic loci) showed a higher level of variability than enzymes (10 polymorphic and 12 monomorphic loci). The slightly higher F_st values for microsatellites than for enzymes could reflect the more recent sampling time for microsatellites. Note that considering polymorphic loci only instead of all loci change substantially the values of the statistics other than F_st (Table 2), illustrating the bias that could occur if only polymorphic markers were selected. In the following analyses (except for the test of CORNUET and LUIKART 1996 ), only summary statistics computed from all markers were considered.

For both source populations, no significant deviations from mutation-drift equilibrium were detected using heterozygosity or allele size variance as test statistics (Table 3).

View this table: In this window In a new window   Table 3. Test for mutation-drift equilibrium using the method of CORNUET and LUIKART 1996

On the other hand, for the introduced populations, all 10 populations for enzymes and 3 of 5 for microsatellites showed a significant excess of heterozygosity, and all 5 populations showed a significant excess of microsatellite allele size variance (Table 3). These deviations from equilibrium are in agreement with at least two demographic events: a bottleneck associated with founding and the mixing of individuals from two differentiated source populations in Barbados. No signal of recent population expansions (deficit of heterozygosity and/or allele size variance) that could have reflected the boom-and-bust phenomenon following each introduction could be detected in any introduced population.

Our Bayesian hypothesis test shows that all posterior support was on our complex introduction scenario (p_m2 > 0.998 and 0.996 for microsatellites and enzymes, respectively), when opposed to the model of constant demography. The constant demography model was still rejected in favor of our complex introduction scenario when applying a correction taking into account the different number of parameters in both models (G > 12.84, d.f. = 5, p < 0.025). Moreover, opposing our complex demographic model to the same model without admixture shows that the posterior support was much stronger for the complex introduction scenario with admixture for both microsatellites (G = 6.87, d.f. = 2, p = 0.032) and enzymes (G > 12.44, d.f. = 2, p < 0.002).

Inferences from posterior distributions:
The density curves as well as the mean and the 5 and 95% quantiles are given in Fig 2 and Fig 3 for the prior and the posterior distributions of each of the demographic and marker parameters.

View larger version (26K): In this window In a new window Download PPT slide   Figure 2. Posterior and prior density curves for the demographic parameters N1, N3, and FR for microsatellite and enzyme data, together with corresponding mean values, and 5, 50, and 95% quantiles. The solid lines correspond to the prior density curves. The dashed, dotted, and dot-dashed lines correspond to the posterior density curves for microsatellites in Australian populations and for enzymes in Australian and Hawaiian populations, respectively. All prior and posterior densities are based on 20,000 and 1000 values, respectively. A graph showing the relationship between N1 and N3 values (100 values randomly sampled from the 1000 posterior values) is also given. Diamonds, crosses, and solid circles are for microsatellite data and Australian and Hawaiian enzyme data, respectively.

View larger version (74K): In this window In a new window Download PPT slide   Figure 3. Posterior and prior density curves for the demographic parameters N2, DB, BR, AR, and TS, and for the marker parameters µ, , and 2, together with corresponding mean values, and 5, 50, and 95% quantiles. The solid lines correspond to the prior density curves. The dashed, dotted, and dot-dashed lines correspond to the posterior density curves for microsatellites in Australian populations and for enzymes in Australian and Hawaiian populations, respectively. All prior and posterior densities are based on 20,000 and 1000 values, respectively.

In describing our results, we regard the posterior median (50% quantile) of each parameter as a point estimate (e.g., WILSON and BALDING 1998 ). Since the data may not be informative about the demographic parameters individually, a graph illustrating the relationship between N₁ and N₃ values is also given in Fig 2. Due to space limitations, similar graphs are not shown for other parameters. With the exception of the parameters dealing with the boom-and-bust events, the posterior density curves differ noticeably from the priors. This means that the genetic data contained a substantial amount of information for at least some demographic parameters. The posterior distributions were very similar for the Australian and Hawaiian enzymes data so that both data sets are hereafter referred to as the enzyme data set.

For both microsatellites and enzymes the posteriors for the stable effective population size (N₁) were sharper than the priors, with diminished support for high and low prior values resulting in a much narrower 95% interval (Fig 2). The posterior distribution was slightly narrower and with a lower median for microsatellite than for allozyme loci. Highest support was for N₁ values of a few hundred effective individuals (point estimates 320–420). These estimates are similar to the effective population size estimated by EASTEAL 1985 but are more realistic as mutation and population size fluctuations are considered here (Fig 1). The distribution of posteriors relative to the prior for N₃, the effective founder size, was markedly narrower with reduced support for high values for microsatellites and broader with more support for larger values for enzymes (Fig 2). However, for both types of marker, the upper 95% quantile for N₃ is less than the lower 5% quantile for N₁, confirming that a bottleneck has occurred in association with introductions. Accordingly, 5% quantiles for the founder ratio () were larger than the prior one. Nevertheless, 95% quantiles for F_R were much lower for the posteriors than the prior, indicating that the bottleneck is of moderate intensity, especially for enzymes. The graph in Fig 2 also shows a different relationship between N₁ and N₃ for microsatellites and enzymes: The correlation between N₁ and N₃ appears to be negative for microsatellites and positive for enzymes, and the range of N₃ values associated with a given range of N₁ values is larger for enzymes than for microsatellites. Posterior distributions for admixture rate and time since isolation were shifted toward higher values, giving support to the occurrence of genetic admixture in Barbados between two substantially differentiated source populations (Fig 3). The shift of posterior distributions was considerably stronger for enzymes than microsatellites.

The posterior distributions for N₂ and the boom duration were similar to the corresponding prior distributions. Moreover, posterior distributions for each of those parameters changed substantially when using different priors, with posterior shapes tracking those of the priors (results not shown). Thus, the genetic data contained very little information on the parameters dealing with the boom-and-bust events.

Regarding marker parameters, posteriors for the mutation rate µ (as well as for the mean mutation rates ) were shifted toward larger values for enzymes (point estimates around µ = 4.5 x 10^-6 vs. 3.3 x 10^-6 for the prior). In contrast a shift toward smaller µ values was observed for microsatellites (point estimates 3.1 x 10^-4 vs. 4.3 x 10^-4 for the prior), although this shift was not apparent for (Fig 3). The combination of summary statistics is expected to contain some information about the variance of the geometric distribution of the GSM assumed for microsatellites (²). However, our data did not allow any inference on this parameter, the prior and posterior distributions being very similar (Fig 3).

Sensitivity to marker and demographic priors:
The results of simulations testing the effect of the mutation rates and models of markers on historical inferences are summarized in Table 4 for four demographic parameters for which significant information could be previously obtained.

A first major conclusion is that the demographic trends previously observed still hold for all additional simulations, indicating that our inferences are relatively robust to prior beliefs on the mutation parameters. However, both stronger allele size constraints and especially lower mutation rates tended to shift posterior support toward a less intense bottleneck and a more balanced admixture. The posterior shifts observed with different mutation rates are larger for microsatellites than enzymes. The effect of allele size constraints at microsatellites is more pronounced on parameters dealing with the admixture event (A_R and T_S) than on those dealing with the founder events, while the effect is substantial on both events for the mutation rate. The assumption of a SMM rather than an IAM for enzymes and a GSM for microsatellites had very little effect on posterior distributions. The only exception concerns the separation time (T_S) between the two source populations, which was larger for microsatellites under a SMM than a GSM. Including a positive mutation bias in the microsatellite mutation model did not affect any posterior distributions.

The robustness of the main conclusions of this study to the prior distribution of the demographic parameters was evaluated using entirely flat priors (i.e., uniform distributions) for all demographic parameters: a U[2–150] for N₃, a U[2–2000] with N₁ N₃ for N₁, a U[2–50,000] with N₂ N₁ for N₂, a U[0–50] for D_B, a U[T_i–60,000] with T_i the number of generations since introduction in Barbados for T_S, and a U[0–0.5] for A_R. As expected, posterior distributions were affected by prior changes (curves not shown but see Table 4). However, the posterior shifts supported again the occurrence of genetic admixture in Barbados and of bottlenecks of moderate intensity in each island. Those additional simulations also confirmed the stronger support for a larger bottleneck at introductions for microsatellites than for enzymes and a more balanced genetic admixture for enzymes than for microsatellites. Note that acceptance rates for microsatellites were two to three times smaller with the initial priors, making an already lengthy treatment even more time consuming.

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Resolution of different demographic events:
The absence of support for transitory population explosions does not necessarily imply that the large increase of density observed in new cane toad populations does not translate into a large increase in the effective number of individuals. Change of the prior distribution results in similar change in the posterior distribution of N₂, such that the posterior tracks the prior, indicating that the data contain very little information on this demographic parameter. Poor resolution is expected for large values of effective population size during a population flush of short duration (D_B << N₁ and N₂), since different (large) N₂ values generate a similar effect in terms of drift or number of coalescent events per generation (i.e., reduced drift and rate of coalescence). The pattern of new allelic variants arising through mutation should contain information about the population size; however, observing this pattern requires a sufficiently large number of generations to allow new variants to be generated. It is likely that the ages of the introduction events (167 generations) and the duration of each population explosion are insufficiently large even for highly mutable markers such as microsatellites. Moreover, the occurrence of founder events after each expansion period makes it difficult to maintain low frequency new variants or, rather, allows the quick emergence of those new variants at high frequency.

This limit should not hold for founder (or bottleneck) events since the detection and quantification of such demographic events does not require the occurrence of new mutations. The effect of bottleneck(s) is essentially a rapid erosion of the existing allelic diversity through the loss of alleles, a feature that should be quickly detectable for any category of markers. In agreement with this, both microsatellite and enzyme data support the occurrence of effective population size reduction(s). In a similar way, the detection of admixture between differentiated taxa relies on differentiation between source populations accumulated as mutations during their independent evolution, rather than the retention of new allelic variants, so that its detection and quantification should be possible early in the history of populations, a prediction also confirmed by our results.

Discrepancies between microsatellite and enzyme markers:
The discrepancies observed for posterior distributions for microsatellite and enzyme markers were particularly large for the demographic parameters related to founding and admixture events. This may reflect the occurrence of particular evolutionary features at one or both categories of markers not taken, or wrongly taken, into account by our model.

For microsatellites, test simulations showed that our perception of the intensity of the founder and admixture events may be biased, using a prior on mutation rate that tended to overestimate mutation rates and/or by the occurrence of allele size constraints. Strong allele size constraints and especially lower mutation rates both tend to shift posterior support toward a less intense bottleneck and a more balanced admixture. This is because lower mutation rates and (strong) allele size constraints both reduce the number of alleles as well as the allele size variance so that less intense bottleneck and more intense admixture are required to explain an excess of heterozygosity and of allele size variance relative to the number of alleles in the gene samples. The cooccurrence of lower mutation rate and strong allele size constraints could hence explain, at least partly, the discrepancy between some posterior distributions for microsatellites and enzymes. The mutation model itself (i.e., SMM vs. GSM) does not explain those discrepancies except, to a certain extent, for the time since separation of the two source populations.

Our simulations also showed that strong founder events may be undersupported and admixture oversupported by the use of a prior on enzyme mutation rate that underestimates mutation rates. Designing mutation rate priors for cane toad enzymes from Drosophila data is far from optimal due to potential variation between species. An additional issue may be that, to estimate mutation rates, amidon gels were used to reveal new variants in Drosophila (MUKAI and COCKERHAM 1977 ; VOELKER et al. 1980 ), whereas EASTEAL 1985 used higher resolution cellulose acetate gels for cane toad.

An alternative but not exclusive explanation for the discrepancies observed between microsatellites and enzymes could be the occurrence of different selective pressure on both categories of markers. If selection in favor of "hybrid" individuals (hybrid vigor, e.g., INGVARSSON and WHITLOCK 2000 ) acted more intensively at enzyme than at microsatellite loci, a selectively neutral model applied to enzyme data would overestimate support for large admixture rate and separation time of the source population and underestimate support for intense founder events that would eliminate alleles.

Limits of our model and the rejection method:
Fully realistic models for the introduction and demographic history of cane toad remain inevitably outside our grasp. We had to assume an equal stable effective population size and number of founding individuals in all islands to make our model tractable with the present data set. The number of demographic parameters in our model is already large and considering different number of founders for each island would introduce compensatory effects between N₁ and/or N₃ values. A second potential drawback of our model is that the exact location of introduction is unknown in several cases. Since the colonization process of cane toad is likely to involve founder events (LEBLOIS et al. 2000 ), the collection of samples outside the release points may artificially increase the differentiation between introduced populations (SLATKIN 1993 ; LECORRE and KREMER 1998 ) and hence underestimate N₁. However, an absence of evidence for isolation by distance was found in a relatively old Australian population (25–35 generations) of cane toad (LEBLOIS et al. 2000 ). This suggests that the differentiation observed at a local scale between recently founded populations by EASTEAL 1985 is transitory and is unlikely to persist in populations of age 48 or 63 generations such as those analyzed in our study.

We tried to take into account some of the complexity (and uncertainty) of the marker mutation processes, especially the variance of mutation rate between loci and the occurrence of multistep mutations for microsatellites. The true mutation processes for microsatellites are probably more complex than the model used here (reviewed in ESTOUP and CORNUET 1999 ). For instance, there is growing evidence for mutational bias with gains in repeats outnumbering losses (e.g., PRIMMER et al. 1998 ; ELLEGREN 2000 ). However, test simulations showed that this particular mutation feature did not have any significant effect on our inferences (Table 4).

The large number of summary statistics used in our rejection algorithm (four for microsatellites and three for enzymes) was necessary to tackle the complexity of our demographic scenario. These statistics are differentially sensitive to different aspects of the sorts of population fluctuation in our demographic model, and the balance between some of them is also potentially informative, as shown by previous studies [e.g., the heterozygosity vs. the number of alleles (CORNUET and LUIKART 1996 ); the allele size variance vs. the heterozygosity (KIMMEL et al. 1998 )]. Running simulations without one of the summary statistics (for instance the allele size variance or the F_st among populations) but with a lower threshold did not increase the precision level of any demographic parameters and actually decreased the precision level of some of them (results not shown).

The large number of sampled genes and markers for microsatellites and especially for enzymes (Table 2) limited the tractability of our treatments. Getting 1000 posterior values for each parameter took 3 days for microsatellites and 7 days for enzymes, when using a standard 350-Mhz PC platform. These limits of tractability were particularly obvious when we attempted to treat the microsatellite and enzyme data together (i.e., as a single data set of 32 loci and 8000 genes). In addition to the problem of dealing with the large dimension of this data set, the large number of conditional statistics (i.e., seven), and the discrepancies between microsatellites and enzymes, posteriors generated such a low acceptance rate that this treatment turned out to be intractable.

Conclusion:
This study is among the first to tackle a biologically realistic scenario incorporating the increased complexities that are usually ignored. Including prior information on demography and history has allowed a realistic Bayesian analysis of recent demographic events and cut the number of parameters to a manageable level. Although substantial uncertainties remain on several parts of our demographic model and discrepancies between microsatellite and enzyme markers were observed, our treatments allowed some important historical inferences to be made with reasonable confidence. In this context, continued reliance on equilibrium assumptions for blatantly nonequilibium systems (e.g., declining or expanding populations) should become less acceptable. Computational resources may provide a barrier to extending the type of analysis done in this article to more sophisticated modeling assumptions and larger sample size. One possible methodological solution that could represent a considerable improvement both in terms of precision and speed would consist of using computational techniques such as Markov chain Monte Carlo (MCMC; e.g., WILSON and BALDING 1998 ; BEERLI and FELSENSTEIN 1999 ). However, although this is theoretically possible, implementing MCMC algorithms on demographic scenarios as complex as the one in our study remains a difficult challenge. Studies evaluating potential evolutionary consequences of complex founder/boom and admixture events in any species with historical and demographic information clearly would benefit from methods such as those applied in this article.

	ACKNOWLEDGMENTS

We thank D. Tikel, C. Dutech, R. Leblois, and V. Olson for assistance in toad sampling; S. Easteal for providing enzyme data and useful discussion on cane toad population genetics; Stuart Baird for stimulating discussions and help in the simulation of gamma distributions; J. Pritchard for advice on various theoretical aspects; and M. Slatkin and two anonymous referees for constructive comments on the manuscript. This work was supported by a Special Investigator Award to C.M. from the Australian Research Council and grants to A.E. from the Institut National de Recherche Agronomique and the French Ministere de l'Environnement on bioinvading systems.

Manuscript received April 3, 2001; Accepted for publication September 21, 2001.

	LITERATURE CITED

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

ALFORD, R. A., M. P. COHEN, M. R. CROSSLAND, M. N. HEARNDEN, D. JAMES et al., 1995 Population biology of Bufo marinus in Northern Australia, pp. 173–181 in Wetland Research in the Wet-Dry tropics of Australia, edited by M. FINLAYSON. Office of the Supervising Scientist Report 101, Canberra, Australia.

BAHLO, M. and R. C. GRIFFITH, 2000  Inference from gene trees in a subdivided population. Theor. Popul. Biol. 57:79-95[Medline].

BEAUMONT, M., 1999  Detecting population expansion and decline using microsatellites. Genetics 153:2013-2029[Abstract/Free Full Text].

BEERLI, P. and J. FELSENSTEIN, 1999  Maximum-likelihood estimation of migration rates and effective population numbers in two populations using a coalescent approach. Genetics 152:763-773[Abstract/Free Full Text].

COOPER, G., N. BURROUGHS, D. A. RAND, D. C. RUBINSZTEIN, and W. AMOS, 1999  Markov chain Monte Carlo analysis of human Y-chromosome microsatellites provides evidence of mutation bias. Proc. Natl. Acad. Sci. USA 96:11916-11921[Abstract/Free Full Text].

CORNUET, J.-M. and G. LUIKART, 1996  Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics 144:2001-2014[Abstract].

DIB, C., S. FAURE, C. FIZAMES, D. SAMSON, and N. DROUOT et al., 1996  A comprehensive map of the human genome based on 5,264 microsatellites. Nature 380:152-154[Medline].

DONNELLY, P. and S. TAVARÉ, 1995  Coalescents and genealogical structure under neutrality. Annu. Rev. Genet. 29:410-421.

EASTEAL, S., 1981  The history of introductions of Bufo marinus (Amphibia: Anura); a natural experiment in evolution. Biol. J. Lin. Soc. 16:93-113.

EASTEAL, S., 1985  The genetics of introduced populations of the giant toad Bufo marinus: III. Geographical patterns of variation. Evolution 39:1065-1075.

EASTEAL, S. and R. D. FLOYD, 1986  The ecological genetics of introduced populations of the giant toad, Bufo marinus (Amphibia: Anura): dispersal and neighborhood size. Lin. Soc. Lond. 27:17-45.

ELLEGREN, H., 2000  Heterogeneous mutation processes in human microsatellite DNA sequences. Nat. Genet. 24:400-402[Medline].

ESTOUP, A., and B. ANGERS, 1998 Microsatellites and minisatellites for molecular ecology: theoretical and empirical considerations, pp. 55–86 in Advances in Molecular Ecology, edited by G. CARVALHO. NATO Press, Amsterdam.

ESTOUP, A., and J.-M. CORNUET, 1999 Microsatellite evolution: inferences from population data, pp. 50–65 in Microsatellites: Evolution and Applications, edited by D. B. GOLDSTEIN and C. SCHLÖTTERER. Oxford University Press, Oxford.

ESTOUP, A., C. L. LARGIADER, E. PERROT, and D. CHOURROUT, 1996  Rapid one-tube DNA extraction for reliable PCR detection of fish polymorphic markers and transgenes. Mol. Mar. Biol. Biotech. 5:295-298.

ESTOUP, A., F. ROUSSET, J.-M. CORNUET, and R. GUYOMARD, 1999  Juxtaposed microsatellites as diagnostic markers for genetic introgression: theoretical aspects. Mol. Biol. Evol. 16:898-908.

FELDMAN, M. W., A. BERGMAN, D. D. POLLOCK, and D. B. GOLDSTEIN, 1997  Microsatellite genetic distances with range constraints: analytic description and problems of estimation. Genetics 145:207-216[Abstract].

FREELAND, W. J., 1986  Population of the cane toad, Bufo marinus, in relation to time since colonization. Aust. Wildl. Res. 13:321-330.

FU, Y.-X. and W.-H. LI, 1999  Coalescing into the 21^st century: an overview and prospects of coalescent theory. Theor. Popul. Biol. 56:1-10[Medline].

GOLDSTEIN, D. B., A. R. LINARES, M. W. FELDMAN, and L. L. CAVALLI-SFORZA, 1995  Genetic absolute dating based on microsatellites and the origin of modern humans. Proc. Natl. Acad. Sci. USA 92:6723-6727[Abstract/Free Full Text].

GONSER, R., P. DONNELY, G. NICHOLSON, and A. DI RIENZO, 2000  Microsatellite mutations and inferences about human demography. Genetics 154:1793-1807[Abstract/Free Full Text].

HILBORN, R., and M. MANGEL, 1997 The Ecological Detective: Confronting Models With Data. Princeton University Press, Princeton, NJ.

HUDSON, R. R., 1991 Gene genealogies and the coalescent process, pp. 1–44 in Oxford Surveys in Evolutionary Biology, edited by D. J. FUTUYAMA and J. ANTONOVICS. Oxford University Press, Oxford.

INGVARSSON, P. K. and M. C. WHITLOCK, 2000  Heterosis increases the effective migration rate. Proc. R. Soc. Lond. Ser. B 267:1321-1326[Medline].

KIMMEL, M. and R. CHAKRABORTY, 1996  Measures of variation at DNS repeat loci under a general stepwise mutation model. Theor. Popul. Biol. 50:345-367[Medline].

KIMMEL, M., R. CHAKRABORTY, J. P. KING, M. BAMSHAD, and W. S. WATKINS et al., 1998  Signatures of population expansion in microsatellite repeat data. Genetics 148:1921-1930[Abstract/Free Full Text].

KIMURA, M. and J. F. CROW, 1964  The number of alleles that can be maintained in a finite population. Genetics 49:725-738[Free Full Text].

KINGMAN, J. F. C., 1982  The coalescent. Stoch. Proc. Appl. 13:235-248.

LEBLOIS, R., F. ROUSSET, D. TIKEL, C. MORITZ, and A. ESTOUP, 2000  Absence of evidence for isolation by distance in an expanding cane toad (Bufo marinus) population: an individual-based analysis of microsatellite genotypes. Mol. Ecol. 9:1905-1909[Medline].

LE CORRE, V. and A. KREMER, 1998  Cumulative effects of founding events during colonisation on genetic diversity and differentiation in an island and stepping-stone model. J. Evol. Biol. 11:495-512.

MARUYAMA, T. and P. A. FUERST, 1985  Population bottlenecks and nonequilibrium models in population genetics. II. Number of alleles in a small population that was formed by a recent bottleneck. Genetics 111:675-689[Abstract/Free Full Text].

MUKAI, T. and C. C. COCKERHAM, 1977  Spontaneous mutation rates at enzyme loci in Drosophila melanogaster. Proc. Natl. Acad. Sci. USA 74:2514-2517[Abstract/Free Full Text].

NEEL, J. V., C. SATOH, K. GORIKI, M. FUJITA, and N. TAKAHASHI et al., 1986  The rate with which spontaneous mutation alters the electrophoretic mobility of polypeptides. Proc. Natl. Acad. Sci. USA 83:389-393[Abstract/Free Full Text].

NEI, M., 1987 Molecular Evolutionary Genetics. Columbia University Press, New York.