Blackwell Publishing LtdColonization history of the Swiss Rhine basin by the bullhead (Cottus gobio): inference under a Bayesian spatially explicit framework
SAMUEL NEUENSCHWANDER,*† CARLO R. LARGIADÈR,*‡ NICOLAS RAY,* MATHIAS CURRAT,*§PASCAL VONLANTHEN*¶ and LAURENT EXCOFFIER**CMPG, Zoological Institute, University of Bern, Baltzerstrasse 6, 3012 Bern, Switzerland, †Department of Ecology and Evolution, Biophore, UNIL-Sorge, 1015 Lausanne, Switzerland, ‡Institute of Clinical Chemistry, University Hospital, University of Bern, Inselspital, 3010 Bern, Switzerland, §Laboratory of Anthropology, Genetics and Peopling history (AGP), Department of Anthropology and Ecology, University of Geneva, 1227 Geneva, Switzerland, ¶Aquatic Ecology, Zoological Institute, University of Bern, Baltzerstrasse 6, 3012 Bern, Switzerland
Abstract
The present distribution of freshwater fish in the Alpine region has been strongly affectedby colonization events occurring after the last glacial maximum (LGM), some 20 000 yearsago. We use here a spatially explicit simulation framework to model and better understandtheir colonization dynamics in the Swiss Rhine basin. This approach is applied to theEuropean bullhead (Cottus gobio), which is an ideal model organism to study fish pastdemographic processes since it has not been managed by humans. The molecular diversityof eight sampled populations is simulated and compared to observed data at six microsatelliteloci under an approximate Bayesian computation framework to estimate the parameters ofthe colonization process. Our demographic estimates fit well with current knowledgeabout the biology of this species, but they suggest that the Swiss Rhine basin was colonizedvery recently, after the Younger Dryas some 6600 years ago. We discuss the implication ofthis result, as well as the strengths and limits of the spatially explicit approach coupled tothe approximate Bayesian computation framework.
Received 27 May 2007; revision received 6 October 2007; accepted 26 October 2007
Introduction
The spatial distribution of populations is deeply influencedby climatic fluctuations, such as periods of range con-tractions during glaciations as well as range shifts andexpansions during more favourable interglacial periods.Such spatial and demographic expansions greatly affectthe genetic variation and structure of populations andspecies (Avise 1998; Hewitt 2000; Ray et al. 2003; Excoffier2004). Previous studies have shown theoretically (Ibrahimet al. 1996; Austerlitz et al. 1997) and empirically (e.g. forsalmonid fish species, Bernatchez & Wilson 1998) thatpostglacial expansions lead to a global loss of geneticvariation. Most phylogeographical studies on European
freshwater fishes have been carried out on economicallyimportant salmonid species such as brown trout (Bernatchez2001), whitefish (Hansen et al. 1999), Atlantic salmon (Potvin& Bernatchez 2001), or grayling (Weiss et al. 2002). However,these species have been subject to massive artificial trans-plantations for more than a century (Largiadèr & Scholl1996; Largiadèr et al. 1996; Brunner et al. 2001; Douglas &Brunner 2002). It is therefore difficult to distinguish, at thegenetic level, between the effects of a natural colonizationand an artificial introduction. To circumvent this problem,several studies have used the European bullhead (Cottusgobio) as a model organism for investigating the colonizationhistory of European freshwater fish species (Hänfling &Brandl 1998; Englbrecht et al. 2000; Kontula & Vainola 2001;Hänfling et al. 2002; Volckaert et al. 2002; Slechtova et al.2004). Compared to other species, C. gobio has been of littleeconomic importance, and it was therefore rarely managed
by humans (Mann 1971), which suggests that its presentdistribution is mainly the result of a natural colonizationprocess (Hänfling & Brandl 1998; Englbrecht et al. 2000).Cottus gobio is thus a very valuable organism to better under-stand past demographic processes that occurred in theAlpine river systems since the last glacial maximum (LGM).
During the last Würm glaciation, large parts of the Alpineregion were covered by ice (Jäckli 1962; Hantke 1980).Freshwater fish species were therefore either completelyabsent from this area or their presence was very restricted.The retreat of the glaciers following LGM, some 20 000 yearsago (Jäckli 1962; Hantke 1980) allowed these species torecolonize the Alps. It is therefore likely that the currentdistribution of freshwater fishes in the Swiss Rhine riversystem must largely result from these early colonizationevents. Recently, Vonlanthen et al. (2007) reported highlevels of differentiation between populations of C. gobio atthe nuclear level in the Swiss Rhine river system, a strongpopulation genetic structure that has also been reported inother parts of Europe (Englbrecht et al. 2000; Hänfling et al.2002; Volckaert et al. 2002).
The aim of the present study is to better define the colo-nization history of C. gobio in the Swiss Rhine basin.Conventional population genetic inferential methods arebased on very simple demographic processes, which mayoften lack important sources of realism. New approaches,such as the reconstruction of genetic variation usingsimulations and the subsequent estimation of appropriateparameters, lead to more detailed conclusions about pastdemographic processes (Wilson & Balding 1998; Beaumont1999; Estoup & Clegg 2003; Estoup et al. 2004; Hamiltonet al. 2005). Because likelihood-based methods are sometimesextremely challenging or even impossible to compute for
complex models (Wilson & Balding 1998; Beaumont 1999),a new Bayesian approach, the approximate Bayesiancomputation framework (ABC, Beaumont et al. 2002;Marjoram et al. 2003) has been developed. This new methodavoids the calculation of likelihoods and thus allows one tostudy complex models. The potential of the ABC methodfor analysing complex demographic scenarios has beendemonstrated recently in several papers (Estoup & Clegg2003; Estoup et al. 2004; Hamilton et al. 2005). We thereforeuse this framework to estimate important demographicparameters of C. gobio colonization process, such as thetiming of the range expansion into the Swiss Rhine basin,migration rates and local effective sizes. Finally, we alsoassess the performance of our estimation procedure usingtest data sets, and discuss some limitations of the currentapproach.
Materials and methods
Genetic data
We based our estimation on a data set analysed in acompanion study (Vonlanthen et al. 2007), investigatingwhether the watershed between the Rhine and Rhône riversystem had been crossed by Cottus gobio. Vonlanthen et al.(2007) collected and genotyped several C. gobio populationsfrom the Rhine and Rhône basin at eight short tandemrepeat (STR) loci and they analysed the data using standardstatistical methods. In the present study, we use a subset ofthe genotyped populations to estimate parameters of thecolonization history of C. gobio in the Rhine river system.Our data set consists of eight populations of the Rhine riversystem (Fig. 1), for which 19–20 individuals per population
Fig. 1 Swiss Rhine river system withsample locations. The grey river segmentsrepresent the entire river system, whereasthe black segments are those used in thisstudy. The structured area illustrates themaximal extension of the glaciers duringthe LGM. The black dots represent thesampled locations (1, Simme; 2, Suze; 3,Birse; 4, Grenet; 5, Emme; 6, Saar; 7,Torneresse; 8, Orbe). For each location, 20fish were sampled except at location 1 withonly 19 sampled fish. The magnifying glassshows a schematic view of the vectorizedriver system subdivided into segments(demes). Each segment may be connectedto one (terminal segment), two, or several(river branching) neighbouring segments.For a given simulation, all segments havean identical length.
A B C E S T I M AT I O N O F C O L O N I Z AT I O N H I S TO RY 3
were genotyped. While all eight loci were found to bepolymorphic at the European level (Englbrecht et al. 1999;Vonlanthen et al. 2007), two are found to be almost orcompletely monomorphic in the Swiss Rhine data set. Sinceour estimation procedure relies on the simulation of geneticdata, we discarded these two loci from our analysis, as thislack of diversity at STR loci was not informative. Wetherefore restricted our analysis to the following six loci:Cgo1033PBBE Cgo56MEHU, Cgo18ZIM, Cgo33ZIM,Cgo34ZIM, and Cgo42ZIM (Englbrecht et al. 1999). Thesimulated data sets were compared to the observations onthe basis of the summary statistics listed in the secondcolumn of Table 4.
Simulation framework
The simulated data sets were generated using the programaquasplatche (Neuenschwander 2006a), which is derivedfrom the splatche program (Currat et al. 2004). Whilesplatche allows the simulation of genetic variation ofspecies living in a two-dimensional habitat (Ray et al. 2003;Currat & Excoffier 2004; Ray et al. 2005) aquasplatchewas specifically developed to simulate the genetic variationof species living in linear branching habitats, such as riversystems. The simulation framework is divided into twoparts: first the demography of a population is simulated ina forward process based on a vectorized river systemsubdivided into demes (subpopulations). Second, theinformation from the demographic simulation is used togenerate genetic samples from the population in a backwardcoalescent simulation process. For this study, we generated1 million simulations using parameter values drawn fromprior distributions.
Digital representation of the river system
In order to efficiently simulate the colonization history of theSwiss Rhine basin by C. gobio, we generated a vectorizedmap of the river system (Fig. 1). Contrary to raster maps,vector maps are ideal to represent linear features like riversystems, because object topology is more easily maintained,and it results in smaller file size. However, one difficultywe encountered in this study was to vectorize lakes that arenot linear features and which must therefore be transformed.We chose here to simply represent lakes by a stream locatedin the centre of the lake. The tributaries flowing into a lakewere therefore elongated up to this central stream. Thisrepresentation of the lakes seems adequate for C. gobio,since this species probably colonized lake tributaries byfollowing lake shores, because the deep waters of a lake areunsuitable for this species requiring running water. As abase map, we used a vectorized river system data setpublished by the Swiss Federal Institute of Topology(swisstopo 2001, scale 1:200 000). The river systems were
subdivided into segments. Each segment was characterizedby its upper and end node, as well as by several verticesdefining the form of the segment. Each segment can beeither connected to one (terminal segment), two, or more(river branching) neighbouring segments (Fig. 1). The riversystem of the resulting map was controlled for overall flowconnectivity, which means that each segment was part of ariver connected to the network. Using the gis softwarearcinfo (ESRI 2005), we exported information on thenodes and vertices of all segments into the format used byaquasplatche (Neuenschwander 2006a). We simplified theriver system to the area where we had sampled populations,and we removed parts of the river system where nopopulations had been sampled. We then subdivided thevectorized river system into segments of equal length.Because the size of a segment may have a great influenceon demographic parameters such as the migration rate(Barton & Wilson 1995), we first assessed the influence ofsegment size on parameter estimation. To this aim, wequalitatively contrasted four models differing in theirsegment size. For each model, we performed 1 millionsimulations. The tested segment sizes were 200 m (4348segments), 1000 m (868 segments), and 2000 m (434segments). Furthermore, we performed simulations withsegment size considered as a free parameter. In that case,for each simulation, a segment size was drawn from auniform distribution between 200 m and 2000 m. Then,we estimated for each of the four models the posteriorprobability using the rejection-sampling method (Pritchardet al. 1999; Estoup et al. 2004). This method estimates theposterior probability of different models by comparing thedistribution of distances between observed and simulatedsummary statistics obtained under various models. Underthe assumption that all models have the same priorprobabilities, the posterior probability of the i-th model,Pr(Mi), is simply obtained by ranking simulations accordingto their associated distances, and counting how manysimulations (ni) from the i-th model are found among somed smallest distances, and simply calculating Pr(Mi) as ni/d(Pritchard et al. 1999; Estoup et al. 2004). Here, we reportthe posterior probabilities for d = 100 closest simulations.We then used the best model (i.e. that with a segment size of2000 m) for later simulations. The computations necessaryto explore these four models required more than 630processor days on 2.4 GHz Linux machines.
Forward demographic simulation
During the last ice age, most parts of the Swiss river basinswere covered by ice, implying that Swiss rivers have beencolonized by fish located in open water after the LGM(Fig. 1). We used as the colonization source, a fish populationof 100 individuals. For each simulation, its geographicallocation was randomly assigned to an ice-free segment.
After the onset of the colonization, the range of thepopulation increased progressively because of a combinationof local demographic growth and migrations from theoccupied demes to empty neighbours. At each generation,the density of each deme is first logistically regulated, as
where Nt is the population density of the deme at time t, Kis the carrying capacity (i.e. the maximum number ofindividuals that can be sustained by local resources), and ris the intrinsic rate of growth per generation. Note that wesimulate haploid individuals, and N and K are thereforeexpressed in number of effective genes. The demographicregulation step is followed by a migration step. In order totake into account the fact that individuals may migrate atdifferent rates during the colonization phase and when thehabitat is already fully colonized (Saether et al. 1999), weintroduced a density-dependent migration rate m thatchanges smoothly between low local densities mLow andhigh local densities mHigh (Fig. 2). We used the generalizedlogistic curve (Richards 1959) computed as
where Dt is the current occupancy index defined as Nt/K.A defines the smoothness of the slope, and it is arbitrarilyset to 1000 (Neuenschwander 2006a). If mLow is larger thanmHigh, the migration rate is higher during colonization,while the reverse is true if mLow is smaller than mHigh. If thetwo migration rates mLow and mHigh are equal, then themigration rate is constant at all densities. The number ofemigrants is then distributed among the neighbouring
demes taking their densities into account. The probabilityPi of a given emigrant to go to deme i is calculated as
Neighbouring demes with low density will thereforereceive more immigrants than those with high densities,which seems realistic (Downhoer et al. 1990). The effectivenumber Mi of emigrants that reach a neighbouring deme iis thus
At each generation, the density and the number ofimmigrants for all demes are stored in a database, which isused later during the genetic backward simulations. Foreach model, we performed 1 million simulations of thecolonization process.
Backward genetic simulations
In the second phase of the simulation framework, a spatiallyexplicit coalescent approach is used to simulate the genegenealogies and the molecular diversity at microsatelliteloci for an arbitrary number of samples. Since the geneticvariation of samples of neutral genes depends only on thedemography of the demes (e.g. Hudson 1990; Nordborg2001), we can generate genetic data using the densities andmigration rates simulated during the demographic phase.Starting at the present generation, the genealogy of sampledgenes is simulated backward in time. At each generation,genes have a given probability to move to a different demeor to coalesce with another gene present in the same deme.The coalescent process at a given locus stops when a singlegene lineage remains, which means that the most recentcommon ancestor (MRCA) of the sampled genes has beenreached. A genealogy is computed independently for eachmicrosatellite locus. Mutations are then sprinkled over thebranches of the genealogy according to a Poisson processwith rate μit, where μi is the mutation rate of locus i and tis the length (in generations) of a given branch. We assumedthat the locus-specific mutation rates μi follow a gammadistribution with mean equal to μ (Gamma (μ, 2/μ),Excoffier et al. 2005). We used a multistep mutation model(generalized stepwise mutation model, GSM, Zhivotovskyet al. 1997; Estoup et al. 2002), where the coefficient parameterPi characterizes the geometric distribution of the number ofrepeats by which new mutant alleles differ from the ancestralallele. The individual coefficient parameters Pi wereassumed to follow a beta distribution around P, defined asfollows: if P ≥ 0.001 then Pi = Beta(a, b) with a = 0.5 + 199Pand b = a(1−P)/P, otherwise Pi = 0 (Excoffier et al. 2005).
N N rK N
Kt tt
+ = +−⎡
⎣⎢⎤⎦⎥1 1 ,
Fig. 2 Density dependent migration rates implemented in ourmodel. In this example the migration rate at low densities mLow ishigher than the migration rate at high densities mHigh, implyingthat the migration rates are larger during the colonization phasethan when the habitat is already fully colonized.
m mm m
AHigh Dt= +
−
+ −Low High
1 1 2,
P K NKNi i i
n
nn
Neighbours
=⎡
⎣⎢⎢
⎤
⎦⎥⎥=
−
∑1
1
M P N mi i i= .
A B C E S T I M AT I O N O F C O L O N I Z AT I O N H I S TO RY 5
We used the Approximate Bayesian Computation (ABC)approach (Beaumont et al. 2002) to estimate the parametersof our model, which are listed in Table 1. The ABC approachallows one to obtain the posterior distributions of theparameters, from which estimators and credible intervalscan be extracted. One major advantage of the Bayesianinference is that some prior information can be incorporatedinto the estimation procedure. Our ABC procedure consistedof three steps: First, we simulated a large set of 1 milliongenetic data sets using the same sample sizes and locationsas the observed data set and using parameter valuesrandomly drawn from previously defined prior distributions(see Table 1). In the second step, the simulated data sets arecompared to the observed data set on the basis of summarystatistics (see below), and the 5000 simulations that areclosest to the observations are retained, while the others arerejected. In the third step, the posterior distributions ofparameters of interest are estimated by a multiple andlocally weighted linear regression of the summary statisticscomputed from the retained simulations (Beaumont et al.2002). In this study, we computed the posterior distributionsaccording to equation 9 of Beaumont et al. (2002). We usedthe median of this posterior distribution as the point estimatebecause it performed best compared to other point esti-mates (mode, mean, and y-intercept of the regression linecomputed on the target summary statistics) when weassessed the accuracy of the estimation (unpublished data).Following Excoffier et al. (2005), all parameters weretransformed before the regression as y = log[tan(x)–1] inorder to ensure that the posterior distributions of theparameters remain within the prior ranges.
Prior distributions
Prior distributions of the parameters were defined usingavailable information from the literature.
1 Colonization time T: the colonization of the Swiss Rhinebasin necessarily occurred after the LGM (Jäckli 1962;Hantke 1980), when most of Switzerland was covered byice. Therefore, the onset of the colonization T was drawnfrom a uniform distribution between 600 and 7000generations, which corresponds approximately to theperiod between 2100 and 24 500 years bp, assuming ageneration time of 3.5 years for C. gobio (Mills & Mann1983).
2 Carrying capacity K: the observed population densityof C. gobio is generally highly variable between rivers.Utzinger et al. (1998) recorded census population sizes of0.002–0.4 individuals per square metre in Swiss alpinerivers, and higher densities of up to one individual persquare metre were recorded for lowland streams (Crispet al. 1974; Mills & Mann 1983). To cover this broad range,we used a log uniform prior for K ranging from 20 to5000 genes (haploid individuals) per kilometre, whichcorresponds to a density between 10 and 2500 diploidfish. A log uniform distribution allowed us to put moreweight on small values.
3 Migration rate m: the occurrence of highly differentiatedpopulations in the Rhine river system (Vonlanthen et al.2007) suggests that C. gobio did not migrate much in therecent past and/or that its effective populations were ingeneral small. On the other hand, C. gobio is also widelydistributed in headwaters, suggesting that the species isable to rapidly colonize new habitats. Downhoer et al.
Table 1 Prior distributions of the model parameters. For each prior distribution the mean, mode and several quantiles [0% (min), 5%, 50%(median), 95%, 100% (max)] are given
Notes: T is the onset of expansion, K is the carrying capacity, r is the growth rate, mLow is the migration rate at low density, mHigh is the migration rate at high density, μ is the mean mutation rate, and P is the mean geometric parameter. The estimates are based on 1 million simulated data sets. IR stands for irrelevant. †T is measured in number of generations. ‡K is measured in number of effective genes (haploid individuals).
(1990) indeed showed that the migration distance ofC. gobio is negatively density dependent with larger dis-tances favoured at lower densities. Therefore, we set theprior of the migration rate at high densities (in occupiedhabitats) mHigh to a uniform prior between 0 and 0.5 andthe prior of migration rate at low densities (during thecolonization phase) mLow was assumed to follow auniform distribution between 0.3 and 0.8. This reflectsour belief that migration was potentially higher duringcolonization than at carrying capacity.
4 Growth rate r: the prior of the local population growthrate r follows a beta distribution B(4, 10) that is restrictedto a range between 0.0 and 3.0. The rate r mainly influencesthe speed of the colonization process, and it has littleimpact when population sizes approach carrying capacity.
5 Mutation rate μ: the prior of the average mutation rate μfollows a beta distribution B(7, 7) and is restricted to arange between 10–7 and 10–3.
6 Coefficient P of the geometric distribution by which anew mutant allele differs from its ancestor: the prior ofthe average geometric parameter P was set uniformlybetween 0 and 0.8.
Since genetic variation depends mainly on the combina-tion of demographic parameters and mutation or migrationprocesses, we also estimated some composite parameters,which were defined as: τ = 2 * T * μ, θ = 2 * K * μ, Km_Low =K * mLow, and Km_High = K*mHigh.
Summary statistics
The ABC method compares simulated and observed dataset on the basis of summary statistics. The summarystatistics were chosen in order to capture different aspectsof the data both at the within-population and at the between-population levels. We performed extensive simulationsto choose the most powerful set of summary statistics(Neuenschwander 2006b). This set consisted of eightsummary statistics. We computed the mean number ofalleles over loci k, the expected heterozygosity over loci H,and an average modified M statistic (Garza & Williamson2001; Excoffier et al. 2005) computed as
where L is the number of loci and R is the allelic range overloci, that is the difference in repeat number between thelargest and the smallest allele. These summary statistics,computed for each population separately, were summarizedby their means, resulting in three summary statisticsMean_k, Mean_H, Mean_M. The same summary statisticswere also computed for a global population made up of all
pooled samples, resulting in the three summary statisticsGlobal_k, Global_H, Global_M. Furthermore, we computedthe global coefficient of differentiation FST (Weir &Cockerham 1984) and the correlation coefficient ρbetween geographical distances D and pairwise FSTcomputed as
where , ,and . The coefficient ρ was computedbetween population three (Fig. 1) and all other populations.The choice of using this statistic was motivated by thefact that, in expanding populations, we expect to finda correlation between the genetic differentiation andthe distance from the onset of the expansion (similar toPrugnolle et al. 2005). Such a correlation (ρ = 0.669) betweenpopulation three (Fig. 1) and all other populations wasindeed observed in our C. gobio data set (Vonlanthen et al.2007).
Assessment of the quality of the estimation
We applied several methods to assess the quality of ourestimation. First, we assessed the potential for a parameter tobe correctly estimated using the coefficient of determination(i.e. the proportion of parameter variance explained by thesummary statistics, R2) computed across all simulations(Lefebvre 1983; Neuenschwander 2006b). This informationis important as some parameters of the model may nothave the same potential to be correctly estimated, and thisinformation is not included in the posterior distribution ofa parameter. Note that a similarity between the posteriorand the prior distributions does not necessarily imply thatthe parameter has a poor potential to be correctly estimated,as it could simply indicate that the prior distribution waswell chosen.
One advantage of the ABC method is that it is possible toassess the performance of the estimation with a reasonablysmall amount of extra computations. The principle is togenerate test data sets (1000 in this study) based on a fixedparameter set (‘true values’). Each of these test data sets isthen used as a pseudo-observed data set, from which theparameter values are known. Using the previously simu-lated data sets (1 million) based on the priors, it is possibleto estimate the parameters using the same ABC frameworkfor each test data set. The estimates can then be comparedto the parameter values (‘true values’) used to generate thetest data sets, which allows us to assess the quality of theestimation. The estimation of the parameters is performedlike in the real estimation under the ABC framework. Thequality of the point estimators were assessed by computingthe relative bias defined as
Mk
R
ll
L
ll
L=+
=
=
∑∑
1
11( )
,
ρ =SS
SS SSFD
FF DD
,
SS F DFD ii= ∑ − −( )( )ST STF D SS FFF i
= ∑ −( )ST STF 2
SS DDD i= ∑ −( )D 2
A B C E S T I M AT I O N O F C O L O N I Z AT I O N H I S TO RY 7
and the relative root mean square error defined as
where qi is the estimator of the parameter θ, and n is thenumber of test data sets (1000 in our case). If the relativebias is positive, the true value is overestimated, whereas itis underestimated when the relative bias is negative. Becausebias and RMSE are relative, a value of 1 corresponds to abias or RMSE equal to 100% of the true value. We alsocomputed the factor 2 statistic, defined as the proportion ofestimated values lying in an interval bounded by 50% and200% of the ‘true’ value. To further assess the quality ofthe posterior distributions, we computed their 50% and90% coverage, defined as the proportion of simulations inwhich the true value lies within the 50% (respectively 90%)credible interval around the estimate. Note that the factor2 gives information on the absolute precision of the estimator,which is not provided by the coverage properties.
Furthermore, we assessed the performance of theestimation by investigating how well simulated datasets based on the posterior distributions fit the observeddata set (Gelman et al. 2004). For this, 1000 genetic data setswere generated using parameters randomly drawn fromthe posterior distributions. The eight summary statisticsdefined above were used to compare the simulated data
sets with the observed data set. The fit of these simulateddata sets was compared to the fit of simulated genetic databased on the priors. The fit of the summary statisticswas also quantified using the relative bias and the relativeRMSE.
We finally evaluated the robustness of our estimationsby altering parameters and priors. In addition to the use ofdifferent segment sizes (see above), we investigated theeffect of an alternative prior for the growth rate r (using auniform instead of a log normal prior), and we modifiedour model to use only a single migration rate instead of twomigration rates for low and high densities. The robustnessof the estimates was quantified in the same way as theeffect of the segment size (see above), by comparing theposterior probabilities of models differing by their priors(Pritchard et al. 1999; Estoup et al. 2004) and the estimatesof the different models to the standard model.
Results
Importance of river segment size
The posterior probability of models differing by the use ofdifferent river segment sizes was computed. The modelwith a segment size of 2000 m is the best supported model(P = 0.42), while the model with a segment size of 200 m isthe model with lowest support (P = 0.12). The models witha segment size of 1000 m and that with a variable seg-ment size performed similar, obtaining both a posteriorprobability of 23%. This result is also illustrated by Fig. 3a
biasn
ii
n=−
=∑11
q θθ
,
RMSEn ii
n= −=∑1 1 21θ
θ( ) ,q
Fig. 3 Effect of different river segment size models on the fit to observed statistics, simulations and on parameters estimation. The left paneshows the Euclidean distance distributions obtained under different models. Unsuccessful simulations, that is simulations where thesampled locations were not colonized at the end of the demographic simulations, were assigned an infinite Euclidean distance. Therefore,the surface under the curves corresponds to the ratio of the number of successful to the total number of simulations (e.g. 0.608 for the 200-m model). The right pane shows the effect of the segment size on the estimation of the time of expansion (T).
showing the distribution of the Euclidean distances betweenobserved and simulated statistics for the four segment sizemodels. The high support for large segment sizes suggeststhat Cottus gobio subpopulations occupy a relatively largeportion of the rivers. In keeping with this result, simulationsperformed with a segment size of 200 m were not successfulin about 40% of the cases, since some sampled locations werenot colonized at the end of the demographic simulations.In contrast, the simulations performed under the other threemodels (1000 m, 2000 m, and variable size) were successfulin most cases (94–99%). It is also worth noting that thesethree models lead to very similar estimates of the expansiontime (T, Fig. 3b). Based on these results, we estimatedparameters based on simulations using a segment size of2000 m in the rest of the paper.
ABC estimations of the colonization history of Cottus gobio
The posterior distributions of all parameters are shown inFig. 4 and corresponding point estimates are listed inTable 2. The onset of the population expansion is estimatedat 1892 generations (median), which corresponds to 6600years bp assuming a generation time of 3.5 years for C. gobio(Mills & Mann 1983). The corresponding 90% credibleinterval ranged from 967 to 3017 generations (i.e. 3400–10 600years bp). The estimated carrying capacity K per deme is 65effective genes or 32 diploid individuals. The posterior
distribution of the carrying capacity is quite narrow asdemonstrated by the relative small 90% confidence intervalof 29–101 genes, compared to the prior range of 20–5000genes. The estimated migration rate at high densities mHighis 0.30 and the migration rate at low densities (i.e. duringcolonization) mLow is 0.42. However, the low coefficient ofdetermination of the estimate of mLow (R2 = 0.1%, Table 3)suggests that this estimation is not very reliable. The sameis true for the growth rate r which was estimated to be 0.61,but with a determination coefficient of only R2 = 0.3%. Theestimated mutation rate μ is 3.0 × 10–4 per generation andthe estimated geometric parameter P is as high as 0.65. Thecomposite parameters scaled by the mutation rate τ and θwere estimated as 1.55 and 0.05, respectively. The compositeparameters scaled by the migration rate KmLow and KmHigh,representing the number of immigrant genes per genera-tion at low and high densities, were estimated as 33 and26, respectively. A value of KmHigh = 26 corresponds tothe exchange of 13 diploid individuals per generationbetween adjacent populations at high populationdensities.
Performance evaluation
The performance of the estimation procedure is detailed inTable 3. On average, the estimated values deviate from thetrue values only by 18.0%. The parameter with the smallest
Table 2 Estimated parameters of the colonization history of the Rhine river system by Cottus gobio. Four point estimates are listed (median,mode, mean, and Regression, which is the y-intercept of the linear regression line computed on the target summary statistics), as well asthe 90% credible intervals
Parameters
Estimates90% credible intervalMedian Mode Mean Regression
Parameters scaled by the mutation rate:τ 1.55 1.51 1.57 1.80 1.07–2.04θ 0.05 0.05 0.05 0.06 0.03–0.07
Parameters scaled by the migration rate:KmLow 33 32 35 45 10–56KmHigh 26 27 25s 31 11–39
Notes: Estimations are based on 1 million simulations. Prior distribution are listed and defined in Table 1 and in the section method. The posterior distributions are shown on Fig. 4.
A B C E S T I M AT I O N O F C O L O N I Z AT I O N H I S TO RY 9
relative bias was θ (3.0%). The parameter KmHigh showedthe largest bias (41.6%), but this deviation represents only0.2% of the range of its prior distribution, which shows thatthis parameter is still quite well estimated. The relativeRMSEs were surprisingly small, ranging between 0.066and 0.431, suggesting that the ABC estimation procedure isoverall quite accurate. The width of the 50% and 90% credibleintervals are on average well estimated with averagevalues of 50% and 81%. However, we note that the coverage
of the single parameters may deviate strongly from theexpectation. Finally, except for mHigh and KmHigh, most ofthe parameter-estimated values lie within a factor 2 aroundthe ‘true’ value, and can therefore be considered as beingsatisfactorily estimated.
The coefficients of determination R2 of the parameters bysummary statistics are listed in Table 3. Parameters with aR2 below 10% are considered to be unreliably estimatedsince summary statistics explain little of their variability
Fig. 4 Posterior distribution (solid line) of the parameters of the colonization history by Cottus gobio. Posterior distributions were estimatedby performing a locally weighted regression adjustment (Beaumont et al. 2002) on 5000 retained simulations (out of 1 million). Priordistributions are shown as dashed lines.
(Neuenschwander 2006b). This is the case for r and mLow(R2 < 1%), as well as for mHigh (R2 < 10%). Also the R2 of themutation rate μ (R2 = 11.8%) was not much higher than ourthreshold of 10%, but the variability of the other parametersby the summary statistics is very satisfactory, with valuesbetween 59.0% and 87.2%. The composite parameters aregenerally better explained by the summary statistics thansimple parameters. For example, the R2 of τ is 79.3%, whileits constituent simple parameters had values of 11.8% (μ)and 59.0% (T). The same is true for the other three compositeparameters θ, KmLow, and KmHigh. It is noteworthy that eventhough the migration rate during colonization mLow has avery low associated R2, its composite parameters KmLow isassociated with a very high R2 of 82.6%. We thus considerthat our estimation procedure is able to correctly estimateparameters, except for r, mLow, and mHigh.
As expected, simulated data sets based on the posteriordistributions are considerably more similar to the observeddata set than those based on the prior distributions. Figure 5shows this fit for each summary statistic separately. The fit(relative RMSE) of simulations based on posterior distribu-tions of parameters is on average two times better than thatbased on the priors (Table 4).
Effect of alternative priors
In general, our estimates are very robust to changes in thepriors of the parameters. As previously shown (Fig. 3),different priors on river segment size did not have anoticeable effect on the estimated time of expansion(except for a very small segment size of 200 m), nor on anyother parameter including the migration rates (results notshown). The comparison of the standard model (with riversegment length of 2000 m) with two alternative modelswith different prior distributions for growth and migrationrates showed that the standard model is the overall bestsupported model, with a posterior probability of 38%.However, the model with a uniform prior distribution(instead of a log normal distribution) for the growth rate,and the model with a single migration rate (instead of twomigration rates), obtained only slightly worse posteriorprobabilities of 36% and 26%, respectively. The similarityof the three models shows that these alternative priorshave little impact on the fit of the data. This is also shownby the comparison of the estimates of the different models(Fig. 6). Despite the lower support of the alternative models,the values of the estimated parameters did not change
Table 3 Accuracy of the estimated parameters assessed by simulation of 1000 test data sets with known parameter values
The parameter set (‘true’ values) used to generate the test data sets are the point estimates (median) of the data set of Cottus gobio (see text). Note, that the ‘composite’ parameters (τ, θ, KmLow, and KmHigh) are computed as the product of their constituent parameters, and are thus not identical to the estimates.The bias and RMSE are expressed in relative units. Factor 2 is defined as the proportion of estimated values which are found in an interval bounded by 50% and 200% of the ’true’ value. Coverage is the proportion of simulation in which the true value is within the 50% (respectively 90%) credible interval around the estimate. The coefficient of determination (R2) characterizes how well the parameter is explained by the summary statistics.*The estimates of these three parameters are poorly explained by the eight summary statistics (R2 < 10%) and their estimates are thus considered as unreliable.
A B C E S T I M AT I O N O F C O L O N I Z AT I O N H I S TO RY 11
Fig. 5 Comparison of posterior (solid line) and prior (dashed line) distributions of summary statistics. The observed value is representedby a vertical line. Posterior distributions were obtained by performing 1000 simulations using parameter values drawn from their posteriordistributions.
Note: The measurements are based on 1000 simulations using parameters drawn from the prior distributions and the posterior distributions, respectively. Bias and RMSE are expressed in relative units. The corresponding distributions are shown in Fig. 5.
Table 4 Reproducibility of the observedsummary statistics using simulations basedon the posterior distributions (posterior)and simulations based on the priordistributions (prior), respectively. The fit isquantified by the relative bias and therelative RMSE. The second column containsthe observed summary statistics
considerably with alternative priors, except for the growthrate r and the carrying capacity K and its related compositeparameters KmHigh, and KmLow, which all have posteriordistributions shifted towards lower values. The posterior
distributions obtained under a model with a single migrationrate leads to migration rates similar to those estimates formHigh in the model with two migration rates. This resultsupports the view that different migration rates at low
Fig. 6 Effect of different priors on estimated parameters. The solid line represents the posteriors obtained under our standard model (riversegment of 2000 m, log normal prior for the growth rate, and two distinct migration rates, as on Fig. 4). The dashed lines represent theposteriors obtained under a model where a uniform prior is used for the growth rate. The dotted lines represent the posteriors under a modelwithout density-dependent migration. The upper left pane shows the Euclidean distance distributions under the three models. Because themodel without density-dependent migration has a single migration parameter (here shown in the plots of KmHigh and mHigh), the posteriordistributions for the parameters KmLow and mLow are not available for this model. The distributions are very similar for most of theparameters, and in particular for the onset of the colonization (T), but the growth rate and the carrying capacities show contrastingdistributions, like those of the corresponding composite parameters (KmLow and KmHigh).
A B C E S T I M AT I O N O F C O L O N I Z AT I O N H I S TO RY 13
densities (i.e. during colonization) have a small effect onthe simulated data. Nevertheless, the fit to the data isbetter under a model with density-dependent migrationrates.
Discussion
Colonization history of Cottus gobio
The expansion time for Cottus gobio, estimated around 6600years bp is surprisingly recent, which suggests that C. gobiocolonized the Swiss Rhine river system well after the endof the LGM. Even when considering the 90% credibleinterval of 3400–10 600 years bp, it seems that C. gobioexpanded on the Swiss Rhine river system later thanpreviously thought. This result is at odds with previousstudies on salmonid species (Steinmann 1951; Rubin 1990;Largiadèr et al. 1996) and with a recent C. gobio study(Vonlanthen et al. 2007) suggesting that the Swiss Rhineriver system had been recolonized in the wake of retreatingglaciers after the LGM starting some 20 000 years ago.Rather, our results suggest that the Swiss Rhine riversystem was colonized after the end of the Younger Dryassome 10 500 years ago. The Younger Dryas was a brief coldperiod lasting for about a thousand years with an averagedrop in yearly temperature of about 5 °C (Severinghauset al. 1998). This cold episode brought back glacial conditionsin high latitudes of the Northern Hemisphere, but the impactof the Younger Dryas was much weaker on the Swiss glaciers(Joerin et al. 2006). A possible explanation for this recentcolonization could be that better adapted populations ofC. gobio have replaced other populations leading to multiplecolonizations of the Rhine river system (Riffel & Schreiber1995; Nolte et al. 2005). Although this hypothesis lacksempirical support and warrants further investigation, itwould imply that our estimate of a late colonization timewould point to the last major colonization event. Further-more, we cannot exclude the possibility that temporalenvironmental fluctuations had a strong impact on thegenetic differentiation of C. gobio populations, which wedid not consider in our simulations. For example, thetemperature of rivers and lakes was obviously colderduring the retreat of the glaciers after the LGM than today,and it is thus possible that the ageing of the fish was slowerand thus the mean generation time longer than today(Mills & Mann 1983).
We estimated C. gobio density as about 32 individuals perkilometre of river segment. If we assume a mean riverwidth of 5 m, this density corresponds to 0.0064 effectiveindividuals per square metre. This density is low but inagreement with the estimation of Utzinger et al. (1998)reporting census population sizes of C. gobio in Swiss alpinerivers of 0.002–0.4 individuals per square metre. A previousanalysis of the genetic structure of C. gobio (Vonlanthen
et al. 2007) revealed highly structured populations in theSwiss Rhine basin, with an overall FST of 0.35. Based on ourown results, this high degree of population differentiationcan be therefore attributed to a very low effective populationsize. In agreement with this result, the analysis of 796 bp ofthe mitochondrial DNA revealed an even higher level ofgenetic structure (overall FST = 0.86) in the Swiss Rhinebasin (Reckeweg 2003). Such a high FST value was due tothe fact that many populations were fixed for a singlemtDNA haplotype. The present model unfortunatelydoes not allow us to reliably estimate the migration ratesdirectly, as their associated coefficients of determination(R2) are very low (Table 3). However, based on the scaledparameters KmLow and KmHigh, it seems that the colonizationability of C. gobio was similar to the migration ability atcarrying capacity.
The estimated mutation rate μ (3.0 * 10–4) is found close tothe assumed mean mutation rate for microsatellites (i.e.5 * 10–4, Estoup et al. 2002). The mean geometric parameter P(0.65) is however, found much larger than that previouslyreported (Estoup et al. 2002), suggesting that a large propor-tion (65%) of the mutations lead to nonsingle-step mutations,resulting in a broad distribution of allele lengths in thesample. However, this wide allele length distributioncould also be attributed to admixture between allopatricpopulations, suggesting that the Rhine basin had beencolonized by several source populations (Riffel & Schreiber1995; Nolte et al. 2005), a scenario that we have not con-sidered here. However, several studies have postulated theoccurrence of a temporary connection between the Rhineand Danube river systems in the area of Lake Constance(Volckaert et al. 2002; Behrmann-Godel et al. 2004) andbetween the Rhine and Rhône river systems in the area ofLake Geneva (Steinmann 1951; Rubin 1990; Largiadèr et al.1996; Vonlanthen et al. 2007) at the end of the LGM, allowingfish to migrate into the Rhine river system. Nevertheless,such possible secondary contact scenarios are difficult toreconcile with the recent estimated colonization time.
Quality of the estimates
This study shows that compared to common geneticanalyses (e.g. Vonlanthen et al. 2007), the ABC frameworkis a powerful tool not only to estimate parameters ofcomplex models but also to reliably estimate associatedcredible intervals. All, but three parameters (the migrationrates mLow and mHigh, and the growth rate r), have thepotential to be correctly estimated since the summarystatistics explain a substantial proportion of their variability.The unscaled parameters, such as the time of the expansion,the carrying capacity and the migration rates are less wellexplained by the summary statistics (they have smallerassociated R2) than the corresponding scaled (composite)parameters (τ, θ, KmLow, and KmHigh). The two parameters
being most difficult to estimate (R2 < 1%) belong to thedemographic model (growth rate and the migration rateduring colonization). Both parameters play an importantrole during the colonization phase and influence coloniza-tion speed, but they are not very important once thelocal densities approach their carrying capacity. Thus, thecolonization period may be too short to significantly affectthe current genetic variation, and therefore prevents acorrect estimation of these parameters. In contrast to themigration rate at low densities mLow, the migration rate athigh densities mHigh has a stronger influence on currentpatterns of genetic variation, and should therefore be betterestimated. The comparison of our model with alternativemodels (alternative priors) show that the estimates are veryrobust to changes in the model.
Outlook
Even though our inference procedure brings more informa-tion on the past demography of bullhead in the Swiss Rhinebasin than previous analyses (Vonlanthen et al. 2007), ourestimation procedure could certainly still be improved.First, the available genetic information could be expandedby adding additional populations and/or by increasingthe number of loci, which should improve the quality of theestimation (see Excoffier et al. 2005). Second, increasing thenumber of simulations should result in slightly more preciseestimates with narrower credible intervals (Excoffier et al.2005; Neuenschwander 2006b). Third, additional realismcould be incorporated in our model, which could potentiallylead to better estimations. In the present model, the riversystem is homogenous as all river segments have similarcharacteristics for a given simulation (e.g. same density perriver segment length). The use of environmental hetero-geneity (e.g. fish density that would depend on the altitude)could make the model more realistic (Wegmann et al. 2006),or one could also add information such as water velocity orwater temperature to the model. However, one of the mostchallenging and promising source of realism is the incor-poration of temporal heterogeneity, where several landscapecomponents change over time. For instance, the exactdynamic of the retreat of the glacier could be modelled tobetter understand its impact on current genetic variation.The availability of more reliable environmental informationat key time periods (e.g. LGM, Younger Dryas) wouldalso certainly be beneficial for any future studies on thecolonization history of various Alpine fish species. The in-corporation of such additional realism would raise thecomplexity of our model. However, as the ABC methodcan deal in principle with arbitrarily complex models, itsmain limitation is available computational power. It seemsreasonable to anticipate that higher computation powerwill soon be available, making the investigation of suchmore realistic models feasible in the near future.
Acknowledgements
We thank William Meikle and Peter Neuenschwander for commentson earlier versions of the manuscript. We have found the sugges-tions of three anonymous reviewers very useful to improveseveral aspects of this work. We are also grateful to Sabine Finkand Daniel Wegmann for valuable discussions and contributionsto this manuscript. This work was supported by Swiss NSF grantsno. 3100A0-100800 and 3100A0-112072 to LE.
References
Austerlitz F, Jung-Muller B, Godelle B, Gouyon P-H (1997)Evolution of coalescence times, genetic diversity and structureduring colonization. Theoretical Population Biology, 51, 148–164.
Avise JC (1998) The history and purview of phylogeography: apersonal reflection. Molecular Ecology, 7, 371–379.
Barton NH, Wilson I (1995) Genealogies and geography. PhilosophicalTransactions of the Royal Society of London. Series B, Biological Sciences,349, 49–59.
Beaumont MA (1999) Detecting population expansion and declineusing microsatellites. Genetics, 153, 2013–2029.
Beaumont MA, Zhang W, Balding DJ (2002) Approximate Bayesiancomputation in population genetics. Genetics, 162, 2025–2035.
Behrmann-Godel J, Gerlach G, Eckmann R (2004) Postglacialcolonization shows evidence for sympatric population splittingof Eurasian perch (Perca fluviatilis L.) in Lake Constance. MolecularEcology, 13, 491–497.
Bernatchez L (2001) The evolutionary history of brown trout(Salmo trutta L.) inferred from phylogeographic, nested clade,and mismatch analyses of mitochondrial DNA variation.Evolution, 55, 351–379.
Bernatchez L, Wilson CC (1998) Comparative phylogeographyof Nearctic and Palearctic fishes. Molecular Ecology, 7, 431–452.
Brunner PC, Douglas MR, Osinov A, Wilson CC, Bernatchez L(2001) Holarctic phylogeography of Arctic charr (Salvelinus alpinusL.) inferred from mitochondrial DNA sequences. Evolution, 55,573–586.
Crisp DT, Mann RHK, McCormack JC (1974) The populations offish at Cow Green, Upper Teesdale, before impoundment. Journalof Applied Ecology, 11, 969–996.
Currat M, Excoffier L (2004) Modern humans did not admix withNeanderthals during their range expansion into Europe. PLoSBiology, 2, e421.
Currat M, Ray N, Excoffier L (2004) splatche: a program to simulategenetic diversity taking into account environmental heterogeneity.Molecular Ecology Notes, 4, 139–142.
Douglas MR, Brunner PC (2002) Biodiversity of Central AlpineCoregonus (Salmoniformes): impact of one-hundred years ofmanagement. Ecological Applications, 12, 154–172.
Downhoer JF, Lejeune P, Gaudin P, Brown AJL (1990) Movementsof the chabot (Cottus gobio) in a small stream. Polskie ArchiwumHydrobiologii, 37, 119–126.
Englbrecht CC, Largiader CR, Hanfling B, Tautz D (1999) Isolationand characterization of polymorphic microsatellite loci in theEuropean bullhead Cottus gobio L-(Osteichthyes) and theirapplicability to related taxa. Molecular Ecology, 8, 1966–1969.
Englbrecht CC, Freyhof J, Nolte A et al. (2000) Phylogeography ofthe bullhead Cottus gobio (Pisces: Teleostei: Cottidae) suggests apre-Pleistocene origin of the major central European populations.Molecular Ecology, 9, 709–722.
A B C E S T I M AT I O N O F C O L O N I Z AT I O N H I S TO RY 15
ESRI (2005) ARCINFO. Environmental System Research Institute,Inc., Redlands, California.
Estoup A, Clegg SM (2003) Bayesian inferences on the recentisland colonization history by the bird Zosterops lateralis lateralis.Molecular Ecology, 12, 657–674.
Estoup A, Jarne P, Cornuet JM (2002) Homoplasy and mutationmodel at microsatellite loci and their consequences for populationgenetics analysis. Molecular Ecology, 11, 1591–1604.
Estoup A, Beaumont MA, Sennedot F, Moritz C, Cornuet JM(2004) Genetic analysis of complex demographic scenarios:spatially expanding populations of the cane toad, Bufo marinus.Evolution, 58, 2021–2036.
Excoffier L (2004) Patterns of DNA sequence diversity and geneticstructure after a range expansion: lessons from the infinite-island model. Molecular Ecology, 13, 853–864.
Excoffier L, Estoup A, Cornuet JM (2005) Bayesian analysis of anadmixture model with mutations and arbitrarily linked markers.Genetics, 169, 1727–1738.
Garza JC, Williamson EG (2001) Detection of reduction in populationsize using data from microsatellite loci. Molecular Ecology, 10, 305–318.
Gelman A, Carlin JB, Stern HS, Rubin DB (2004) Bayesian DataAnalysis. CRC Press, Boca Raton, Florida.
Hamilton G, Currat M, Ray N et al. (2005) Bayesian estimation ofrecent migration rates after a spatial expansion. Genetics, 170,409–417.
Hänfling B, Brandl R (1998) Genetic variability, population sizeand isolation of distinct populations in the freshwater fish Cottusgobio L. Molecular Ecology, 7, 1625–1632.
Hänfling B, Hellemans B, Volckaert FA, Carvalho GR (2002) Lateglacial history of the cold-adapted freshwater fish Cottusgobio, revealed by microsatellites. Molecular Ecology, 11, 1717–1729.
Hansen MM, Mensberg KLD, Berg S (1999) Postglacial recolon-ization patterns and genetic relationships among whitefish(Coregonus sp.) populations in Denmark, inferred from mito-chondrial DNA and microsatellite markers. Molecular Ecology, 8,239–252.
Hantke R (1980) Letzte Warmzeiten, Würm-Eiszeit, Eisabbau undNacheiszeit der Alpen-Nordseite vom Rhein- zum Rhone-System. Ott,Thun, Switzerland.
Hewitt G (2000) The genetic legacy of the Quaternary ice ages.Nature, 405, 907–913.
Hudson RR (1990) Gene genealogies and the coalescent process.In: Oxford Surveys in Evolutionary Biology (eds Futuyma D,Antonovics J), pp. 1–44. Oxford University Press, Oxford, UK.
Ibrahim KM, Nichols RA, Hewitt GM (1996) Spatial patterns ofgenetic variation generated by different forms of dispersal duringrange expansion. Heredity, 77, 282–291.
Jäckli H (1962) Die Vergletscherung der Schweiz im Würm-maximum. Ecologae Geologicae Helvetiae, 55, 385–294.
Joerin UE, Stocker TF, Schluchter C (2006) Multicentury glacierfluctuations in the Swiss Alps during the Holocene. Holocene, 16,697–704.
Kontula T, Vainola R (2001) Postglacial colonization of NorthernEurope by distinct phylogeographic lineages of the bullhead,Cottus gobio. Molecular Ecology, 10, 1983–2002.
Largiadèr CR, Scholl A (1996) Genetic introgression betweennative and introduced brown trout Salmo trutta L. populationsin the Rhône River Basin. Molecular Ecology, 5, 417–426.
Largiadèr CR, Scholl A, Guyomard R (1996) The role of naturaland artifical propagation on the genetic diversity of brown trout
(Salmo trutta L.) of the upper Rôhne drainage. In: Conservation ofEndangered Freshwater Fish in Europe (eds Kirchhofer A, Hefti D),pp. 181–197. Birkhäuser-Verlag, Basel, Switzerland.
Lefebvre J (1983) Introduction Aux Analyses Statistiques Multidimen-sionnelles, pp. 102–104. Masson, Paris.
Mann RHK (1971) Populations, growth and production of fish in4 small streams in southern England. Journal of Animal Ecology,40, 155–190.
Marjoram P, Molitor J, Plagnol V, Tavare S (2003) Markov chainMonte Carlo without likelihoods. Proceedings of the NationalAcademy of Sciences, USA, 100, 15324–15328.
Mills CA, Mann RHK (1983) The bullhead Cottus gobio, a versatileand successful fish. Annual Reports of the Freshwater BiologicalAssociation, 51, 76–88.
Neuenschwander S (2006a) aquasplatche: a program to simulategenetic diversity in populations living in linear habitats. MolecularEcology Notes, 6, 583–585.
Neuenschwander S (2006b) Reconstruction of the post-glacial re-colonization of the Swiss Alps by the bullhead (Cottus gobio) basedon spatially explicit computer simulations PhD Thesis, Universityof Bern, Bern, Switzerland.
Nolte AW, Freyhof J, Stemshorn KC, Tautz D (2005) An invasivelineage of sculpins, Cottus sp (Pisces, Teleostei) in the Rhine withnew habitat adaptations has originated from hybridizationbetween old phylogeographic groups. Proceedings of the RoyalSociety B: Biological Sciences, 272, 2379–2387.
Nordborg M (2001) Coalescent theory. In: Handbook of StatisticalGenetics (ed. Balding DJ). John Wiley & Sons, Ltd, Chichester, UK.
Potvin C, Bernatchez L (2001) Lacustrine spatial distributionof landlocked Atlantic salmon populations assessed acrossgenerations by multilocus individual assignment and mixed-stock analyses. Molecular Ecology, 10, 2375–2388.
Pritchard JK, Seielstad MT, Perez-Lezaun A, Feldman MW (1999)Population growth of human Y chromosomes: a study of Ychromosome microsatellites. Molecular Biology and Evolution, 16,1791–1798.
Prugnolle F, Manica A, Balloux F (2005) Geography predictsneutral genetic diversity of human populations. Current Biology,15, R159–R160.
Ray N, Currat M, Excoffier L (2003) Intra-deme molecular diversityin spatially expanding populations. Molecular Biology and Evolu-tion, 20, 76–86.
Ray N, Currat M, Berthier P, Excoffier L (2005) Recovering thegeographic origin of early modern humans by realistic andspatially explicit simulations. Genome Research., 15, 1161–1167.
Reckeweg G (2003) Postglacial colonisation history of the Lake Genevabasin by the bullhead Cottus gobio L. Diploma work, University ofBern, Bern, Switzerland.
Richards FJ (1959) A flexible growth function for empirical use.Journal of Experimental Botany, 10, 290–300.
Riffel M, Schreiber A (1995) Coarse-grained population structurein Central European sculpin (Cottus gobio L.): secondary contactor ongoing genetic drift? Journal of Zoological Systematics andEvolutionary Research, 33, 173–184.
Rubin JF (1990) Biologie de l’omble chevalier Salvelinus alpinus (L.)dans le Lac Léman (Suisse). University of Lausanne, Lausanne,Switzerland.
Saether BE, Engen S, Lande R (1999) Finite metapopulationmodels with density-dependent migration and stochastic localdynamics. Proceedings of the Royal Society B: Biological Sciences,266, 113–118.
Timing of abrupt climate change at the end of the YoungerDryas interval from thermally fractionated gases in polar ice.Nature, 391, 141–146.
Slechtova V, Bohlen J, Freyhof J, Persat H, Delmastro GB (2004)The Alps as barrier to dispersal in cold-adapted freshwaterfishes? Phylogeographic history and taxonomic status of thebullhead in the Adriatic freshwater drainage. Molecular Phylo-genetics and Evolution, 33, 225–239.
Steinmann P (1951) Monographie der schweizerischen Koregonen(Teil 2.). Swiss Journal of Hydrology, 13, 54–191.
Swisstopo (2001) Vector200, Swisstopo, Bern, Switzerland.Utzinger J, Roth C, Peter A (1998) Effects of environmental para-
meters on the distribution of bullhead Cottus gobio with particularconsideration of the effects of obstructions. Journal of AppliedEcology, 35, 882–892.
Volckaert FAM, Hanfling B, Hellemans B, Carvalho GR (2002)Timing of the population dynamics of bullhead Cottus gobio(Teleostei: Cottidae) during the Pleistocene. Journal of EvolutionaryBiology, 15, 930–944.
Vonlanthen P, Excoffier L, Bittner D et al. (2007) Genetic analysisof potential postglacial watershed crossings in CentralEurope by the bullhead (Cottus gobio L.). Molecular Ecology, 16,4572–4584.
Wegmann D, Currat M, Excoffier L (2006) Molecular diversity after
a range expansion in heterogeneous environments. Genetics,174, 2009–2020.
Weir BS, Cockerham CC (1984) Estimating F-statistics for theanalysis of population structure. Evolution, 38, 1358–1370.
Weiss S, Persat H, Eppe R, Schlotterer C, Uiblein F (2002) Complexpatterns of colonization and refugia revealed for Europeangrayling Thymallus thymallus, based on complete sequencing ofthe mitochondrial DNA control region. Molecular Ecology, 11,1393–1407.
Wilson IJ, Balding DJ (1998) Genealogical inference from micro-satellite data. Genetics, 150, 499–510.
Zhivotovsky LA, Feldman MW, Grishechkin SA (1997) Biasedmutations and microsatellite variation. Molecular Biology andEvolution, 14, 926–933.
This study forms part of Samuel Neuenschwander’s PhD thesiswork on the reconstruction of the post-glacial re-colonization historyof the Swiss Alp Rivers. Samuel Neuenschwander, Laurent Excoffier,Nicolas Ray, and Mathias Currat are population geneticists withinterest in estimating demographic parameters from genetic dataunder complex evolutionary models. Carlo Largiadèr and PascalVonlanthen are fish population geneticists.
