## Abstract

Reconstructing the demographic history of populations is a central issue in evolutionary biology. Using likelihood-based methods coupled with Monte Carlo simulations, it is now possible to reconstruct past changes in population size from genetic data. Using simulated data sets under various demographic scenarios, we evaluate the statistical performance of Msvar, a full-likelihood Bayesian method that infers past demographic change from microsatellite data. Our simulation tests show that Msvar is very efficient at detecting population declines and expansions, provided the event is neither too weak nor too recent. We further show that Msvar outperforms two moment-based methods (the *M*-ratio test and Bottleneck) for detecting population size changes, whatever the time and the severity of the event. The same trend emerges from a compilation of empirical studies. The latest version of Msvar provides estimates of the current and the ancestral population size and the time since the population started changing in size. We show that, in the absence of prior knowledge, Msvar provides little information on the mutation rate, which results in biased estimates and/or wide credibility intervals for each of the demographic parameters. However, scaling the population size parameters with the mutation rate and scaling the time with current population size, as coalescent theory requires, significantly improves the quality of the estimates for contraction but not for expansion scenarios. Finally, our results suggest that Msvar is robust to moderate departures from a strict stepwise mutation model.

INFERRING past demography is a central concern in evolutionary biology and applied ecology. Characterizing past variations in population size is crucial, *e.g.*, for understanding the impact of past climatic fluctuations on the current distribution of species (Jacobsen *et al.* 2005; Elmer *et al.* 2009; Hu *et al.* 2009) and for the conservation of endangered species (Frankham *et al.* 2002). Characterizing the demographic history of a species by direct approaches requires the monitoring of census data, which can be extremely difficult, not to say impossible, particularly in long-lived species. Yet variations in census numbers of individuals also affect the dynamics of the genes carried by these individuals. A powerful alternative to direct approaches is therefore to use the recent advances in population genetic theory, which allow inferences on past demography from the observed distribution of genetic variation in natural populations (Lawton-Rauh 2008).

Until recently, most of these indirect methods relied on summary statistics calculated from genetic data and tests for departure from their theoretical distribution under a given demographic and mutational model (Cornuet and Luikart 1996; Schneider and Excoffier 1999; Garza and Williamson 2001). For instance, Cornuet and Luikart’s (1996) approach relies on the rationale that rare alleles, which contribute marginally to the heterozygosity, are more likely to be lost following a bottleneck. A transient excess in heterozygosity, compared to that expected at equilibrium given the observed number of alleles in the sample, can therefore be used as a proxy to detect a bottleneck (Luikart and Cornuet 1998). Conversely, a transient heterozygosity deficiency may provide evidence for a population expansion (Cornuet and Luikart 1996; Leblois *et al.* 2006). In the same line of ideas, Garza and Williamson (2001) proposed a test to detect past population declines, based on the ratio (*M*) of the number of alleles to the range in allele size observed at microsatellite loci. Because they are easy to implement and do not require time-consuming computations, these moment-based methods have been used in many empirical studies (see, *e.g.*, Spencer *et al.* 2000; Comps *et al.* 2001; Colautti *et al.* 2005). However, these methods suffer from a limited statistical power because they do not make full use of the data. Furthermore, they do not provide any estimate of the severity and the duration of the bottleneck.

Likelihood-based methods coupled with Monte Carlo sampling offer a powerful alternative to these moment-based methods (Felsenstein 1992; Griffiths and Tavaré 1994; Emerson *et al.* 2001). They rely upon the computation of the likelihood of a sample configuration, *i.e.*, the probability to observe the allele counts or the DNA polymorphic sites in that sample, given a demographic and mutational model. The parameters of interest of the underlying model are then estimated by maximizing the likelihood of the observed data. Likelihood-based methods that have been developed for inferring past demographic changes from the observed distribution of genetic variation include, *e.g.*, Batwing (Wilson and Balding 1998; Wilson *et al.* 2003), Beast (Drummond and Rambaut 2007), IM and IMa (Hey and Nielsen 2004, 2007), Lamarc (Kuhner 2006), and Msvar (Beaumont 1999). These methods differ not only with respect to the underlying demographic model, but also with respect to the markers used (microsatellites, DNA sequences, etc.). However, because the computational burden required to evaluate statistical power and accuracy is particularly high, only few studies have attempted to test these methods (Wilson *et al.* 2003; Abdo *et al.* 2004; Rousset and Leblois 2007; Chikhi *et al.* 2010; Strasburg and Rieseberg 2010).

Among those methods, the one developed by Beaumont (1999), implemented in the software package Msvar and further improved by Storz and Beaumont (2002) and Storz *et al.* (2002), has been increasingly used in the past few years to infer past demographic changes (supporting information, Table S1). Msvar assumes a demographic model consisting of a single isolated population, which has undergone a linear or exponential change in effective population size at some time in the past. This method is designed to analyze multilocus microsatellite data that evolve according to a stepwise mutation model (SMM) (Ellegren 2004). Msvar uses a Markov chain Monte Carlo (MCMC) method to sample from the posterior distribution of the model parameters (*i.e.*, the current effective population size, the ancestral effective population size before the demographic change, the time at which the latter occurred, and the mutation rate of microsatellite loci).

Although Msvar has been widely used, the statistical performance of the method has never been extensively evaluated. In his original article, Beaumont (1999) simulated a handful of data sets with known mutational and demographic parameters and then evaluated the performance of the method for detecting demographic events and its sensitivity to the shape (linear or exponential) of the demographic change. However, the precision of the estimation of the model parameters was not evaluated. Furthermore, the performance of Msvar with respect to the severity of demographic change, the time since the population started changing in size, and the mutation model has not been studied yet.

Here, we therefore aimed at evaluating the statistical performance of Msvar (i) in detecting population declines and expansions and (ii) in providing accurate estimates of the model parameters, as a function of the severity of the demographic change, the time since it occurred, and the mutation model. To that end, we performed stochastic simulations to generate microsatellite data sets under different demographic scenarios and mutation models and then analyzed these simulated data with Msvar. In light of our results, we comment upon the published empirical studies that used Msvar and provide some guidelines for future studies.

## METHODS

#### Demographic model:

The demographic model implemented in Msvar (Beaumont 1999; Storz and Beaumont 2002) considers an isolated panmictic population of size *N*_{0} at sampling time (*t* = 0). Going backward in time, the population size *N*(*t*) changes deterministically (either linearly or exponentially) to an ancestral size *N*_{1} at time *t* = *T*_{a} and then remains constant at *N*_{1} for *t* > *T*_{a} (Beaumont 1999). In the following, we will consider only an exponential change in population size, with*t* < *T*_{a}, and *N*(*t* ≥ *T*_{a}) = *N*_{1}. For simplicity, the time is measured in units of generations, and population sizes are expressed as numbers of diploid individuals.

#### Simulation study:

To test how Msvar performed depending upon the nature of the demographic change (decline or expansion), its strength, and its time of occurrence, we simulated population declines and expansions for a range of parameter values for the current population size *N*_{0}, the ancestral population size *N*_{1}, and the time *T*_{a}. The computational burden of the method prevented an exhaustive exploration of the parameter space. In a first set of simulations, we therefore concentrated on a set of parameter values that represented a range of situations characterized by weak, moderate, and strong changes in population size, with varying time of occurrence. For population declines, we fixed the current population size *N*_{0} = 100 in all scenarios and varied the ancestral population size *N*_{1} = {1000; 10,000; 100,000} and the time since the demographic change *T*_{a} = {10; 50; 100; 500}. For population expansions, we fixed the ancestral population size *N*_{1} = 100 in all scenarios and varied the current population size *N*_{0} = {1000; 10,000; 100,000} and the time since the demographic change *T*_{a} = {10; 50; 100; 500}. A total of 24 sets of demographic parameters were therefore considered. For this first set of simulations, each locus evolved according to a strict SMM, as assumed in Msvar. The mutation rate μ was set at 10^{−3}, which is in agreement with estimates from the literature (Ellegren 2004).

Then, to test how Msvar performed depending upon the mutation model, we performed a second set of simulations. The mutation process of microsatellites is complex and highly heterogeneous across loci and organisms (Ellegren 2000, 2004). While some observations of spontaneous mutations support a strict SMM, others suggest that multistep mutations occur, with a frequency of multistep changes *p* varying from 0.04 to 0.74 (Ellegren 2000, 2004). Apart from a strict SMM, we thus simulated microsatellite data under a generalized stepwise model (GSM) with *p* = 0.22, an average value found in the literature (Dib *et al.* 1996; Ellegren 2000; Estoup *et al.* 2001; Ellegren 2004), and with *p* = 0.74, the most extreme value reported ever (Fitzsimmons 1998). The mutation rate μ was set at 10^{−3}. For that second set of simulations, we considered a population decline scenario (with *N*_{0} = 100, *N*_{1} =10,000, and *T*_{a} = 500), a population expansion scenario (with *N*_{0} = 10,000, *N*_{1} =100, and *T*_{a} = 500), and a stable population scenario, taking the (constant) population size as the harmonic mean of the population size change from 100 to 10,000 for *T*_{a} = 500 generations; *i.e.*, *N*_{0} = *N*_{1} = 464. This second set of simulations therefore consisted of seven sets of parameters: three mutation models were considered for the stable population scenario (the SMM and the two GSMs), and two mutation models were considered for each of the declining and expanding population scenarios (the two GSMs).

Microsatellite data were simulated with Simcoal2 (Laval and Excoffier 2004), which generates samples of genes under various demographic models, using a discrete-generation coalescent algorithm. Discrete-generation algorithms produce simultaneous and multiple coalescences, which are canceled out in the continuous-time approximation of the coalescent. There might therefore be a slight discrepancy between the coalescence rate in the discrete-generation algorithm and the continuous-time approximation that is assumed in Msvar, particularly for large sample sizes and small effective population size (see, *e.g.*, Figure S2 in Cornuet *et al.* 2008). However, we find it more relevant to simulate the data without relying on approximations. Each data set consisted of a sample of 50 diploid individuals, genotyped at 10 unlinked microsatellite loci. This sampling scheme is consistent with empirical studies that inferred past demographic changes using Msvar: from an exhaustive survey of the literature (Table S1), we found that the median numbers of microsatellite loci and sampled individuals across data sets were 11.5 and 30, respectively.

For each set of parameters, we simulated five microsatellite data sets to have replicates from the same underlying demographic and mutation model. We therefore obtained a total of 120 simulated data sets for the first set of simulations and 35 data sets for the second set. For each set, we calculated the mean and standard deviation over the five replicates of the expected heterozygosity *H*_{e} (Nei 1978), the observed number of alleles *N*_{a}, the range in allele size *A*_{r}, and the variance of allele range *V*_{a}, using Arlequin (Excoffier *et al.* 2005).

#### Parameterization of Msvar:

In Msvar, the posterior distribution of the model parameters is computed by means of a MCMC method using the Metropolis–Hastings algorithm (Metropolis *et al.* 1953; Hastings 1970). The likelihood is calculated from the genealogical history of the sample of genes, represented as a sequence of events (coalescences and mutations, see Beaumont 1999).

We used version 1.3 of Msvar, which provides separate estimates for *N*_{0,} *N*_{1,} μ, and *T*_{a} (Storz and Beaumont 2002). This implementation of Beaumont’s (1999) method, available at http://www.rubic.rdg.ac.uk/∼mab/stuff/, relies upon a hierarchical model where demographic and mutational parameters are allowed to vary among loci. The extent of interlocus variation is set by priors, as described in File S1. The parameter values reported in the following correspond to the mean of *N*_{0}, *N*_{1}, *T*_{a}, and μ, across loci. To test whether the method could retrieve information from the data, we chose relatively flat priors on the mean parameters, including the mutation rate.

#### Implementation:

Analyses were run on a Beowulf cluster made of 19 computer nodes, with CPUs ranging from biprocessors AMD Opteron monocore running at 1.8 GHz to biprocessors Intel Xeon quadcore running at 2.0 GHz. For each of the simulated data sets, three independent Msvar analyses were performed, with different starting values of the model parameters and different sets of seeds for the random number generator. For the first set of simulations (strict SMM with population size change), each Markov chain was initially run for 10^{9} steps and was thinned to 40,000 output lines by recording parameter values every 25,000 steps. In the absence of convergence, longer chains were run (see results). For the second set of simulations (strict SMM with stable population and GSM with stable, expanding, and declining populations), each Markov chain was run for 3 × 10^{9} steps, with parameter values recorded every 30,000 steps. The first 10% of steps of the chains were discarded as burn-in. For each data set, convergence was assessed by computing the multivariate extension of Gelman and Rubin’s diagnostic (Brooks and Gelman 1998) on the three independent Markov chains. Gelman and Rubin’s diagnostic is based on the computation of the ratio of the pooled-chains variance over the within-chain variance, which should be close to 1 if the chains converge to the target distribution. The multivariate Gelman and Rubin’s diagnostic was calculated from the means **M** ≡ {*M*_{N0}, *M*_{N1}, *M*_{Ta}, *M*_{μ}} and standard deviations **V** ≡ {*V*_{N0}, *V*_{N1}, *V*_{Ta}, *V*_{μ}} of Msvar parameters across loci, using the CODA package (Plummer *et al.* 2006) implemented in the statistical software R (R Development Core Team 2009). Although it might be recommended to run more chains to compute Gelman and Rubin’s diagnostic (*e.g.*, Nettel *et al.* 2009 used eight independent chains), the computational burden prevented us from running more than three chains in the present study.

#### Analysis of Msvar outputs:

Msvar outputs were analyzed by focusing on two issues: (i) the performance of Msvar at detecting past demographic changes and (ii) the precision of Msvar estimates of the model parameters. For each of the simulated data sets, we combined the three Markov chains before running the following analyses. The strength of evidence of population expansion *vs.* population decline (and vice versa) was evaluated using Bayes factors (Jeffreys 1961; Kass and Raftery 1995), as suggested by Beaumont (1999) and Storz and Beaumont (2002). The Bayes factor is a ratio where the numerator is the posterior probability of one model divided by its prior probability and the denominator is the posterior probability of an alternative model divided by its prior probability (Gelman *et al.* 1995). With identical priors for the population decline and the population expansion models (*i.e.*, identical priors for *N*_{0} and *N*_{1}), the Bayes factor for, *e.g.*, a population decline is the ratio of the posterior probability of a population decline divided by the posterior probability of a population expansion. This ratio can be estimated by counting the number of states in the chain in which the population has declined (*i.e.*, *N*_{0}/*N*_{1} < 1) and then dividing this by the number of states in which the population has expanded (*i.e.*, *N*_{0}/*N*_{1} > 1) (see Storz and Beaumont 2002).

We estimated the marginal posterior distributions of the model parameters using the LOCFIT package (Loader 1999) implemented in R (R Development Core Team 2009). Point estimates of natural parameters *N*_{0}, *N*_{1}, *T*_{a}, and μ were computed from the mode of their marginal posterior distribution. The 90% highest probability density (HPD) intervals were computed with the CODA package. We also estimated the marginal posterior distributions of the scaled parameters θ_{0} ≡ 4*N*_{0}μ, θ_{1} ≡ 4*N*_{1}μ, and *t*_{f} ≡ *T*_{a}/(2*N*_{0}), and we computed point estimates and 90% HPD intervals for these scaled parameters. For each demographic scenario considered, we calculated the absolute value of the bias for both natural and scaled parameters over the five replicated data sets.

#### Detection of population size change with Bottleneck and the *M*-ratio test:

Finally, for the first set of simulations (strict SMM with population size change), we compared the performance of Msvar to detect genetic signatures of demographic changes with the two most widely used moment-based methods available for microsatellite data. First, we analyzed the data sets using the method developed by Cornuet and Luikart (1996) and implemented in the software package Bottleneck v.1.2 (Cornuet and Luikart 1996). Wilcoxon signed-rank tests were performed to determine if a data set exhibited a significant number of loci with heterozygosity excess as expected in bottlenecked populations (Luikart *et al.* 1998) or with heterozygosity deficiency as expected in expanding populations (Cornuet and Luikart 1996). Second, we calculated Garza and Williamson’s (2001) *M* ratio on the 60 data sets corresponding to population declines. We compared empirical values of the *M* ratio to 95% critical values (*M*_{c}) derived from 10,000 simulations of stable populations using the program Critical_M. Simulations were performed using the true value of θ_{1} (θ_{1} = 4, 40, and 400 in the scenarios considered) and assuming a strict stepwise mutation model. We considered that an *M* ratio below the critical value *M*_{c} was indicative of a population decline.

## RESULTS

#### Genetic diversity of the simulated data sets, under a strict SMM:

For the first set of simulations (strict SMM with population size change), the expected heterozygosity *H*_{e}, the number of alleles *N*_{a}, and the range in allele size *A*_{r} are reported in Table 1. For contraction scenarios, *H*_{e} ranged from 0.24 to 0.94. *N*_{a} ranged from 2.3 to 23.7 and *A*_{r} varied from 1.3 to 39.2. In agreement with theoretical expectations, *H*_{e}, *N*_{a}, and *A*_{r} increased with *N*_{1}, the genetic diversity in the current population being sustained by large ancestral populations. Furthermore, genetic diversity decreased with increasing *T*_{a}, the loss of genetic diversity being more pronounced for long contraction events. For expansion scenarios, *H*_{e} ranged from 0.30 to 0.58. *N*_{a} ranged from 2.7 to 4.9 and *A*_{r} varied from 1.8 to 4.0. In agreement with theoretical expectations, *H*_{e}, *N*_{a}, and *A*_{r} increased with increasing *T*_{a} since the number of mutations that segregate in the population increases with the age of the expansion event. We also observed a tendency for genetic diversity to increase with increasing *N*_{0}, although this trend was not clear cut.

#### MCMC convergence:

In the following, we used Gelman *et al.*’s (2004) rule of thumb, which suggests that values of the multivariate Gelman and Rubin’s convergence diagnostic between 1.0 and 1.1 indicate reasonable convergence, whereas values >1.1 indicate poor convergence. Of the 120 analyses of the first set of simulations, 67 converged after 10^{9} steps (Table S2). The average computational time of these chains was 1.5 days for expansions and 3 days for contractions. The 53 nonconverged analyses were run again for 3 × 10^{9} steps and were thinned to 120,000 output lines by recording parameter values every 25,000 steps. Of these 53 analyses, 20 converged after 3 × 10^{9} steps, which took on average 20 days per chain. Finally, the last 33 nonconverged analyses were run for 1.5 × 10^{10} steps, which took 60 days per chain on average. Of these 33 analyses, 16 converged after 1.5 × 10^{10} steps. Therefore, a total of 17 analyses of 120 (14.2%) did not converge after 1.5 × 10^{10} steps. Most of these nonconverged analyses corresponded to recently and severely bottlenecked populations (*T*_{a} < 500 and *N*_{0}/*N*_{1} = 0.001; Table S2). However, visual inspection of the three chains in the nonconverged analyses, as well as the similarity of the marginal posterior distributions, suggested that the chains were close to equilibrium. Therefore, we included the 17 nonconverged analyses in our results. The cumulative computation time for the completion of all the analyses included in our study exceeded 276 × 10^{3} hr (33.5 years). There was a significantly positive correlation between the time of convergence and the average range in allele size *A*_{r} in the sample for both contractions (Spearman’s ρ = 0.82; *P* < 0.001) and expansions (Spearman’s ρ = 0.49; *P* < 0.001).

#### Detection of demographic events with Msvar, under a strict SMM:

Bayes factors (BF) were computed for each of the 120 analyses of the first set of simulations and interpreted following Jeffreys (1961): BF ≥ 10 indicate strong support, BF ranging from 3 to 10 indicate substantial support, BF ranging from 0.33 to 3 indicate no support, and values <0.33 indicate false detection of contraction or expansion. In 85 analyses of 120 (70.8%), Bayes factors indicated a change in population size consistent with the simulated scenario with substantial to strong support (BF ≥ 3 and BF ≥ 10, respectively; see Figure 1). Of the 60 Markov chains corresponding to contraction scenarios, 41 (68.3%) indicated a population decline (BF ≥ 3), of which 40 (97.6%) showed strong support (BF ≥ 10). Fifteen of these 40 analyses (37.5%) did not converge. Of the 60 analyses corresponding to expansion scenarios, 44 (73.3%) indicated a population expansion (BF ≥ 3), of which 34 (77.3%) showed strong support (BF ≥ 10). Two of these 34 analyses (5.9%) did not converge. Overall, all the ancient (*T*_{a} ≥ 50) and severe demographic changes (*N*_{0}/*N*_{1} ≤ 0.01 for contractions and *N*_{0}/*N*_{1} ≥ 100 for expansions) were detected with substantial to strong support (Figure 1). By contrast, recent declines and expansions (*T*_{a} = 10) were largely undetected (BF < 3), except for strong contractions (*N*_{0}/*N*_{1} = 0.001). Moreover, weak contractions (*N*_{0}/*N*_{1} = 0.1) were largely undetected whatever their time of occurrence and one false expansion was even detected for an ancient and weak bottleneck (BF < 0.33, *T*_{a} = 500, *N*_{0}/*N*_{1} = 0.1).

#### Comparison of Msvar with moment-based methods:

Because Bayes factors cannot be formally compared to *P*-values, we were not able to use the same criterion for detecting population size change with Msvar, Bottleneck, and the *M*-ratio test. Therefore, we reported in Figure 2 the criteria that are generally used in empirical studies: BF ≥ 3 for Msvar, the result of the Wilcoxon signed-rank tests at the α = 0.05 level for Bottleneck (Cornuet and Luikart 1996), and an *M* ratio below the critical value *M*_{c} (Garza and Williamson 2001) for the *M*-ratio test.

Bottleneck detected a significant excess of heterozygosity in only 5 of the 60 data sets corresponding to contraction scenarios (8.3%). Ancient events (*T*_{a} = 500) were never detected whatever their severity. Moreover, there was no clear relationship between the rate of detection of population decline and the severity of the event. Finally, 4 data sets corresponding to ancient contractions (*T*_{a} ≥ 100) showed significant heterozygote deficiency, hence supporting population expansions. Contrastingly, Bottleneck detected a significant deficiency in heterozygosity in 35 of 60 data sets corresponding to expansion scenarios (58.3%). Garza and Williamson’s (2001) *M*-ratio method correctly detected a signal of contraction in 32 of 60 data sets (53.3%). The rate of detection was higher for ancient (*T*_{a} ≥ 50) and moderate-to-severe population declines (*N*_{0}/*N*_{1} ≥ 0.01), with 26 significant tests of 30. Recent (*T*_{a} = 10) and/or weak declines (*N*_{0}/*N*_{1} = 0.1) were barely detected (6 significant tests of 30).

#### Estimation of demographic and mutational parameters with Msvar, under a strict SMM:

Demographic and mutational parameters were estimated for all data sets, by combining the three Markov chains run for each data set. We assessed the quality of the estimates by examining the marginal posterior distributions of the parameters, compared to their prior distributions. We summarized these results by calculating the modes and the 90% HPD intervals for each data set (Figure S1, Figure S2, Figure S3, and Figure S4) as well as the absolute value of the bias and the average HPD range over the five replicate data sets for each of the 24 demographic scenarios (Figures 3 and 4).

##### Estimates of the natural parameters N_{0}, N_{1}, T_{a}, and μ:

Overall, the marginal posterior distributions of the demographic parameters *N*_{0}, *N*_{1}, and *T*_{a} were wide and departed only slightly from the priors (*e.g.*, Figure 5, A and B). The estimated 90% HPD limits were therefore broad (Figure 3), ranging from −4 to 8 in log_{10} scale (Figure S1 and Figure S2).

For contractions, replicated data sets tended to provide more consistent results for old and severe events, compared to recent events (*T*_{a} = 10) or events of low severity (*N*_{0}/*N*_{1} = 0.1) (Figure S1). The precision of the demographic parameter estimates tended to increase with increasing severity of the demographic change (measured by the ratio *N*_{0}/*N*_{1}) and the time of the event: (i) the 90% HPD range of the demographic parameter estimates decreased with increasing *N*_{0}/*N*_{1} (Figure 3 and Figure S1); (ii) for moderate to strong contractions (*N*_{0}/*N*_{1} > 0.1), the 90% HPD range decreased with increasing *T*_{a}; and (iii) for *N*_{0}/*N*_{1} = 0.1, the 90% HPD range was the lowest for intermediate values of *T*_{a}. The absolute value of the bias of *N*_{0} estimates tended to be lower than that of *N*_{1} and was maximized for recent events (*T*_{a} = 10).

The quality of the estimates of *N*_{0}, *N*_{1}, and *T*_{a} was poorer for expansions, compared to contractions. The marginal posterior distributions were not sharply peaked and did not depart markedly from the priors. The 90% HPD limits were wide and the absolute value of the bias was high, overall (Figure 3). This was true whatever the severity of the event and its time of occurrence. It is noteworthy that, with few exceptions, all demographic parameter estimates differed markedly across replicate data sets (Figure S2). We noted that for a given expansion severity, estimates of *N*_{0} increased with *T*_{a}, while estimates of *N*_{1} decreased with *T*_{a} (Figure S2).

For contractions and expansions, the marginal posterior distributions of μ departed only slightly from the prior distributions, whose mean was set at α_{μ} = −4 on a log_{10} scale. Because the true mutation rate μ of the simulated data sets was set at −3 on a log_{10} scale, the mutational parameter μ was therefore systematically underestimated, as already pointed out by Milton *et al.* (2009). The 90% HPD intervals of the marginal posterior distributions of μ were wide (data not shown).

Finally, we examined the patterns of correlation between natural parameters to assess the performance of Msvar to estimate natural parameters separately. We observed strong correlations between natural parameters of the model. Overall, both *N*_{0} and *N*_{1} were negatively correlated with the mutational parameter μ and there was a positive correlation between *N*_{0} and *T*_{a} (Figure S5). The correlations were stronger for more severe events and more ancient events. Furthermore, the correlations were more pronounced for contractions than for expansions.

##### Estimates of the scaled parameters θ_{0}, θ_{1}, and t_{f}:

Scaled parameters were overall much more precisely estimated than the natural parameters for contractions, whereas they were poorly estimated for expansions. As with the natural parameters, the quality of the estimates depended upon the severity of the demographic change and its time of occurrence.

For contractions, the marginal posterior distributions of the scaled parameters θ_{0}, θ_{1}, and *t*_{f} were very peaked and departed markedly from the prior distributions (*e.g.*, Figure 5, C and D), except for contractions of low severity (*N*_{0}/*N*_{1} = 0.1). The precision (low bias, narrow 90% HPD interval) increased with increasing severity of the event and time of occurrence (Figure 4 and Figure S3). In particular, estimates of θ_{1} and *t*_{f} were overall very precise for moderate to severe bottlenecks (*N*_{0}/*N*_{1} < 0.1), except for very recent events (*T*_{a} = 10). Although θ_{0} was also well estimated for ancient declines (*T*_{a} > 50) from moderate to strong severity, the bias and the range of 90% HPD intervals were larger compared to those of θ_{1}. Replicate data sets provided consistent results for θ_{1} and *t*_{f}, for moderate and strong contractions (*N*_{0}/*N*_{1} < 0.1) that occurred >10 generations ago (*T*_{a} > 10). Larger variation across replicate data sets was observed for θ_{0}.

For expansions, the marginal posterior distributions of the scaled parameters θ_{0} and θ_{1} were peaked and departed markedly from the prior distributions (data not shown). However, the mode of the marginal posterior distributions for θ_{0} departed markedly from the true simulated value, resulting in severe biases (Figure 4 and Figure S4). By contrast, the scaled parameter θ_{1} exhibited low bias in all scenarios (Figure 4), although the 90% HPD intervals were wide, especially for weak and recent expansions (*N*_{0}/*N*_{1} = 10; *T*_{a} > 10). We noted that the marginal posterior distributions of θ_{0} and θ_{1} were skewed, respectively, to upper and lower values (Figure 4 and Figure S4). The 90% HPD intervals of the marginal posterior distributions for the time parameter *t*_{f} were wide in most conditions (Figure 4 and Figure S4) and estimates were severely biased, except for *T*_{a} = 500. Both the 90% HPD interval and the bias decreased with increasing *T*_{a}. We observed large variations across replicate data sets for all scaled parameters in almost all situations, particularly for *t*_{f}.

#### Influence of the mutation model in Msvar:

Of the 30 analyses presented in Figure 6 for the GSM, 12 (40%) did not converge after 3 × 10^{9} steps. Out of these, 9 analyses (75%) concerned the data sets generated with the strongest GSM (*p* = 0.74). With data generated under the moderate GSM (*p* = 0.22), Msvar successfully detected a population decline for the five simulated data sets and a population expansion for four simulated data sets of five. However, Msvar detected a false signal of population decline for two data sets of five that were simulated under a stable population scenario (Figure 6). Under a strong GSM (*p* = 0.74), Msvar detected a signal of population decline with strong support (BF ≥ 10), whatever the simulated scenario (Figure 6).

The quality of Msvar estimates of scaled parameters for the moderate GSM (*p* = 0.22) was very similar to that observed for the strict SMM, with very precise estimates of θ_{1} and *t*_{f}, a slightly larger bias, and 90% HPD intervals for θ_{0} compared to θ_{1} in contraction scenarios and poorer estimates, with large variations across replicate data sets, in expansion scenarios (Figure S6). For stable population scenarios, both the strict SMM and the moderate GSM (*p* = 0.22) produced unbiased estimates of θ_{0} and θ_{1}, but with very large 90% HPD intervals. Note that in the absence of population size change, estimates of *t*_{f} are meaningless. Very consistently, Msvar produced biased estimates of the model parameters, with very narrow 90% HPD intervals, for all the data sets generated under the strong GSM (*p* = 0.74) (Figure S6).

## DISCUSSION

#### Comparing Msvar, Bottleneck, and the *M*-ratio test:

Bottleneck performed poorly in detecting population declines from our simulated data sets under a SMM, with only 5 significant tests of 60. The statistical power of Bottleneck for population declines is much lower when microsatellite loci evolve under a strict SMM than under an infinite-allele model (Cornuet and Luikart 1996) or a GSM (Leblois *et al.* 2006). This may partly explain the low performance of Bottleneck in our comparative study. Our results for weak population declines (*N*_{1}/*N*_{0} = 10) are in agreement with previous simulation-based evaluations, given the set of demographic and mutational parameters considered here (see, *e.g.*, Figure 3B in Cornuet and Luikart 1996). For moderate to severe population declines (*N*_{1}/*N*_{0} ≥ 100), however, the rate of detection was lower in our study than in Cornuet and Luikart (1996). Two possible reasons may explain this discrepancy. First, the average heterozygosity in our simulated data sets was overall higher than in Cornuet and Luikart (1996) who considered a variable mutation rate across loci and simulations, to cover a range of heterozygosities per set of parameters. Second, we simulated an exponential change of population size, whereas Cornuet and Luikart (1996) assumed an instantaneous reduction of population size in their simulation-based tests. The impact of the shape of the demographic change on the performance of Bottleneck has not been studied yet. Consistent with our study, the simulation-based evaluation of Bottleneck by Leblois *et al.* (2006) also showed a low statistical power of the method. Interestingly, Bottleneck performed largely better for expansions (58.3%) than for contractions (8.3%), given the model parameter values of our study.

We found that the *M*-ratio test was more efficient than Bottleneck, which is consistent with Leblois *et al.* (2006), for retrieving signals of population declines from our simulated data sets (32 significant tests of 60). The rate of detection was higher for ancient and moderate-to-severe declines, while recent and weak declines were barely detected. These results are consistent with previous simulation-based studies that have shown that the *M*-ratio test has low statistical power for small θ_{1} values (here, θ_{1} = 4, see Garza and Williamson 2001; Williamson-Natesan 2005) and for recent population declines (see Williamson-Natesan 2005; Leblois *et al.* 2006). Here, we applied the *M*-ratio test by comparing the statistic *M* estimated from the data with the expected distribution of that statistic, conditionally on the true value of θ_{1}. This procedure increased the statistical power of the *M*-ratio test. Since the true value of the parameter θ_{1} is generally unknown in real situations, Garza and Williamson (2001) recommended the use of *M*_{c} = 0.68 as a conservative threshold for the critical value. The reanalysis of our data sets with *M*_{c} = 0.68 (*e.g.*, as in Leblois *et al.* 2006) resulted in a lower rate of detection (22 significant tests of 60), but in similar qualitative trends (higher rate of detection for severe and ancient population declines).

Given the set of demographic and mutational parameters used in our study, and using the decision criteria recommended by the developers of each method, Msvar clearly outperformed the *M*-ratio test and Bottleneck for detecting population size change. While Msvar correctly detected 68.3% of the declines, the *M*-ratio test and Bottleneck detected only 53.3% and 8.3%, respectively, of the declines. Any population decline detected by the *M*-ratio test and Bottleneck was also recovered by Msvar, apart from one case of weak recent decline that was identified only by Bottleneck (*T*_{a} = 10 and *N*_{0}/*N*_{1} = 0.1). Therefore, our study does not support the previous claims that the *M*-ratio test and Bottleneck are best suited to detect recent population declines, whereas Msvar is more appropriate to detect ancient contractions (Garza and Williamson 2001; Williamson-Natesan 2005). Moreover, while Msvar detected 73.3% of the population expansions, Bottleneck detected only 58.3% of the expansions. Any expansion detected by Bottleneck was also recovered by Msvar.

#### Performance of Msvar: What does coalescent theory tell us?

Not surprisingly, we found that the performance of Msvar to infer past demography strongly depended on the information available in the data, which may be inferred from coalescent theory. Coalescent theory indeed predicts that variations in population size strongly affect the shape of gene genealogies, which are star shaped with long terminal branches in expanding populations and shallower in declining populations (Figure 7 and Hein *et al.* 2005).

Coalescent theory further shows that only scaled parameters can be directly estimated from the data (Tavaré *et al.* 1997; Nordborg 2007). Indeed, all parameters in coalescent models are scaled, and the likelihood function in Msvar makes no exception (Beaumont 1999). Hence, inference of unscaled quantities such as population size, or time measured in generations, requires external information. In our study, unscaled parameters were therefore much less precisely estimated than the scaled ones (Figure 5) and were also highly correlated (Figure S5; see also Figure 5 in Storz *et al.* 2002). We deliberately chose poorly informative priors, to test the capacity of Msvar to retrieve information from the data only. In empirical studies, more informative priors of the natural parameters are usually specified. We acknowledge that Msvar offers a principled approach for providing prior information on the mutation rate, to recover posterior densities for natural parameters. Yet it should be borne in mind that precise estimates of unscaled parameters may then largely stem from the specification of the priors. Imagine that analyses were performed using a prior distribution for the mutation rate with very low standard deviation (*i.e*., σ_{μ} close to zero). We would then necessarily recover the same level of precision for the natural parameters and the scaled parameters. Yet this improved precision may come at the expense of accuracy, if the prior distribution for the mutation rate departs from its true distribution.

##### Scenarios of population decline:

Msvar was very efficient for detecting population declines. However, its performance for detecting change in population size and accurately estimating the model parameters was lowest for recent events (*T*_{a} = 10) of low-to-moderate severity (*N*_{0}/*N*_{1} ≥ 0.01), as well as for events of low severity (*N*_{0}/*N*_{1} = 0.1). This is expected since, for very recent declines, the gene genealogy can barely be distinguished from that expected in a stable population with population size *N*_{1} (Figure 7A). Interestingly, for *N*_{0}/*N*_{1} = 0.1, the performance of Msvar was maximized for intermediate values of *T*_{a}, particularly with respect to the precision of θ_{1} estimates. This might be easily understood by considering that for ancient events (*T*_{a} = 500) most coalescence events occur while *N*(*t*) is close to the current population size (see, *e.g.*, Figure 7B). This might be further quantified by calculating the expected number *j* of ancestral lineages at (scaled) time τ, from which a sample of *n* genes is descending. This number has a known distribution in a constant-size population (Tavaré 1984), and Leblois and Slatkin (2007) extended Tavaré’s (1984) formula in the case of an exponentially growing or declining population. Using their Equation 2 that gives an expression for the probability Pr(*j* | *n*) that a sample of *n* genes has *j* ancestors *T*_{a} generations ago, we may compute the expected number *m* of lineages at the time of the population size change as*a*_{(i)} ≡ *a*(*a* + 1) … (*a* + *i* − 1), *a*_{[i]} ≡ *a*(*a* − 1) … (*a* − *i* + 1), and *T*_{a} = 500 and *N*_{0}/*N*_{1} = 0.1, we get *m* = 1.43, which confirms that most coalescence events are expected to occur in the current population with this set of parameter values.

For moderate to severe contractions (*N*_{0}/*N*_{1} ≤ 0.1), both the bias and the 90% HPD range of θ_{0} decreased with increasing *T*_{a}. Using Equation 1, we found that the expected number *m* of lineages at the time of the event varies between 48.49 and 2.20 for *T*_{a} varying from 10 to 500 and for *N*_{0}/*N*_{1} = 0.01. This indicates that more coalescence events are expected to occur in the declining population when the event is older (see also Figure 7C). In contrast, θ_{1} was overall precisely estimated (see Figure 4 and Figure S3). This is so because, for the scenarios considered here, a large part of the genealogy depends upon the ancestral history, with several lineages coalescing in the ancestral population (see, *e.g.*, Figure 7, B and C) at a rate that depends upon θ_{1}. Had we considered older events (*T*_{a} > 500), though, thereby decreasing the number of lineages in the ancestral population, it is likely that the precision of θ_{1} estimates would have declined.

In summary, most scenarios of population decline result in gene genealogies with large times to the most recent common ancestor (TMRCAs). With the set of model parameters considered here, since a large part of gene genealogies depends upon θ_{1}, this latter parameter is generally precisely and accurately estimated. Contrastingly, θ_{0} can be precisely and accurately estimated only if the demographic event is severe and ancient. If the change in population size is too recent, provided that it is not too pronounced, θ_{0} estimates tend to converge to the true value of θ_{1}, and no change of population size is detected. If the difference in population size is weak, then the difference in coalescence rates before and after the event is not sufficient for Msvar to detect a population size change and to provide precise estimates of θ_{0} and θ_{1}.

##### Scenarios of population expansion:

Msvar was also very efficient for detecting expansions. Nevertheless, the estimates of the scaled current population size θ_{0} were more severely biased and less precise, compared to scenarios of population decline, for the same relative severity of the event. This may be explained by the fact that expanding populations result in young genealogies with short TMRCAs (compare Figure 7, A–C, to 7, D–F) and hence rare mutation events. We found that the absolute value of the bias increased with *N*_{0}. We further found that both the 90% HPD range and the absolute value of the bias of θ_{0} decreased with increasing *T*_{a} (Figure 4 and Figure S4). Using Equation 1, we found that the expected number *m* of lineages at the time of the event varies between 48.49 and 2.20 for *T*_{a} varying from 10 to 500 and for *N*_{0}/*N*_{1} = 100. This indicates that the number of coalescence events in the expanding population is expected to increase as *T*_{a} increases (see also Figure 7F). More generally, the likelihood surface for expanding populations is complex (Beaumont 1999). In particular, as the genealogies become more star shaped, the joint posterior distribution of θ_{0} and *t*_{f} reduces to a ridge along a line log_{10}(2μ*T*_{a}) = *k* independent of θ_{0} (Figure 4b in Beaumont 1999), which suggests that Msvar provides information on 2μ*T*_{a}, rather than on θ_{0} in expanding populations. It is worth noting that Wakeley *et al.* (2001) found similar results in inferring demographic history from single-nucleotide polymorphims (see a comparison in Beaumont 2004).

Estimates of θ_{1} had a low bias but a large 90% HPD range (Figure 4). Although the marginal posterior distributions of θ_{1} were generally peaked, they were flat tailed on the left. This is so, because large ancestral population sizes are not compatible with the low polymorphism observed in the data. Instead, a large range of small values of θ_{1} may be equally likely, provided the genealogy is star shaped.

#### Influence of the underlying demographic model:

Athough Msvar equally detected population declines and expansions, inferences of the demographic parameters were in general more accurate for declines than for expansions. In addition to the above argument from coalescent theory, the exponential model assumed for population size change may partly explain this pattern. For declines, the size of the declining population *N*(*t*) decreases sharply at *T*_{a} and converges rapidly to *N*_{0} (Figure 7, A–C). Therefore, most coalescence events occurring in the declining population take place while *N*(*t*) is close to *N*_{0}. For expansions, instead, the size of the expanding population *N*(*t*) increases smoothly at *T*_{a} before it converges rapidly to *N*_{0} (Figure 7, D–F). Therefore, a large proportion of the coalescence events that occur in the growing population take place while *N*(*t*) is close to *N*_{1} (compare, *e.g.*, Figure 7C to 7F). This can be expressed more formally by considering the harmonic mean of population sizes, which provides the coalescent rate during the change in population size (Hein *et al.* 2005). For, say, *T*_{a} = 500, the harmonic mean of an exponentially declining population with *N*_{0} = 100 and *N*_{1} = 10,000 is 464, which is strictly equal to the harmonic mean of an exponentially growing population with *N*_{0} = 10,000 and *N*_{1} = 100. Hence, the harmonic mean of a declining population is closer to its current size (*N*_{0}) than its ancestral size (*N*_{1}), while the reverse is true for expanding populations. Therefore, given the exponential model of population growth, one might expect poor statistical properties of θ_{0} estimates in expanding populations, compared to declining populations.

#### Robustness of Msvar to the misspecification of the mutation model:

Most importantly, our results suggest that Msvar is robust to moderate departures from a strict SMM, *e.g.*, a GSM with *p* ≤ 0.22, typical of those observed in the literature (see Figure 6 and Figure S6). However, severe departures from a strict SMM (here, a GSM with *p* = 0.74) led Msvar to detect a signal of population decline with strong support (BF ≥ 10), even in expanding populations (Figure 6 and Figure S6). This is not surprising since it has been recognized that violation of the assumptions of the SMM might induce severe bias in the inference of demographic history (Gonser *et al.* 2000). Indeed, mutations that arise under a strong GSM involve large changes in allele length, which produce some gaps in the distribution of allele types. The large resulting variance of allele range *V*_{a} is reminiscent of that observed with population decline (Storz and Beaumont 2002), even in expanding populations (compare Table 1 to Table S3).

#### Insights from empirical studies:

The better performance of Msvar compared to the *M*-ratio test and Bottleneck also emerged from the empirical studies that inferred past demographic changes from microsatellite data using Msvar and at least one of the *M*-ratio or Bottleneck methods (Table S1). We found indeed that Msvar detected a population decline whenever one of the moment-based methods provided a significant test. By contrast, a large number of population declines that were not detected with any of the moment-based methods were detected with Msvar. Unfortunately, the scarcity of expansion events detected in the literature (Table S1) prevented any empirical comparison of Msvar and Bottleneck for growing populations. Importantly, the average genetic diversity measured from our simulated data sets was not substantially different from that observed in empirical studies (compare Table 1 to Table S1).

Because of the large heterogeneity of the published results, we did not attempt to analyze the quality of Msvar estimates in empirical studies. Some studies used Msvar 0.4 (Beaumont 1999), hence providing estimates for scaled parameters, and some studies used Msvar 1.3 (Storz and Beaumont 2002), hence providing estimates for unscaled parameters. Only a handful of studies used both methods, and few provided estimates of the scaled parameters using Msvar 1.3, as in the present study. Finally, credibility intervals were often not reported or calculated using different methods, which hampered any comparison among studies.

#### Recommendation guidelines and conclusions:

Our simulation tests as well as an exhaustive survey of the literature clearly demonstrate that Msvar outperforms both the *M*-ratio test and Bottleneck for detecting population declines. Our study further shows that Msvar is also very efficient to detect population expansions and outperforms Bottleneck in that respect. However, to our knowledge, Msvar has only scarcely been applied on presumably expanding populations (see, *e.g.*, Hufbauer *et al.* 2004; Bonhomme *et al.* 2008; Wirth *et al.* 2008). Hence, we confidently recommend the use of Msvar for detecting past population size variation, even if this method is computationally demanding.

Most importantly, in contrast to the *M*-ratio test and Bottleneck, Msvar provides estimates of the parameters that characterize the population demographic history and the mutational model. Using Msvar 1.3 (Storz and Beaumont 2002), we have shown that the scaled parameters are more precisely estimated than the natural parameters. Although the latter are easier to interpret, our results clearly advocate drawing conclusions from inferences of θ_{0}, θ_{1}, and *t*_{f}. These parameters were precisely estimated for population declines, provided that the change in population size was neither too recent nor too weak, given the scenarios considered in our study. For expansions instead, both unscaled and scaled parameters were poorly estimated, although the method was efficient for detecting increase in population size. Hence, our results suggest that Msvar estimates in presumably expanding populations should be taken cautiously. We did not compare the performance of Msvar 0.4 (Beaumont 1999), which provides estimates for scaled parameters, to that of Msvar 1.3 (Storz and Beaumont 2002), which provides estimates for unscaled parameters, since the two versions differ by a number of other aspects. Msvar 0.4 assumes a basic (nonhierarchical) model, where the parameters are not allowed to vary among loci. The parameterization of Msvar 1.3 by means of a hierarchical model allows for some variation of the parameters among loci, which may provide a means of identifying aberrant loci. However, with broad priors on interlocus variation of the model parameters (**V**), Msvar 1.3 does not fully “borrow strength” from the different loci, *e.g.*, by simply pooling information (multiplying the likelihoods) across loci (Beaumont and Rannala 2004). This generally results in broader posterior distributions in the hierarchical Msvar 1.3 model, compared to the basic Msvar 0.4 model. Hence, although our results suggest that the use of scaled parameters should be preferred, further analyses are required to compare the performances of Msvar 0.4 and 1.3.

Finally, we recommend that inferences about population demographic change with Msvar should be interpreted cautiously in light of potential departures from the model assumptions. First, Msvar assumes that microsatellites evolve according to a strict SMM. Although a moderate departure from this mutation model, as classically measured with observations of spontaneous mutations (Ellegren 2000, 2004), does not seem to undermine Msvar performance (see Figure 6 and Figure S6), loci that evolve under a strong GSM may invalidate the approach by detecting false signals of population decline whatever the true demographic history. However, the hierarchical model implemented in Msvar 1.3 allows for variations in mutation and demographic parameters across loci and may thus limit the potential biases due to misspecifications of the mutation model. Indeed, Storz *et al.* (2002) argue that loci that strongly depart from the strict SMM shall be given less weight, thereby minimizing their impact on the inference made. Second, Msvar assumes that populations are isolated. Real populations, however, are in general connected by gene flow. It is now acknowledged that population structure and/or isolation by distance may result in incorrect inference of population demographic history (Pope *et al.* 2000; Leblois *et al.* 2006; Nielsen and Beaumont 2009; Chikhi *et al.* 2010; Peter *et al.* 2010). Finally, further work is needed to evaluate how Msvar performs when the demography is more complex, *e.g.*, with successions of population declines and expansions.

## Acknowledgments

We are grateful to L. Chikhi for sharing an unpublished manuscript, as well as to Mark A. Beaumont and two anonymous reviewers for their constructive criticism of this article. This work was supported by the Plan Pluri-Formation “Evolution et Structure des Ecosystèmes” (2004–2008) from the Muséum National d’Histoire Naturelle (MNHN), by a Nouragues Research grant from the Centre National de la Recherche Scientifique (CNRS) programme Amazonie I (2007), by the CNRS programme Amazonie II (2008–2011), and by the Agence Nationale de la Recherche programme blanc Études de Méthodes Inférentielles et Logiciels pour l'Évolution (EMILE) NT09-611697. This work is part of Christophe Girod’s Ph.D., who was supported by a grant from the French Ministry of Research (2007–2010). Part of this work was carried out by using the resources of the Computational Biology Service Unit from the MNHN (CNRS Unité Mixte de Service 2700).

## Footnotes

Supporting information is available online at http://www.genetics.org/cgi/content/full/genetics.110.121764/DC1.

- Received August 4, 2010.
- Accepted February 21, 2011.

- Copyright © 2011 by the Genetics Society of America