Influence of Spatial and Temporal Heterogeneities on the Estimation of Demographic Parameters in a Continuous Population Using Individual Microsatellite Data

Corresponding author: Raphael Leblois, Institut des Sciences de l'Evolution, UMR 5554-CC065, Université des Sciences et Techniques du Languedoc, Pl. E. Bataillon, 34095 Montpellier, France., leblois{at}isem.univ-montp2.fr (E-mail)

Communicating editor: L. EXCOFFIER

	ABSTRACT

TOP ABSTRACT MODELS AND METHODS RESULTS DISCUSSION LITERATURE CITED

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² from individual microsatellite data, where D is the density of adults and ² the average squared axial parent-offspring distance in a continuous population evolving under isolation by distance. A coalescence-based simulation algorithm was used to study the effect on D² estimation of temporal and spatial fluctuations of demographic parameters. Estimation of present-time D² 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² 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 present-time 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., mark-recapture 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 ; KOENIG et 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 present-time population dynamics of a species. WHITLOCK and MCCAULEY 1999 recently concluded that estimates of the number of migrants between subpopulations from F-statistics 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 present-time 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 Wright-Fisher 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 LEBLOIS et al. 2003 ). A method of analysis adapted to a "continuous" population evolving under isolation by distance was developed to estimate the product D², where D is the density of adults and ² the average squared axial parent-offspring 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 4D² (ROUSSET 1997 ). The method is valid for leptokurtic distributions of dispersal distance (ROUSSET 2000 ; LEBLOIS et 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² (ROUSSET 2000 ; SUMNER et al. 2001 ).

As for any population genetics method of demographic parameter estimation, the quality of the estimation of D² using this method may be affected by local and temporal spatial heterogeneities in demographic parameters. In this study, we adapted the coalescence-based simulation algorithm of LEBLOIS et al. 2003 to study the effect of temporal and spatial fluctuations of demographic parameters on the estimation of present-time D². 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 present-time D² 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 high-density zone within a lower-density area.

	MODELS AND METHODS

TOP ABSTRACT MODELS AND METHODS RESULTS DISCUSSION LITERATURE CITED

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 (MALECOT 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., MALECOT 1967 ; SAWYER 1977 ), it is known that they do not follow a well-defined set of biological assumptions (MARUYAMA 1972 ; FELSENSTEIN 1975 ; see BARTON et al. 2002 for an alternative approach for continuous populations). To avoid edge effects, a two-dimensional 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 (LEBLOIS et 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 large-N approximation of the n-coalescent 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 generation-by-generation algorithm leads to less efficient simulations in terms of computation time than do those based on the n-coalescent theory (KINGMAN 1982A , KINGMAN 1982B ; 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 generation-by-generation 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., MALECOT 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 parent-offspring axial distance, backward in time, expressed in number of steps on the lattice. Under a two-dimensional 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 ; KOT et 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

(1)

with parameters M and n controlling the total dispersal rate and the kurtosis, respectively. This distribution corresponds to a truncated variant of the discrete Pareto, or , distribution (see, e.g., PATIL and JOSHI 1968 ). By suitable choice of the two parameter values, large kurtosis can be obtained with high migration rates (ROUSSET 2000 ). For some distributions, the first p terms were arbitrarily fixed:

(2)

Dispersal was assumed to be independent in each direction, so that f_dx,dy = f_dx x 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 x 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² estimation. Therefore, microsatellite markers were simulated in the present study. On the basis of direct observations of mutations at human microsatellite loci (DIB et 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 (PRITCHARD et al. 1999 ; ESTOUP et al. 2001 ). The variance of the geometric distribution was fixed at 0.36 (ESTOUP et al. 2001 ), a value computed from the mutation data in DIB et al. 1996 . The mutation rate was equal to 5 x 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 ; SCHLOTTERER 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² estimations.

Method of analysis:
Each simulation iteration gives the genotypes at 10 polymorphic loci of 100 (i.e., 10 x 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 x 10) individuals taken every two nodes from an area of (20 x 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

(3)

where SS_{b[etween]() i,j}(X_i.:u - X_i.:u)² measures divergence between genes taken from two different individuals and SS_{w[ithin]() i,j,u}(X_ij:u - X_i.:u)² measures divergence between genes within the same individual (X_ij:u is an indicator variable taking the value 1 if gene i of individual j is of allelic type u and the value 0 otherwise; ROUSSET 2000 ). Thus, â compares the genetic divergence of individuals at distance r (numerator) to the divergence of the two-gene copy within the individual (denominator), which is essentially what the parameter a_r does. Because stepwise mutations occur at microsatellite loci, a statistic taking into account the allele size might appear to be attractive. However, LEBLOIS et al. 2003 have shown that incorporation of allele size into the estimate of a_r gives unreliable results due to the high variance of the estimates. Therefore, only the parameter a_r described in Equation 3 was used in this study.

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² is related to the parameter

(4)

where Q_w is the probability of identity of genes within individuals, Q_r is the probability of identity of two genes in different individuals at distance r, and Q₀ is the probability of identity of two genes in different individuals in the same node (ROUSSET 2004 , Equation 8.12). Without random selfing, â_r is not the most relevant statistic. Rather one should estimate not only Q_w but also Q₀. Since there is only one adult per node of the lattice, Q₀ cannot be estimated directly from adults: it must be approximated as the identity between close adults or (better) between close juveniles before competition (see ROUSSET 2004 , Chap. 8, for further discussion). In this way, it is easy to adapt the methods considered in this article, but this is not considered further.

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 4D² (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, r << /. Beyond this limit, the linear relationship between a_r and the logarithm of the distance holds less well (see ROUSSET 1997 for theoretical details). Thus, if the analysis is done on a small geographical scale, the use of loci with high mutation rates such as microsatellites does not bias the estimation. This is illustrated by LEBLOIS et al. 2003 , using simulations.

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². Therefore, we present the bias and the MSE for the slope values of the regression lines and not for D² 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/(4D²) [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)²]. 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/(4D²) was within a factor of two from the expected value (i.e., in the interval [expected slope/2; 2 x expected slope]).

Spatial and temporal heterogeneities:
One important advantage of the generation-by-generation 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 isolation-by-distance 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 parent-offspring distance (²). Two different dispersal distributions with different ² 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 two-dimensional lattice of (500 x 500) nodes with one individual per node. A first dispersal distribution, given in expression (2) with parameters

(5)

has a moderate ² value (² = 4 in lattice units) and is the dispersal distribution from the present until the time of change, G_c. A second dispersal distribution, with parameters M = 0.187 and n = 1.246 for 0 < k 48 corresponds to a very high ² value (² = 100) and is the dispersal distribution from the time of change G_c until the time of the most recent common ancestor (TMRCA). Four simulations were run with G_c = 10, G_c = 20, G_c = 100 generations (going backward in time), and G_c infinite as baseline (i.e., no change in dispersal features over time).

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.

View this table: In this window In a new window   Table 1. Models used to study the effects of density variation in time on the estimation of 1/(4D2)

For the model with 1 individual every 9 nodes, we adapted the dispersal distribution to keep a constant ² = 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 ² = 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 90-fold density increase (from 1/9 to 10 individuals per node) and a dispersal distribution adapted to keep a constant ² gave similar results (results not shown).

For each case of density change considered, four simulations were run, using a two-dimensional habitat of (500 x 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 low-density models (1 and 1/9 individuals per node for weak and strong variations, respectively) were implemented from sampling time to G_c and the high-density models (10 individuals per node) from G_c to the TMRCA. In the case of density flush, the high-density models (1 individual per node) were implemented from sampling time to G_c and the low-density 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 (Fig 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 two-dimensional 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 x 10) area, which were sampled from a (10 x 10) area in the source population, which itself evolved on a (160 x 160) lattice. From generation G_c to present, the introduced population spreads over a lattice of (160 x 160) nodes. The entire (160 x 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 x 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 ² = 4. Four simulations were run with G_c = 10, G_c = 20, G_c = 100, and G_c infinite as baseline.

View larger version (26K): In this window In a new window Download PPT slide   Figure 1. Schema of a demographic expansion with constant density as modeled in this study. (a) The source population from which a subpopulation (dark gray grid) is introduced in an empty habitat (dotted arrow). (b) The empty habitat on which the introduced population spreads within a few generations (solid arrows). In our simulations, two-dimensional habitats are represented on a torus and not on a plane square as in this figure.

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 high-density areas. Hence, we considered a lattice model with homogenous density except on a squared area where the density of individuals is higher (Fig 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

(6)

View larger version (18K): In this window In a new window Download PPT slide   Figure 2. Schema of the spatial density heterogeneities as modeled in this study. (a) A small high-density zone (dark gray grid) strictly corresponds to the sampling area (black grid) on a two-dimensional habitat with a lower density (light gray grid). (b) A large high-density zone (dark gray grid) includes the sampling area (blackgrid) on a two-dimensional habitat with a lower density (light gray grid). (c) A large high-density zone (dark gray grid) is present on a two-dimensional habitat with a lower density (light gray grid); the sampling area (black grid) is located outside the high-density zone. In our simulations, two-dimensional habitats are represented on a torus and not on a plane square as in this figure.

Simulations were run for a sample of 100 individuals taken every two nodes from an area of (20 x 20) nodes evolving in a (160 x 160) lattice. Density is one individual per node, except on a (n x n) zone including the sample area where density is 10 individuals per node. Two cases were considered: (i) a small high-density zone of (20 x 20) nodes, which strictly corresponds to the sample area (Fig 2A), and (ii) a larger high-density zone of (40 x 40) nodes, which includes the (20 x 20) nodes sample area (Fig 2B). We were particularly interested in assessing whether the estimated density corresponds to the density on the sampling area (i.e., the local density) or whether the estimation is influenced largely by the density surrounding the sampling area (i.e., the neighboring density). This was performed by alternatively considering that the expected D² value corresponded to a density of 10 (local density) and 1 (surrounding density) individuals per node. An additional simulation was run with a single large high-density zone of (40 x 40) nodes located outside the sampling area, the distance between the high-density and sampling zones being equal to 50 nodes (Fig 2C).

	RESULTS

TOP ABSTRACT MODELS AND METHODS RESULTS DISCUSSION LITERATURE CITED

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 intermediate-density model (baseline simulation for a weak reduction) and much lower when considering a low-density 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².

View this table: In this window In a new window   Table 2. Effect of a temporal reduction of dispersal on the estimation of 1/(4D2)

View this table: In this window In a new window   Table 3. Effect of a weak (10 times density reduction) and strong (90 times density reduction) bottleneck on the estimation of 1/(4D2)

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². This result is in agreement with a transition from a high D² value (² = 100) during the past generations (i.e., before G_c) to a much lower value after G_c (² = 4). In other words, the method of D² 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 present-time D². Hence, our simulations generally show that the precision of the present-time D² 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²) 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 present-time D² 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 present-time D² 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 present-time D² 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 present-time D² 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 present-time D², 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² 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 100-fold density change and for G_c = 10, the bias reaches 391% and none of the estimations are within a factor of two of D² (Table 4). Hence, although the bias and the MSE decrease when G_c increases, the estimation remains unreliable for both 100- and 10-fold density change. These results contrast sharply with those pertaining to bottlenecks and dispersal changes.

View this table: In this window In a new window   Table 4. Effect of a weak (9 times density increase) and strong (100 times density increase) density flush on the estimation of 1/(4D2)

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 present-time D² 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 present-time D² (Table 5). However, the MSE is very small (10%) and 97% of the estimations are within a factor of two of the expected D² value. Hence, a spatial expansion as modeled here has only a short-term and limited influence on the present-time D² estimation; the method is precise even for very recent expansions.

Spatial heterogeneity in density (sampling within a high-density zone):
Table 6 shows that D² estimation is not robust when the high-density zone is small and strictly corresponds to the sampling area. The bias and MSE values indicate that in this case the low-density area surrounding the sampling area strongly influences the D² 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² 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 high-density 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² value. About 90% of the estimates are within a factor of two of the local D² value, while none of them are within a factor of two of the surrounding D² value. The third case of a large high-density 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

TOP ABSTRACT MODELS AND METHODS RESULTS DISCUSSION LITERATURE CITED

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, 20–50 generation in the past), have a limited influence on the estimation of local and present-time 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., ² 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 mark-recapture 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 present-time 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 F-statistic-based method can be interpreted using arguments from the coalescence theory and analytical treatment available in this field. Values of F-statistics, 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. F-statistics thus depend mainly on differences between the distributions of coalescence probability for different pairs of genes in recent generations. As the sensitivity of F-statistics 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 F-statistics 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 isolation-by-distance 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 present-time and local D² 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, F-statistics 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² 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 present-time D² 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 present-time D². On the other hand, density flushes (and to a much lower extent population bottlenecks) may strongly affect present-time D² 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 ; ESTOUP et 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 present-time D².

Our simulations also show that for conservation biology studies dealing with bottlenecked populations the estimation of D² 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 (NEI et 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² when the sample is taken on a small high-density zone. In this case, the estimates correspond neither to the D² values on the sampling area nor to the surrounding D² values. However, if sampling is done in a sufficiently large high-density 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 high-density 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 high-density 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² 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 2004-007.

Manuscript received July 4, 2003; Accepted for publication October 18, 2003.

	LITERATURE CITED

TOP ABSTRACT MODELS AND METHODS RESULTS DISCUSSION LITERATURE CITED

BARTON, N. H., F. DEPAULIS, and A. M. ETHERIDGE, 2002  Neutral evolution in spatially continuous populations. Theor. Popul. Biol. 61:31-48.[CrossRef][Medline]

BEEBEE, T. and G. ROWE, 2001  Application of genetic bottleneck testing to the investigation of amphibian declines: a case study with natterjack toads. Conserv. Biol. 15:266-270.[CrossRef]

BOILEAU, M. G., P. D. N. HEBERT, and S. S. SCHWARTZ, 1992  Non-equilibrium gene frequency divergence: persistent founder effects in natural populations. J. Evol. Biol. 5:25-39.[CrossRef]

BROOKER, L. and M. BROOKER, 2002  Dispersal and population dynamics of the blue-breasted fairy-wren, Malurus pulcherrimus, in fragmented habitat in the Western Australian wheatbelt. Wildlife Res. 29:225-233.[CrossRef]

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

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

ENDLER, J. A., 1977 Geographical variation, speciation, and clines. Princeton University Press, Princeton, NJ.

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

ESTOUP, A., I. J. WILSON, C. SULLIVAN, J.-M. CORNUET, and C. MORITZ, 2001  Inferring population history from microsatellite and enzyme data in serially introduced cane toads, Bufo marinus.. Genetics 159:1671-1687.[Abstract/Free Full Text]

FELSENSTEIN, J., 1975  A pain in the torus: some difficulties with models of isolation by distance. Am. Nat. 109:359-368.[CrossRef]

HASTINGS, A. and S. HARRISON, 1994  Metapopulation dynamics and genetics. Annu. Rev. Ecol. Syst. 25:167-188.[CrossRef]

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

KINGMAN, J. F. C., 1982b  On the genealogy of large populations. J. Appl. Probab. 19A:27-43.[CrossRef]

KOENIG, W. D., D. VAN VUREN, and P. N. HOOGE, 1996  Detectability, philopatry, and the distribution of dispersal distances in vertebrates. Trends Ecol. Evol. 11:514-517.[CrossRef]

KOT, M., M. A. LEWIS, and P. VAN DEN DRIESSCHE, 1996  Dispersal data and the spread of invading organisms. Ecology 77:2027-2042.[CrossRef]

LEBLOIS, R., A. ESTOUP, and F. ROUSSET, 2003  Influence of mutational and sampling factors on the estimation of demographic parameters in a ‘continuous’ population under isolation by distance. Mol. Biol. Evol. 20:491-502.[Abstract/Free Full Text]

LUIKART, G. and J.-M. CORNUET, 1998  Empirical evaluation of a test for identifying recently bottlenecked populations from allele frequency data. Conserv. Biol. 12:228-237.[CrossRef]

MALÉCOT, G., 1967 Identical loci and relationship, pp. 317–332 in Proceedings of the Fifth Berkeley Symposium on Mathematical statistics and Probability, Vol. 4, edited by L. M. LECAM and J. NEYMAN. University of California Press, Berkeley, CA.

MALÉCOT, G., 1975  Heterozygosity and relationship in regularly subdivided populations. Theor. Popul. Biol. 8:212-241.[CrossRef][Medline]

MARUYAMA, T., 1972  Rate of decrease of genetic variability in a two-dimensional continuous population of finite size. Genetics 70:639-651.[Abstract/Free Full Text]

NEI, M., T. MARUYAMA, and R. CHAKRABORTY, 1975  The bottleneck effect and genetic variability in populations. Evolution 29:1-10.

NORDBORG, M., 2001 Coalescent theory, pp. 179–208 in Handbook of Statistical Genetics, edited by D. A. BALDING, M. BISHOP and C. CANNINGS. John Wiley & Sons, Chichester, UK.

PATIL, G. P., and S. W. JOSHI, 1968 A Dictionary and Bibliography of Discrete Distribution. Oliver & Boyd, Edinburgh.

PORTNOY, S. and M. F. WILLSON, 1993  Seed dispersal curves: behavior of the tail of the distribution. Evol. Ecol. 7:25-44.

PRITCHARD, J. K., M. T. SEIELSTAD, A. PEREZ-LEZAUN, and M. W. FELDMAN, 1999  Population growth of human Y chromosome microsatellites. Mol. Biol. Evol. 16:1791-1798.[Abstract]

ROUSSET, F., 1996  Equilibrium values of measures of population subdivision for stepwise mutation processes. Genetics 142:1357-1362.[Abstract]

ROUSSET, F., 1997  Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145:1219-1228.[Abstract]

ROUSSET, F., 2000  Genetic differentiation between individuals. J. Evol. Biol. 13:58-62.[CrossRef]

ROUSSET, F., 2001 Inferences from spatial population genetics, pp. 239–265 in Handbook of Statistical Genetics, edited by D. A. BALDING, M. BISHOP and C. CANNINGS. John Wiley & Sons, Chichester, UK.

ROUSSET, F., 2002  Inbreeding and relatedness coefficients: What do they measure? Heredity 88:371-380.[CrossRef][Medline]

ROUSSET, F., 2004 Genetic Structure and Selection in Subdivided Populations. Princeton University Press, Princeton, NJ.

SAWYER, S., 1977  Asymptotic properties of the equilibrium probability of identity in a geographically structured population. Adv. Appl. Probab. 9:268-282.[CrossRef]

SCHLÖTTERER, C., 2000  Evolutionary dynamics of microsatellite DNA. Chromosoma 109:365-371.[Medline]

SLATKIN, M., 1993  Isolation by distance in equilibrium and non-equilibrium populations. Evolution 47:264-279.[CrossRef]

SLATKIN, M., 1994 Gene flow and population structure, pp. 3–17 in Ecological Genetics, edited by L. A. REAL. Princeton University Press, Princeton, NJ.

SPONG, G. and L. HELLBORG, 2002  A near-extinction event in lynx: Do microsatellite data tell the tale? Conserv. Ecol. 6(1):15.

STONE, G. N. and P. SUNNUCKS, 1993  Genetic consequences of an invasion through a patchy environment: the cynipid gallwasp Andricus quercuscalicis (Hymenoptera: Cynipidae). Mol. Ecol. 2:251-268.

SUMNER, J., F. ROUSSET, A. ESTOUP, and C. MORITZ, 2001  ‘Neighborhood’ size, dispersal and density estimates in the prickly forest skink (Gnypetoscincus queenslandiae) using individual genetic and demographic methods. Mol. Ecol. 10:1917-1927.[CrossRef][Medline]

WILLIAMSON, M., 1996 Biological Invasions. Chapman & Hall, London.

WHITLOCK, M. C. and D. E. MCCAULEY, 1999  Indirect measure of gene flow and migration: F_ST1/(4Nm+1). Heredity 82:117-125.

This article has been cited by other articles:

				F. Rousset and R. Leblois Likelihood and Approximate Likelihood Analyses of Genetic Structure in a Linear Habitat: Performance and Robustness to Model Mis-Specification Mol. Biol. Evol., December 1, 2007; 24(12): 2730 - 2745. [Abstract] [Full Text] [PDF]