Abstract
Drift and migration disequilibrium are very common in animal and plant populations. Yet their impact on methods of estimation of demographic parameters was rarely evaluated especially in complex realistic population models. The effect of such disequilibria on the estimation of demographic parameters depends on the population model, the statistics, and the genetic markers used. Here we considered the estimation of the product Dσ^{2} from individual microsatellite data, where D is the density of adults and σ^{2} the average squared axial parentoffspring distance in a continuous population evolving under isolation by distance. A coalescencebased simulation algorithm was used to study the effect on Dσ^{2} estimation of temporal and spatial fluctuations of demographic parameters. Estimation of presenttime Dσ^{2} values was found to be robust to temporal changes in dispersal, to density reduction, and to spatial expansions with constant density, even for relatively recent changes (i.e., a few tens of generations ago). By contrast, density increase in the recent past gave Dσ^{2} estimations biased largely toward past demographic parameters values. The method was also robust to spatial heterogeneity in density and estimated local demographic parameters when the density is homogenous around the sampling area (e.g., on a surface that equals four times the sampling area). Hence, in the limit of the situations studied in this article, and with the exception of the case of density increase, temporal and spatial fluctuations of demographic parameters appear to have a limited influence on the estimation of local and presenttime demographic parameters with the method studied.
DISPERSAL rates and population sizes or densities are important demographic parameters in evolutionary processes. Many studies have attempted to estimate those parameters, using direct methods (e.g., markrecapture methods) or indirect methods (genetic markers). Discrepancies between estimations based on direct and indirect methods have often been attributed to inadequacies of the assumptions of the genetic models in indirect methods (Hastings and Harrison 1994; Slatkin 1994; Koeniget al. 1996). The assumptions that have usually been considered inadequate are those related to the modalities of dispersal (e.g., the island model), the mutation rates and processes of genetic markers, the selective neutrality of genetic markers, and the demographic stability in time and space. The latter assumption raises the question of the exact meaning of demographic parameter estimations in biological systems for which temporal and/or spatial fluctuations of demographic parameters have occurred. With a few exceptions (e.g., Stone and Sunnucks 1993; Beebee and Rowe 2001; Spong and Hellborg 2002), population geneticists usually consider that contemporary spatial patterns of diversity reflect the past more than the presenttime population dynamics of a species. Whitlock and McCauley (1999) recently concluded that estimates of the number of migrants between subpopulations from Fstatistics under the assumption of an island model at equilibrium were “likely to be correct within a few orders of magnitude” only because assumptions of the genetic model (i.e., equal migration, no selection, and demographic stability) are often violated in biological systems. This degree of precision is of little value for understanding the presenttime demographic processes of populations. This is particularly worrying in a practical context since reliable estimates of present or at least recent migration rates, dispersal distances, or densities are increasingly demanded as integral elements of applied management and conservation decisions.
The effect of temporal and spatial fluctuations on the estimation of demographic parameters strongly depends on the type and intensity of the fluctuation encountered. However, it also strongly depends on the population models assumed, the statistics computed, and the genetic markers used. Most studies dealing with disequilibrium situations referred to the classical island model or to the WrightFisher population model and only a few of them have considered more sophisticated and realistic models (but see Slatkin 1993). In numerous species, individual dispersal is restricted in space (see references in Lebloiset al. 2003). A method of analysis adapted to a “continuous” population evolving under isolation by distance was developed to estimate the product Dσ^{2}, where D is the density of adults and σ^{2} the average squared axial parentoffspring distance (Rousset 2000). This method uses a regression of estimators of a parameter a_{r} to the geographical distances or the logarithm of the geographical distances in one or two dimensions, respectively. The parameter a_{r}, defined in Rousset (2000), is analogous to the parameter F_{ST}/(1  F_{ST}) but is calculated between individuals (see Method of analysis for details about this parameter and its estimator). The inverse of the slope of the regression line gives an estimate of 4πDσ^{2} (Rousset 1997). The method is valid for leptokurtic distributions of dispersal distance (Rousset 2000; Lebloiset al. 2003), a feature commonly observed in natural populations (review and data in Endler 1977; Portnoy and Willson 1993). Because analysis of genetic differentiation is made at a small (local) geographical scale, heterogeneity of demographic parameters such as dispersal or density is reduced and hence its influence on genetic differentiation is also reduced (Slatkin 1993; Rousset 2001). The good properties of this method have been confirmed by comparisons of direct and indirect estimates of Dσ^{2} (Rousset 2000; Sumneret al. 2001).
As for any population genetics method of demographic parameter estimation, the quality of the estimation of Dσ^{2} using this method may be affected by local and temporal spatial heterogeneities in demographic parameters. In this study, we adapted the coalescencebased simulation algorithm of Leblois et al. (2003) to study the effect of temporal and spatial fluctuations of demographic parameters on the estimation of presenttime Dσ^{2}. Although one can imagine many scenarios dealing with demographic heterogeneities in space and time, we have chosen to focus our study on demographic scenarios often met in empirical surveys in conservation biology and in the study of introduced invading species. In this context, we assessed the effect on the estimation of the presenttime Dσ^{2} of (i) a temporal change of the dispersal feature, (ii) a density reduction (bottleneck) or increase (flush) in time, (iii) a spatial expansion with constant density, and (iv) a sample of individuals taken from a highdensity zone within a lowerdensity area.
MODELS AND METHODS
Spatial model and population cycle: The model that we considered for “continuous” populations is the lattice model with each lattice node corresponding to one diploid individual. This model without demic structure is viewed as an approximation for truly continuous populations with infinitely strong density regulation (Malécot 1975; Rousset 2000). More realistic continuous models would incorporate the feature that individuals could settle in any position in a continuous space. Although such models have been formulated (e.g., Malécot 1967; Sawyer 1977), it is known that they do not follow a welldefined set of biological assumptions (Maruyama 1972; Felsenstein 1975; see Bartonet al. 2002 for an alternative approach for continuous populations). To avoid edge effects, a twodimensional lattice is represented on a torus. Edges and lattice size have little effect on local differentiation when the habitat area (i.e., the lattice size) is large compared to the mean dispersal (Lebloiset al. 2003). Finally, we considered diploid individuals with dispersal through gametes only. The life cycle is divided into five steps: (i) at each reproductive event, each individual gives birth to a great number of gametes and dies; (ii) gametes undergo the effect of mutations; (iii) gametes disperse; (iv) diploid individuals are formed; and (v) competition brings back the number of adults in each deme to N (usually N = 1 but see Spatial and temporal heterogeneities). We assume here random assortment of gametes present after dispersal at a given node. This is akin to random selfing in a population of N diploids without spatial structure, by which selfing occurs with frequency 1/N. How alternative assumptions would affect the analysis is discussed below.
Coalescent algorithm: In this work, we focused on isolation by distance. For this category of models, no analytical treatment of coalescence time or coalescence probabilities has been done for more than two genes. The coalescent algorithm used in this study is thus not based on the largeN approximation of the ncoalescent theory; rather it is an exact algorithm for which coalescence and migration events are considered generation by generation until the common ancestor of the sample has been found. The idea of tracing lineages back in time generation by generation is fundamental in the coalescence theory, and is well described in Nordborg (2001). Such a generationbygeneration algorithm leads to less efficient simulations in terms of computation time than do those based on the ncoalescent theory (Kingman 1982a,b; Nordborg 2001). However, this algorithm is much more flexible when complex demographic and dispersal features are considered. Note that, since multiple coalescent events are taken into account by considering the probability of a coalescence event of k genes in a given parental node (= 1/2^{k}^{1} under the model with one individual per lattice node), it allows us to build an exact coalescent tree under very small population size. The entire generationbygeneration algorithm that gives the coalescent tree for a sample of n genes evolving under isolation by distance, with density and dispersal homogenous in space and time, is detailed in Leblois et al. (2003). The algorithm and the program used in this study were checked at every step during its elaboration by comparing simulated values of probabilities of identity of two genes under models of isolation by distance on finite lattices with their exact analytically computed values (e.g., Malécot 1975 for the lattice model) with adaptation to different mutation models following general methods valid for any assumption about dispersal and density (Rousset 1996). These comparisons show that estimates of identity probabilities from our program and analytical expectations differ by less than one per thousand for sufficiently long runs.
Dispersal functions: Let (dx, dy) be the parentoffspring axial distance, backward in time, expressed in number of steps on the lattice. Under a twodimensional model, the probability distribution of the random variable (dx, dy) is given by b_{dx}_{,}_{dy}, the “backward” dispersal function. The term backward is used because the position of the parental gene is determined knowing the position of its descendant gene. This function is calculated using f_{dx}_{,}_{dy}, the forward dispersal density function describing where descendants go. Biologically realistic dispersal functions often have a high kurtosis (Endler 1977; Kotet al. 1996). As previously explained (Rousset 2000), the commonly used discrete probability distributions for dispersal are not appropriate here because high kurtosis can be achieved only by assuming a low dispersal probability, i.e., that most offspring reproduce exactly where their parents reproduced. Thus we used forward dispersal distributions for which the probability of moving k steps (for 0 < k ≤ K_{max}) in one direction is of the form
Dispersal was assumed to be independent in each direction, so that f_{dx}_{,}_{dy} = f_{dx} × f_{dy}. When density is homogenous in space, backward dispersal functions are equal to forward dispersal functions, so that b_{dx}_{,}_{dy} = f_{dx}_{,}_{dy} = f_{dx} × f_{dy}.
Mutation processes: The number of mutations on each branch of the coalescent tree follows a binomial distribution with parameter (μ, L), where μ is the mutation rate and L the length of the branch. The allelic states of each gene of the sample were obtained starting from the common ancestor of the sample (root of the genealogical tree) from an allelic state determined according to a probability distribution determined by the mutation model and then going forward in time adding mutations one by one on each branch of the tree. The study of Leblois et al. (2003) stressed the interest in using loci with high levels of polymorphism for Dσ^{2} estimation. Therefore, microsatellite markers were simulated in the present study. On the basis of direct observations of mutations at human microsatellite loci (Dibet al. 1996; Ellegren 2000), the generalized stepwise model (GSM) in which the change in the number of repeat units forms a geometric random variable was adopted (Pritchardet al. 1999; Estoupet al. 2001). The variance of the geometric distribution was fixed at 0.36 (Estoupet al. 2001), a value computed from the mutation data in Dib et al. (1996). The mutation rate was equal to 5 × 10^{4}, a value considered as the average mutation rate in many species (reviewed in Estoup and Angers 1998). The GSM does not capture all the complexity of the mutation process at microsatellite loci (reviewed in Ellegren 2000; Schlötterer 2000). However, Leblois et al. (2003) have shown that exact mutation processes, and in particular the occurrence of constraints on allele size and increase of mutation rate with allele length, have little influence on Dσ^{2} estimations.
Method of analysis: Each simulation iteration gives the genotypes at 10 polymorphic loci of 100 (i.e., 10 × 10) individuals characterized by their coordinates on the lattice. Ten loci and 100 individuals were considered as representative of the number of loci and individuals commonly analyzed in empirical studies based on microsatellites. Independent coalescent trees were used to simulate multilocus genotypes at independent loci. In practice it is difficult to sample all individuals in a small area. Simulations were run for a sample of (10 × 10) individuals taken every two nodes from an area of (20 × 20) nodes in the lattice. In this we aimed to roughly mimic a sampling scheme commonly achieved in empirical studies. This process was repeated 1000 times giving 1000 multilocus samples of 100 individuals sharing the same demographic history.
For each simulated multilocus sample, estimates of the parameter a_{r} = (Q_{w}  Q_{r})/(1  Q_{w}) were computed for each pair of individuals, with Q_{w} the probability of identity in state for two genes taken from the same individual and Q_{r} the probability of identity in state for two genes at geographical distance r (Rousset 2000). The parameter a_{r} is a parameter analogous to F_{ST}/(1 F_{ST}) calculated between individuals (not between populations as in Rousset 1997). An estimator of a_{r} for a pair ξ of individuals taken from the P different possible pairs is
The generalized random selfing assumption made in this article implies that the identity within individuals is identical to the identity between juveniles competing for a site. More generally, Dσ^{2} is related to the parameter
For each simulated data set, the value of the slope of the regression line between â and the logarithm of geographical distance was computed. In the limit of low mutation rates, the inverse of the slope is an estimate of the product 4πDσ^{2} (Rousset 1997). High mutation rates should not result in a large sample bias as long as one focuses on local processes involving distances between sampled individuals,
The quality of an estimator is usually assessed through the computation of its bias and its mean square error (MSE). These measures are suitable when estimates have an approximately normal distribution but not when estimates are sometimes infinite. In the present case, a negative slope should be interpreted as an infinite estimate of Dσ^{2}. Therefore, we present the bias and the MSE for the slope values of the regression lines and not for Dσ^{2} estimates. Thus, the following statistics were estimated over all repetitions: (i) the mean relative bias between the value of the slope and the expected value, 1/(4πDσ^{2}) [i.e., (observed slope  expected slope)/expected slope]; (ii) the standard error on this relative bias; and (iii) the mean square error [i.e., MSE = ((observed slope  expected slope)/expected slope)^{2}]. The bias and the MSE are relative values since they are computed from the ratio of the observed to the expected value. We also computed the probability that the estimate of 1/(4πDσ^{2}) was within a factor of two from the expected value (i.e., in the interval [expected slope/2; 2 × expected slope]).
Spatial and temporal heterogeneities: One important advantage of the generationbygeneration algorithm is that virtually any demographic model including those with variations in time and space of demographic parameters can be easily implemented.
Temporal change in dispersal: We first studied the effect of a simple decrease of dispersal capabilities in time. Decrease in dispersal under isolationbydistance models can be modeled in various ways (i.e., changing various parameters in the dispersal distributions). Here we considered a decrease over time of the average squared axial parentoffspring distance (σ^{2}). Two different dispersal distributions with different σ^{2} values were used, while all other parameters of the distribution (i.e., the global shape of the distribution) remained unchanged. This situation corresponds to a change in a landscape (e.g., a fragmentation) resulting in modifying the ability of a species to move within this landscape (e.g., Brooker and Brooker 2002). Simulations were run with a twodimensional lattice of (500 × 500) nodes with one individual per node. A first dispersal distribution, given in expression (2) with parameters
Temporal change in density: A second category of fluctuations is temporal variations in density of individuals. We studied two simple situations: (i) a decrease in density from past to present (population bottleneck) and (ii) an increase in density from past to present (population flush). Such bottleneck or flush events are expected to occur in endangered or invasive populations, respectively. These situations were implemented in our simulations by changing the number of individuals per lattice node over time. Four different lattice models were used: one with 1 individual per node, one with 10 individuals per node, one with 1 individual every 3 nodes in each direction, and one with 1 individual every 10 nodes in each direction. These models correspond to densities of 1, 10, 1/9, and 1/100, respectively. Having less than 1 individual per node avoids the consideration of models with a too high number of individuals per node (i.e. >10) before or after a change in density, which would strongly deviate from the concept of continuous population to which the method of estimation applies. For easier coding, we modeled densities lower than 1 individual per node, considering that a given proportion of nodes of the lattice are always “empty” (e.g., for a density of 1/9, 8/9 of the nodes are empty). This is equivalent to a model with a larger lattice unit (e.g., a lattice unit three times larger in each dimension for a density of 1/9 compared to the lattice unit for a density of 1). A summary of the different density changes studied is presented in Table 1.
For the model with 1 individual every 9 nodes, we adapted the dispersal distribution to keep a constant σ^{2} = 4. Since dispersal may occur only between “nonempty” nodes, the dispersal distribution parameters are then M = 0.299 and n = 4.159 for 0 < k ≤ 48. For the model with 1/100, 1, or 10 individuals per node, the dispersal distribution parameters are those used in the previous section [cf. expression (5)]. We have not adapted the dispersal distribution to keep a constant σ^{2} = 4 for the model with 1 individual every 100 nodes because it was mathematically impossible to adjust this distribution with a too small number of points in the distribution (i.e., in this case, there are only five possible moves in each direction between “suitable” nodes, which are located at 0, 10, 20, 30, and 40 lattice units). However, additional simulations with a 90fold density increase (from 1/9 to 10 individuals per node) and a dispersal distribution adapted to keep a constant σ^{2} gave similar results (results not shown).
For each case of density change considered, four simulations were run, using a twodimensional habitat of (500 × 500) nodes with G_{c} = 10, G_{c} = 20, G_{c} = 100 generations, and G_{c} infinite as baseline. For each bottleneck and flush case, we simulated a weak density variation (10 and 9 times density change, respectively) and a strong density variation (90 and 100 times density change, respectively). In the case of bottleneck, the lowdensity models (1 and 1/9 individuals per node for weak and strong variations, respectively) were implemented from sampling time to G_{c} and the highdensity models (10 individuals per node) from G_{c} to the TMRCA. In the case of density flush, the highdensity models (1 individual per node) were implemented from sampling time to G_{c} and the lowdensity models (1/9 and 1/100 individuals per node for weak and strong variations, respectively) from G_{c} to the TMRCA (Table 1).
Spatial expansion with constant density: The third type of studied situation is a population expansion in space with constant density of individuals (Figure 1). The population introduced into an empty habitat is composed of individuals that have evolved in a source population at equilibrium with some demographic features (i.e., density and dispersal distribution). The introduced population spreads within a few generations on an empty twodimensional habitat with the same demographic features as the source population. This situation corresponds to the case of an introduced species that colonizes a new territory with similar ecological features to that of its native territory. Before expansion (i.e., at generation G_{c}), the introduced population is composed of 100 individuals located on a (10 × 10) area, which were sampled from a (10 × 10) area in the source population, which itself evolved on a (160 × 160) lattice. From generation G_{c} to present, the introduced population spreads over a lattice of (160 × 160) nodes. The entire (160 × 160) matrix is potentially occupied in two generations. At sampling time, as in the previous sections, 100 individuals were taken from an area of (20 × 20) nodes located outside the area of introduction, the distance between the introduction area and the sampling area being equal to 50 nodes. The forward dispersal distribution parameters are those given in expression (5) and correspond to a σ^{2} = 4. Four simulations were run with G_{c} = 10, G_{c} = 20, G_{c} = 100, and G_{c} infinite as baseline.
Spatial density heterogeneities: The situations we choose to study reflect the fact that biologists usually collect individual samples in localities where they are easy to collect, that is, in highdensity areas. Hence, we considered a lattice model with homogenous density except on a squared area where the density of individuals is higher (Figure 2). In such models with density heterogeneities in space, backward and forward dispersal differ. Each lattice node has a backward distribution that depends on the density of each surrounding node (e.g., each node being at a distance less or equal to the K_{max} step). Those surrounding nodes correspond to all locations from which genes could have come in one generation (forward in time). Since those nodes are occupied by different numbers of individuals and because nodes occupied by more individuals contribute potentially more to the number of immigrants that reach a given node, we have to weight each term of the backward dispersal distribution by the number of individuals of the node from where immigrants have come. Let N_{x}_{,}_{y}_{,}_{G} be the number of individuals at node (x, y) at generation G. Then for any node (x, y) the probability b_{dx}_{,}_{dy} for a gene to move backward dx steps in one direction and dy in the other is equal to
RESULTS
Interpretation of observed bias: Observed bias in our simulations might be attributable to (i) a bias, inherent to the method, due to the effect of a high mutation rate on the parameter value (this we call “mutational bias”), (ii) a bias due to the deviation of the estimates relative to the parameter value considering a finite sample of individuals and loci (this we name “small sample bias”), and (iii) a bias introduced by the demographic fluctuations studied. Additional details on the small sample and mutational biases can be found in Leblois et al. (2003). All results in the present study should be interpreted taking into account the small sample and mutational biases that can be observed in the simulations without demographic fluctuations that were included in all situations studied as baseline (G_{c} infinite). For example, in the case of a reduction of density (bottleneck, Table 3), the mutational and small sample bias is large when considering an intermediatedensity model (baseline simulation for a weak reduction) and much lower when considering a lowdensity model (baseline simulation for a stronger reduction). This difference is due partly to the different densities of individuals in the two baseline simulations, which influence the global level of genetic diversity in the sample. Leblois et al. (2003) indeed showed that differences in genetic diversity have a substantial effect on the estimation of Dσ^{2}.
Temporal change in dispersal: Simulation results show that the bias due to a reduction of dispersal is negative (Table 2) and thus corresponds to an overestimation of the present time Dσ^{2}. This result is in agreement with a transition from a high Dσ^{2} value (σ^{2} = 100) during the past generations (i.e., before G_{c}) to a much lower value after G_{c} (σ^{2} = 4). In other words, the method of Dσ^{2} estimation has a memory of temporal changes in dispersal. However, this memory is short term since a reduction of dispersal 100 generations ago gave only a slight negative bias compensated by the positive small sample and mutational biases (cf. first column of Table 2). Moreover, even for a recent reduction of dispersal (G_{c} = 10), the bias is <25% (i.e., <0.25), a relatively low value compared to the high amplitude of the dispersal change. Standard error of the estimation also remains low for all G_{c} values, and for changes older than 20 generations, >95% of the estimations are within a factor of two of the presenttime Dσ^{2}. Hence, our simulations generally show that the precision of the presenttime Dσ^{2} estimation is relatively robust to temporal changes in dispersal.
Temporal reduction of density (bottleneck): The negative bias observed in Table 3 (i.e., overestimation of Dσ^{2}) reflects the higher population density from generation G_{c} until the TMRCA. For a 10 times reduction of density, the method is quite robust when the density change occurred 20 or more generations ago. The bias and the MSE are low (<10%) and almost 99% of the estimations are within a factor of two of the presenttime Dσ^{2} value. For very recent density change (e.g., G_{c} = 10) the bias is substantial. However, the MSE remains low and >90% of the estimations are still within a factor of two of the presenttime Dσ^{2} value.
The effect of reduction of density is more marked for a stronger change in density (i.e., 90 times density reduction). For a very recent density reduction (i.e., 10 generations ago), the negative bias reaches 50% and only 24% of the estimations are within a factor of two of the presenttime Dσ^{2} value. For G_{c} = 100, the bias and the MSE become similar to the baseline. Note that all estimations are within a factor of two of the presenttime Dσ^{2} for G_{c} ≥ 20. Therefore, even for large recent density reductions, the method appears to be relatively robust.
Temporal increase in density (demographic flush): The positive bias observed in Table 4, which corresponds to an underestimation of the presenttime Dσ^{2}, reflects the lower population density from generation G_{c} until the TMRCA. For a small increase in density (10 times), the bias and the MSE are high even for a relatively ancient flush (e.g., G_{c} = 100). The proportion of estimations being within a factor of two of Dσ^{2} remains small (<50%) even for G_{c} = 100. The effect of the flush also increases substantially with the intensity of the density change. For a 100fold density change and for G_{c} = 10, the bias reaches 391% and none of the estimations are within a factor of two of Dσ^{2} (Table 4). Hence, although the bias and the MSE decrease when G_{c} increases, the estimation remains unreliable for both 100 and 10fold density change. These results contrast sharply with those pertaining to bottlenecks and dispersal changes.
Spatial increase in population size with constant density (demographic expansion): All measures (bias, MSE, and proportion of estimates within a factor of two) indicate that the estimation of the presenttime Dσ^{2} is good when the spatial expansion occurred 20 or more generations ago (Table 5). For G_{c} = 10 only, an 8% negative bias is observed, which corresponds to an overestimation of the presenttime Dσ^{2} (Table 5). However, the MSE is very small (10%) and 97% of the estimations are within a factor of two of the expected Dσ^{2} value. Hence, a spatial expansion as modeled here has only a shortterm and limited influence on the presenttime Dσ^{2} estimation; the method is precise even for very recent expansions.
Spatial heterogeneity in density (sampling within a highdensity zone): Table 6 shows that Dσ^{2} estimation is not robust when the highdensity zone is small and strictly corresponds to the sampling area. The bias and MSE values indicate that in this case the lowdensity area surrounding the sampling area strongly influences the Dσ^{2} estimation, which becomes a bad measure of both local density (i.e., the density on the sampling area) and surrounding density (i.e., the density surrounding the sampling area). It can be seen, however, that two times coverage probabilities, although globally low, are higher when referring to the local rather than to the surrounding area density as expected (Dσ^{2} value 0.018 vs. 0.001). This suggests that there is a tendency for the method to measure the local rather than the surrounding density. This trend becomes obvious when looking at results for a larger highdensity zone (Table 6). In this case, the bias and the MSE are much lower when considering the local rather than the surrounding zone for the Dσ^{2} value. About 90% of the estimates are within a factor of two of the local Dσ^{2} value, while none of them are within a factor of two of the surrounding Dσ^{2} value. The third case of a large highdensity zone located outside the sampling area (i.e., 50 nodes away) confirms this result (Table 6). Hence, our simulations generally show that the method estimates local demographic parameters and is robust for such measurement when the density is relatively homogenous around the sampling area (e.g., over an area equal to four times the sampling area).
DISCUSSION
This work is the first one focusing on the study of evolutionary disequilibrium situations in the complex but realistic population model of a continuous population evolving under isolation by distance. Within the limits of the situations studied in this article, and with the exception of the case of a density flush, we found that temporal and spatial fluctuations of demographic parameters, if not too strong and not too recent (i.e., more than, say, 2050 generation in the past), have a limited influence on the estimation of local and presenttime demographic parameters with the method of Rousset (2000). It is worth noting that we are talking about changes on timescales of a few tens of generations in the past, which may be very recent by standards in population genetics, but not for lots of species undergoing demographic changes due to ongoing human impact. Moreover, the numbers of generations defining the time of demographic change in this study should be considered as indicative of only the length of the effect of the demographic changes studied rather than as absolute reference numbers. As a matter of fact, the persistence in time of the effect of demographic fluctuations strongly depends on various features of the demographic model (e.g., σ^{2} values) and disequilibrium situations. It is thus preferable to consider general trends rather than precise numbers for each situation. For clarity, those trends have been summarized in Table 7.
The robustness of the method of Rousset (2000) to several temporal and spatial demographic fluctuations somewhat contradicts previous studies dealing with the study of evolutionary disequilibrium. In their review, Koenig et al. (1996) concluded that estimations of dispersal parameters from genetic data give ideas about past rather than present dispersal and gene flow, so that direct methods, such as markrecapture methods, should give a better estimation of actual dispersal parameters. Boileau et al. (1992) similarly showed that hundreds or thousands of generations are required to erase the effects of colonization processes on “F_{ST}like estimates” from allozyme data in large populations, concluding that estimates of gene flow from genetic data should be taken with care. We fully agree that temporal demographic fluctuations in a population are likely to have a strong and persistent effect on some population genetics statistics and methods. However, the present study shows that some indirect methods and genetic markers give accurate estimations of presenttime density and dispersal features even when the demographic history includes relatively recent demographic changes.
The general robustness to spatial and temporal heterogeneities of the present Fstatisticbased method can be interpreted using arguments from the coalescence theory and analytical treatment available in this field. Values of Fstatistics, under the assumption of low mutation rate, can be deduced by comparing the distributions of coalescence probability for different pairs of genes (e.g., pairs from the same deme and pairs from different demes; e.g., Rousset 2002). These distributions differ essentially by an excess of coalescence probability for the most related genes, this excess being concentrated in a brief period τ in the recent past. Fstatistics thus depend mainly on differences between the distributions of coalescence probability for different pairs of genes in recent generations. As the sensitivity of Fstatistics values to past demographic fluctuations is also related to this recent time period, past demographic fluctuations have less effect when the time period τ is short. This recent time period τ is shorter when high dispersal rates and/or low deme size are considered (Rousset 2004). Hence, if models with small deme size and high migration rates, such as isolation by distance between individuals where each deme is of size two genes, are considered the influence of past demographic fluctuations on the estimation of demographic parameters from Fstatistics is limited. By contrast, under the classical island model with large deme size and low migration rates, the effect of past demographic fluctuations is expected to be more problematic. Moreover, under isolationbydistance models, the more distant the demes are on the lattice, the more the period τ is expanding to the past, increasing the effect of past demographic parameter fluctuations (Slatkin 1994; Rousset 2004). Because the present method focuses on local differentiation and thus on recent evolutionary processes corresponding to a narrow recent past zone, it is again logical that past demographic fluctuations have limited effects on the estimation of the presenttime and local Dσ^{2} with this method. The same reasoning can be used to understand why the method gives estimates of the local demographic parameter values rather than estimates of the surrounding demographic parameter values. As the period τ is short in the models considered, Fstatistics depend mainly on genetic events (migration, coalescence, mutation) that occurred in a recent past and, because dispersal is localized, at a local geographical scale. Therefore, the estimate of Dσ^{2} by the present method should correspond to the local demographic parameter values on the sampling area and should not be much influenced by demographic features of zones that are far away from the sampling area.
Close examination of our results brings up several issues. Our simulations showed that, for the study of invading species, the present method should give precise estimates of the presenttime Dσ^{2} provided that no demographic flush occurred during the expansion process. This is an interesting feature of the method, which makes it appropriate to study invasive organisms for which demographic features are similar in the newly founded population and in the original source population. Our simulations further showed that if a change in dispersal occurred during the invasion process, this new dispersal feature should translate quickly in the estimation of the presenttime Dσ^{2}. On the other hand, density flushes (and to a much lower extent population bottlenecks) may strongly affect presenttime Dσ^{2} estimation. Invading species populations often experience complex demographic fluctuations that may include both bottlenecks (i.e., founder events) and density flushes during their spreading (e.g., Williamson 1996; Estoupet al. 2001). Therefore, it seems necessary to run additional simulations adapted to those complex demographic scenarios to thoroughly evaluate the robustness of the estimation of the presenttime Dσ^{2}.
Our simulations also show that for conservation biology studies dealing with bottlenecked populations the estimation of Dσ^{2} is potentially biased toward past demographic parameter values. However, the memory of past demographic parameter values is short so that this bias is important for only a strong and recent decrease in density. A major genetic consequence of a population bottleneck is that the number of alleles decreases much faster than the heterozygosity (Neiet al. 1975; Luikart and Cornuet 1998). One might have expected the precision of the method to be reduced due to the lower number of alleles in the bottlenecked population. However, standard error on the bias was weak whatever the strength of the bottleneck. One possible explanation for this result is that the method’s precision depends more on the mean heterozygosity level than on the average number of alleles.
Our simulations indicate that surrounding densities considerably influence the estimation of local Dσ^{2} when the sample is taken on a small highdensity zone. In this case, the estimates correspond neither to the Dσ^{2} values on the sampling area nor to the surrounding Dσ^{2} values. However, if sampling is done in a sufficiently large highdensity zone (e.g., on a surface equals to four times the sampling area), the estimates correspond more to the local density (i.e., the density in the sampling area). Our simulations allowed us to study the case of a highdensity zone in the middle of a large homogenous zone with low density. This situation is realistic for various demographic systems and mimics a classical experimental bias (i.e., the fact that biologists generally collect their samples in highdensity areas). However, many biological situations with spatial density heterogeneities would correspond rather to random density fluctuations on each lattice node. It is expected that differentiation in such scenarios will be a function of some “effective” density and dispersal rate. The lack of analytical formulas for these effective parameters limits the interpretation of a simulation study of the performance of estimators. Nevertheless, there is no obvious reason to believe that the estimation of the effective Dσ^{2} would be affected more by such random fluctuations than by previously studied spatial heterogeneities.
Acknowledgments
We thank Thomas Lenormand and Franck Shaw for constructive comments on the manuscript. This work was financially supported by the Action Incitative Programmée no. 00202 “biodiversité” from the Institut Français de la Biodiversité and grant no. D4E/SRP/01118 “biological invasion” from the Ministère de l’Ecologie et du Développement Durable. This is paper ISEM 2004007.
Footnotes

Communicating editor: L. Excoffier
 Received July 4, 2003.
 Accepted October 18, 2003.
 Copyright © 2004 by the Genetics Society of America