1. Introduction
The concept of effective population size (Ne) was introduced by Wright (Reference Wright1931, Reference Wright1938) and can be defined as the number of breeding individuals in an idealized population that would show the same amount of variation of allele frequencies under random genetic drift. This parameter is usually smaller than the absolute population size. It is a key parameter in conservation and management because it affects the degree to which a population can respond to selection. Ne influences the rate of loss of genetic diversity, the rate of fixation of deleterious alleles and the efficiency of natural selection at maintaining beneficial alleles (Berthier et al., Reference Berthier, Beaumont, Cornuet and Luikart2002). If Ne declines too far, the loss of genetic variation resulting from genetic drift may put species or populations at risk of extinction by losing the raw material on which selection can operate.
Demographic data from mating systems or genetic data can be used to infer Ne. Unfortunately demographic data are generally difficult to collect in many wild populations and life-history information is often unavailable with sufficient precision to make a good estimate of Ne (Frankham, Reference Frankham1995). Genetic methods are more widely used. They are increasingly being developed by applying polymorphic molecular markers to resolve taxonomic problems and describe evolutionary and demographic history of species and populations. Currently, many different methods are available to infer Ne from genetic data in one or multiple samples. The most widely used estimator consists of measuring the variance of allele frequencies between generations (Tallmon et al., Reference Tallmon, Luikart and Beaumont2004). However, several studies noted that it is often biased towards high values (Luikart et al., Reference Luikart, Cornuet and Allendorf1999; Wang, Reference Wang2001; Berthier et al., Reference Berthier, Beaumont, Cornuet and Luikart2002). Phylogenetic methods are more efficient than non-phylogenetic methods (Felsenstein, Reference Felsenstein1992), primarily due to additional information provided by the tree structure. Genealogical modelling has greatly facilitated the estimation of demographic and mutational parameters using length variants microsatellite data (Chakraborty & Kimmel, Reference Chakraborty, Kimmel, Goldstein and Schlötterer1999; Feldman et al., Reference Feldman, Kumm, Pritchard, Goldstein and Schlötterer1999; King et al., Reference King, Kimmel and Chakraborty2000 by Storz & Beaumont, Reference Storz and Beaumont2002). In the past, most studies of genetic variation had used protein electrophoresis (Guyomard, Reference Guyomard, Gueguen and Prouzet1994). However, protein polymorphism has a limited potential for discrimination (Norris et al., Reference Norris, Bradley and Cunningham1999).
In this paper, we focus on highly variable microsatellite loci due to their power to detect genetic structures within and among populations. Furthermore, their assumed neutrality avoids the effects of selection, and allows the standard coalescence model to be used to estimate effective population sizes from the distribution of allelic frequencies. Effects of migration are integrated in some models but make the estimation more complicated. In the framework of the coalescence approach, a tree linking the alleles up to their common ancestor describes the relationships among alleles. A coalescence event appears each time two lineages in the tree join into a common ancestor, and the intervals between such events have a distribution that depends on Ne (Kuhner et al., Reference Kuhner, Yamato and Felsenstein1995). The usual approach for estimating this parameter is to compare statistics calculated from empirical datasets with a distribution generated by Monte Carlo simulations of the coalescent process (Wilson & Balding, Reference Wilson and Balding1998; Storz & Beaumont, Reference Storz and Beaumont2002).
Here, we compare estimates of Ne based on the coalescence model, and the performances of such methods based either on a single sample (MSVAR (Beaumont, Reference Beaumont1999) and DIYABC (Cornuet et al., Reference Cornuet, Santos, Beaumont, Robert, Marin, Balding, Guillemaud and Estoup2008)) or on two samples (TM3; Berthier et al., Reference Berthier, Beaumont, Cornuet and Luikart2002). Hence, we also compare a likelihood approach (TM3 and MSVAR) and Approximate Bayesian Computation (DIYABC). Comparative studies on real data that gauge the relative performances of Ne methods are few. Salmonids are an ideal species for assessing population structure's influence on Ne estimations. They show a propensity for structuring into genetically distinct populations due to the combined effects of their rearing habitat and their homing behaviour (Stabell, Reference Stabell1984). Here, we use four wild anadromous (i.e. adults migrate from the sea to breed in freshwater) European Atlantic salmon populations from north-west France (Rivers Oir and Scorff) and north-east Scotland (Rivers Spey and Shin) (Fig. 1). Atlantic salmon are subject to many pressures in Europe, including physical barriers to migration, exploitation by net and rod fisheries, pollution, the introduction of non-native salmon stocks, physical degradation of spawning and nursery habitat, and increased marine mortality. During the last 30 years, the decline of wild salmon on both sides of the North Atlantic (Parrish et al., Reference Parrish, Behnke, Gephard, McCormick and Reeves1998; Jonsson & Jonsson, Reference Jonsson and Jonsson2004) has affected populations to differing degrees (Hawkins, Reference Hawkins and Mills2000). The four populations studied are pressured by different factors and are therefore subject to varying conservation and management strategies. Because their characteristics are well understood (Baglinière & Champigneulle, Reference Baglinière and Champigneulle1986; Baglinière et al., Reference Baglinière, Marchand and Vauclin2005; Butler, Reference Butler2004; Butler et al., Reference Butler, Middlemas, McKelvey, McMyn, Leyshon, Walker, Thompson, Boyd, Duck, Armstrong, Graham and Baxter2008), and they have large differences in abundance, they provide a useful opportunity to design tools for estimating Ne.
Fig. 1. Geographic locations of the four wild Atlantic salmon populations: Rivers Oir and Scorff in north-west France and the Rivers Shin and Spey in north-east Scotland.
For each population, fish from two samples have been genotyped with 37 microsatellite markers chosen for their high genetic quality (Nikolic et al., Reference Nikolic, Fève, Chevalet, Høyheim and Riquet2009). Based on a general analysis of their present genetic variability and structure, and on different subsets of markers, this paper provides an overview on the performance and sensitivity of three different estimators of Ne applied to declining populations of different sizes.
2. Materials and methods
(i). Current status and management
Atlantic salmon provide highly valued ecosystem services. They support rod and net fisheries (e.g. Butler et al., Reference Butler, Radford, Riddington and Laughton2009), and provide source stock for an aquaculture industry that produces over a hundred million Atlantic salmon, and whose biomass already exceeds that of wild populations (Gross, Reference Gross1998). However, wild Atlantic salmon are considered an endangered species. As a consequence of contracting range, a decline in abundance and the modification of demographic structure of adult populations (Anonymous, 2001, 2003; Caron & Fontaine, Reference Caron and Fontaine2003), the species has been placed on the Red List of threatened species in Europe (Porcher & Baglinière, Reference Porcher, Baglinière, Keith and Allardi2001). Recent studies (e.g. McGinnity et al., Reference McGinnity, Prodohl, Ferguson, Hynes, O'Maoileidigh, Baker, Cotter, O'Hea, Cooke, Rogan, Taggart and Cross2003; Finnegan & Stevens, Reference Finnegan and Stevens2008) have shown that the introduction of non-native salmon stocks can result in the alteration of genetic structure of the recipient populations, posing another potential threat to population fitness.
Atlantic salmon conservation and management is complicated by the species' life history. Juvenile salmon spend 1–4 years in freshwater depending on the growth conditions and the latitudinal position of the river (Baglinière, Reference Baglinière1976). They then migrate to the ocean and return as adults to spawn after 1 year (‘grilse’) or multiple years (‘Multi-sea winter’ or MSW) (Klemetsen et al., Reference Klemetsen, Amundsen, Dempson, Jonsson, Jonsson, O‘Connell and Mortensen2003). The Atlantic salmon is a good example of a highly migratory species with complex spatial and temporal patterns demonstrating significant local adaptations (Taylor, Reference Taylor1991), homing behaviour (Hansen & Jonsson, Reference Hansen and Jonsson1994; Saglio, Reference Saglio, Gueguen and Prouzet1994) and reproductive strategies (Fleming, Reference Fleming1996).
Table 1 shows the differences in habitat, census population size, fishery pressure, and conservation and management regimes for the Oir, Scorff, Spey and Shin. Census population sizes were estimated by extrapolating from in-river rod catches using known exploitation rates. The Scorff, Shin and Spey are coastal catchments, whereas the Oir is an estuarine tributary of the river Sélune and the most productive area for salmon within the Sélune catchment. Salmon in the Spey are little affected by artificial barriers to migration, whereas the Oir, Scorff and Shin are more affected. While water quality in the Oir and Scorff is impacted by human pressures (mainly agricultural practices and urban run-off in the estuaries), quality in the Spey and Shin is not acutely affected by anthropogenic activities. The Spey and Shin juvenile populations predominately migrate to sea after 2–4 years, and returning adults include both grilse and MSW fish. In the Oir and Scorff, most juveniles migrate after 1 year, and most returning adults are grilse. Exploitation rates by net and rod fisheries are low (10–20%) for the Shin and Spey, but high (30–50%) for the Oir and Scorff. According to the International Union for the Conservation of Nature's (IUCN, 1994) classification of conservation status, the Spey and Shin populations are of ‘minor concern’, whereas the Oir and Scorff populations are ‘vulnerable’.
Table 1. Biotic and abiotic characteristics of the four studied salmon populations and their catchments
a Spey samples were taken from an upper catchment sub-stock of unknown size, estimated to be less than the minimum total adult population.
NA=not available.
There is a history of stock enhancement in all of the rivers. Hatchery-reared juvenile salmon of unknown Scottish origin have been planted in the Scorff during 1973–1979 (Baglinière, Reference Baglinière1979). In the Oir, some returning adults might have originated from hatchery-reared smolts sourced from native Sélune broodstock in the 1990s. In the Spey, native wild broodstock have been used to source ova and juveniles for enhancing production in areas above impassable obstacles since the 1970s. In the Shin, a hatchery program was established to mitigate the effects of the installation of hydroelectricity dams in the 1950s, based on native broodstock. Escaped farmed salmon occasionally enter the Spey and Shin in small numbers.
The Spey supports one of the largest Atlantic salmon populations in Scotland, with at least 60 000 adults entering the river annually (Butler, Reference Butler2004). Because of the relatively pristine nature of the river's habitat, and the status of the salmon population, the Spey was designated a Special Area of Conservation under the EU Habitats Directive (Council Directive 92·43/EEC) in 1999, with Atlantic salmon as a qualifying species. The objectives of special area of conservations are to avoid deterioration of the habitats of qualifying species or significant disturbance to those species, ensuring that the integrity of the site is maintained and it achieves favourable conservation status of the qualifying features (Anonymous, 2000).
(ii). Samples and DNA extraction
A total of 367 wild adult anadromous Atlantic salmon were sampled during the spawning migration in 2005 and 1988 for the Scorff, Oir and Spey, and in 2005 and 1992 for the Shin. For the Scorff, Oir and Shin, it was assumed that the samples were representative of the annual adult spawning populations, which are small (100–4000 fish; Table 1). However, in the Spey samples were taken from the upper catchment, where spring-running fish are known to originate (Laughton, Reference Laughton1991). Larger populations in rivers such as the Spey are known to contain genetically distinct population units (‘sub-stocks’), which differ in the timing of their return migration (Stewart et al., Reference Stewart, Smith and Youngson2002; Jordan et al., Reference Jordan, Cross, Crozier, Ferguson, Galvin, Hurrell, McGinnity, Martin, Moffett, Price, Youngson and Verspoor2005). Consequently, it was assumed that the Spey samples were taken from a sub-stock of spring-running fish of unknown size, but of less than 60 000. In all rivers, the individuals came from the same cohort (Table 2) to avoid biases in effective population sizes (Jorde & Ryman, Reference Jorde and Ryman1996).
Table 2. Demographic and genetic parameters of the four studied salmon populations
(FIS, FST, FIT) Wright's F-statistics; (Nm) number of migrants per generation with the private allele method of Barton & Slatkin (Reference Barton and Slatkin1986) ; (Hnb) unbiased expected heterozygosity; (Hobs) observed heterozygosity; (Nall) mean numbers of alleles per locus; (Tf) estimation of time (in years) since population started to decline with the interval at 95%; (Ancestral Ne) estimation of median ancestral population size with the interval at 95%; (Current Ne) estimation of median current population size with the interval at 95%; (a) MSVAR method; (b) DIYABC method; (c) TM3 method.
Sample sizes from each river ranged between 89 and 96 individuals. Pectoral or caudal fin clips were conserved in 95% ethanol, and scale samples were placed and dried in paper envelopes. All samples were kept at ambient temperature. Genomic DNA from fin and scale samples was extracted by boiling samples in 230 μl solution (proteinase K, TE buffer (tris/EDTA) and chelex) at 55°C for several hours, with a final period of 105°C for 15 min (Estoup et al., Reference Estoup, Largiadèr, Perrot and Chourrout1996; S. Launey, personal communication). After one night at 4°C, the supernatant was diluted in 200 μl of chelex and then stored at −20°C.
(iii). Genotyping
In the present study, a set of 37 salmon microsatellite loci from the Salmon Genome project (http://www.salmongenome.no/cgi-bin/sgp.cgi) previously identified by Nikolic et al. (Reference Nikolic, Fève, Chevalet, Høyheim and Riquet2009) were screened in all the individuals with the M13 labelling method (Schuelke, Reference Schuelke2000). PCR for all 367 individuals was handled in 10 μl within a 384 plate (TECAN 200): 1·5 mM MgCl2, 200 μM dNTPs, 0·1 μM forward primer, 0·15 μM reverse primer, 0·15 μM of M13-Fluo, 25–50 ng DNA and 0·5 unit Taq DNA polymerase. Precisions on the primers and amplification conditions for each marker are given in Nikolic et al. (Reference Nikolic, Fève, Chevalet, Høyheim and Riquet2009). A 2 μl volume of PCR product was added to 8 μl of deionized formamide and the internal size standard GENESCAN-400HD Rox (Applied Biosystems). Individual genotypes were obtained using ABI 3730 multi-capillary sequencer. Fluorescent DNA fragments size data were labelled by Genescan Analysis Software v3.7 (Applied Biosystems) to assign individuals by Genotyper 3·7 NT software (Applied Biosystems). From 10 random replicates, we evaluated genotyping error as relatively common (2·16%). From the initial set of 37 microsatellite markers, subsets of 28, 20, 10 and 5 markers have been selected according to their highest (H+ subsets) or lowest (H− subsets) observed heterozygosity, as estimated from the 367 individuals (Nikolic et al., Reference Nikolic, Fève, Chevalet, Høyheim and Riquet2009).
(iv). Historical processes
(a). Genetic diversity
The genetic structure of the populations was assessed from the complete genotypic dataset involving the 37 markers. Genetic diversity parameters were estimated using Genetix software version 4.05.2 (Belkhir et al., Reference Belkhir, Borsa, Goudet, Chikhi and Bonhomme1998) and GenAlEx 6 software (Peakall & Smouse, Reference Peakall and Smouse2006). Mean number of alleles per locus, actual heterozygosity and unbiased expected heterozygosity were calculated using the Hardy–Weinberg equilibrium (Nall, Hobs and Hnb, respectively). The inbreeding coefficient F IS per population was estimated with 10 000 bootstraps by Genetix software.
We developed a coalescent model to predict final genetic diversity when population size has undergone past changes, DemoDivMS (available at https://qgp.jouy.inra.fr/). Based on the Stepwise Mutation Model, this allows the present diversity at a microsatellite marker locus to be predicted as a function of present and past population size and mutation rate. Analytical calculations provide ‘Pk distributions’ (Shriver et al., Reference Shriver, Jin, Ferrell and Deka1997), i.e. the expected frequency of pairs of alleles that are alike in state, (i.e. the expected homozygosity at the locus) and the frequencies of pairs of alleles with any given difference of the numbers of the microsatellite motif (Chevalet & Nikolic, Reference Chevalet and Nikolic2010). This program includes an extension allowing for a continuous gene flow from a stable large population.
For the smaller populations (Oir and Scorff), we simulated the evolution of small populations (200–1000 individuals), derived from a larger population (10 000–50 000 effective size) 2000 to 4000 generations ago. We built a scenario in three steps assuming known ancestral and current effective sizes. Between the origin and the present day, we investigated the occurrence of a bottleneck, and checked various values of the effective sizes before and during the bottleneck, and of the times when population size changed and of the rate of immigration (Nm between 0 and 6). We assumed the Single Step Mutation Model and a mutation rate of 3×10−4 and 9×10−4.
(b). Structure and gene flow
Geographic and temporal genetic differentiation were estimated with R ST (Slatkin, Reference Slatkin1995) by GenAlEx software and with F ST according to Weir & Cockerham (Reference Weir and Cockerham1984) by Fstat 2.9.4 software (Goudet, Reference Goudet1995). Statistical significance was calculated from 10 000 permutations. Genetic distances (Nei, Reference Nei1978) were derived from allele frequencies. Ninety-five percent confidence intervals of the mean F ST, F IS and F IT estimates were obtained by bootstrapping (1000 replicates) over loci by Genetix software. Partition of genetic variance by Euclidean genetic distances among and within populations was calculated according to analysis of molecular variance (AMOVA) resulting in R ST estimations by 9999 permutations. The repartition of individual populations was graphically represented using factorial corresponding analysis (FCA) with GENETIX software.
Average effective numbers of migrants per generation (Nm) were derived using the four salmon populations sampled in 2005 and 1988/1992 by applying the private allele method (Barton & Slatkin, Reference Barton and Slatkin1986) using GENEPOP v3.4 (Raymond & Rousset, Reference Raymond and Rousset1995), and from F ST according to the relationship Nm=(1−F ST)/(4* F ST) using GENETIX v 4.03. To give a better indication of the migration rate per population, we used IM software (Nielsen & Wakeley, Reference Nielsen and Wakeley2001; Hey & Nielsen, Reference Hey and Nielsen2004) between populations within each country (Oir-Scorff and Shin-Spey), with 500 000 burning steps and 5 000 000 records period. Genetic assignment was performed using a Bayesian method (Rannala & Mountain, Reference Rannala and Mountain1997) as in GENECLASS 2·0 (Piry et al., Reference Piry, Alapetite, Cornuet, Paetkau, Baudouin and Estoup2004), and using the STRUCTURE software (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000) with 20 000 iterations and a burn-in of 5000.
(c). Evolutionary scenarios
Past evolution of populations was analysed using four algorithms and the 37 available markers using LAMARC (version 2.1, http://evolution.gs.washington.edu/lamarc, Kuhner, Reference Kuhner2006), MSVAR (Beaumont, Reference Beaumont1999), BOTTLENECK (Piry et al., Reference Piry, Luikart and Cornuet1999) and DIYABC (Cornuet et al., Reference Cornuet, Santos, Beaumont, Robert, Marin, Balding, Guillemaud and Estoup2008). The average mutation rate over loci (u), ancestral time (Tf) and ancestral Ne were estimated using MSVAR, running 2×107 steps per Markov chain Monte Carlo (MCMC) chain and DIYABC with 500 000 simulations. LAMARC was run only with the set of 37 markers to derive the growth rate (g) and provide the global evolution of populations (growing if g is positive, shrinking if g is negative). The g was determined by running 10 initial chains of 1000 steps each, discarding the first 500 and for 2 final chains of 10 000 steps discarding the first 1000. The recent effective size decline was tested using BOTTLENECK with the Two Phase Step Mutation model with 10 000 iterations. Evolutionary scenarios were compared using DIYABC software to derive the most likely outcomes.
(v). Estimators of current effective population sizes
(a). Current Ne from short-term temporal analysis
In each population, a temporal method was used to estimate the harmonic mean of current Ne between samples. The method was based on short-term allelic frequency changes between sampling periods (Oir, Scorff and Spey 1988 versus 2005; Shin 1992 versus 2005). Four generations were assumed to have elapsed between the 1988 and 2005 samples, and three generations between the 1992 and 2005 samples. The TM3 temporal method (Berthier et al., Reference Berthier, Beaumont, Cornuet and Luikart2002) was chosen for its higher efficiency compared to the classical F-statistic estimator (Nei & Tajima, Reference Nei and Tajima1981; Waples, Reference Waples1989), because it shows a narrower credible interval (CI) and greater accuracy when genetic drift is strong (Berthier et al., Reference Berthier, Beaumont, Cornuet and Luikart2002). The method was run using the eight subsets of markers and all the 37 markers with 500 000 iterations. Bayesian priors on the maximum size were based on demographic population data.
(b). Current Ne from long-term analysis
Two methods, DIYABC and MSVAR, were used to assess the current Ne from the distribution of alleles and coalescence processes. Long-term current Ne refers to the analysis with these two methods. For each population, both samples (1988/1992 and 2005) were analysed separately. DIYABC is an Approximate Bayesian Computation (ABC) that simulates data sets from priors, and then only data sets that are closest to the observed set are retained. The parameter values used to simulate these selected data provide an approximate posterior distribution by local linear regression. A second difference is the simulation of coalescence. Traditionally, it has been assumed that Ne is large enough to discount the probability that two or more coalescence events occurred in the same generation. However, population size can be very small and multiple coalescence events can occur in the same generation. The alternative is to reconstruct the lineages one generation at a time. DIYABC swaps between these two algorithms according to Ne. The last algorithm is taken when the effective size is very small or when the first algorithm overestimates the number of lineages (Cornuet et al., Reference Cornuet, Santos, Beaumont, Robert, Marin, Balding, Guillemaud and Estoup2008).
MSVAR version 1.3 was run assuming a linear trend of population size, and setting as prior a starting value compatible with demographic data. We assumed linear evolution because the average changes in population over long periods are more likely to be linear than exponential (Beaumont, Reference Beaumont1999). As for TM3, a Bayesian prior was set on the maximum size for DIYABC. Also, an Ne of 50 000 (which could represent a maximum ancestral size) was set for all four populations in order to evaluate the influence of prior information on posterior distributions.
3. Results
(i). Historical processes
(a). Genetic diversity
Inbreeding coefficient (F IS) quantifies the difference between observed and expected heterozygosity and evaluates the reduction of heterozygosity due to non-panmictic reproduction (Weir & Cockerham, Reference Weir and Cockerham1984). F IS was very low in all populations (Table 2) with the highest values for the Shin in 2005 and the Oir in 1988. Average numbers of alleles per locus (Nall) and observed heterozygosity (Hobs) across the four Atlantic salmon populations for the two sampling periods (1988/1992 and 2005), ranged from 11 (Oir in 1988 and Scorff in 1988) to 15 (Spey in 1988) and from 0·76 (Scorff) to 0·82 (Spey), respectively (Table 2). Nall was stable for Spey, increased for the Oir and Scorff and decreased for the Shin (Table 2). Unbiased expected heterozygosity (Hnb) were very close to Hobs in all samples. Hobs levels across the 37 microsatellites markers were very high with a maximum at 0·93 (SsaD144, BHMS331 and BHMS230) and few loci showed a significant deviation from Hardy–Weinberg equilibrium (Nikolic et al., Reference Nikolic, Fève, Chevalet, Høyheim and Riquet2009).
Our calculation of the expected evolution of diversity under variable population sizes by DemoDivMs (section 2.4.1) was qualitatively explained by a recent bottleneck 25–100 generations ago, assuming an effective size of 2500–5000 before the bottleneck. Even though results on differentiation and assignment results showed weak migration, we checked its impact on smaller populations. From IM values (0·004 and 0·005), the mean number Nm of immigrants in a population of size 200 is about 1. Accounting for possible migration was modelled by assuming immigrants come from a large population in equilibrium that is representative of the ancestral population of Atlantic salmons. The analysis indicates that the high observed heterozygosity cannot be maintained in the smallest populations (Ne=200 or 1000) by gene flow alone, except if it would take place at a higher rate (Nm>5) than observed. Assuming a recent bottleneck seems necessary to account for the observed diversity, even if the joint effects of genetic drift and gene flow are required, in such a way that the higher the migration rate is, the more ancient the bottleneck must be assumed. This suggested that this population was derived from a larger one that underwent a recent bottleneck.
(b). Structure and gene flow
Temporal genetic differentiations (R ST and F ST) were significant (P<0·05) only for the Scorff which had the highest values (0·022 and 0·009), indicating that this population underwent a significant genetic change during this period (Appendix A). The lowest values of R ST and F ST between samples were found for the Shin (0 and 0·001) (Appendix A). Similar trends were observed with Nei distance. We have not presented the pairwise R ST in Appendix A because they were of similar magnitude to the F ST. In spite of the low F ST values, the populations studied were well differentiated by the factorial correspondence analysis (Fig. 2) and using the unsupervised Bayesian approach in STRUCTURE. The set of 367 individuals was clearly split into four clusters indicating the presence of mixed stock introduction in the Oir 2005 sample and non-native introduction in the 1988 Scorff sample, which was absent in the 2005 sample.
Fig. 2. Factorial correspondence analysis of 367 wild salmon from the four populations (GENETIX software version 4.05.2, Belkhir et al., Reference Belkhir, Borsa, Goudet, Chikhi and Bonhomme1998). The doped circles represent the migrants detected by GENECLASS 2·0 (Piry et al., Reference Piry, Alapetite, Cornuet, Paetkau, Baudouin and Estoup2004).
An AMOVA (based on Euclidean R ST distances) across the four populations for the two sampling times revealed that the largest proportion of variation (90%) was found within populations. AMOVA between samples from the same populations showed the highest variance (2%) for the Scorff and lowest (0%) for the Shin, while the Oir and Spey were 1%.
The pairwise Nm estimated from the private allele method of Barton & Slatkin (Reference Barton and Slatkin1986) were lower than the pairwise from F ST of Weir & Cockerham (Reference Weir and Cockerham1984) (Table 3) with the highest (8·07) between Spey 1988 and Scorff 1988. The migration rates by IM were at 0·0030 and 0·0038 for Oir, 0·0047 and 0·038 for Scorff, 0·0026 and 0·0025 for Shin and 0·005 and 0·0067 for Spey, respectively, in 2005 and past samples. Four individuals were assigned outside their population of origin with GENECLASS 2·0: two from the Oir in 2005 were assigned to the Scorff 2005 population (probability (P) were equal to 0 and 0·002), one individual from the Oir in 1988 was assigned to the Spey 1988 population (P=0), and one individual from the Scorff in 1988 was assigned to the Spey 1988 population (P=0·001). These four migrants were all males and can be seen in the factorial correspondence analysis (Fig. 2, dotted circles).
Table 3. Pairwise numbers of migrant (Nm) from F ST of Weir & Cockerham (Reference Weir and Cockerham1984) by GENETIX software on the superior half matrix and number of migrants using private alleles (Barton & Slatkin, Reference Barton and Slatkin1986) by GENEPOP software on the inferior half matrix
(c). Evolutionary scenarios
Average mutation rates for the 37 microsatellites were estimated at 3×10−4 by MSVAR and slightly higher by DIYABC at 9×10−4. The first result is concordant with previous studies on Salmo salar by O'Reilly et al. (Reference O'Reilly, Herbinger and Wright1998). A negative posterior distribution of log10(r) values (with r equal Current Ne/Ancestral Ne) was revealed by MSVAR 0·4 and a negative growth rate (g) for overall populations by Lamarc software, which is consistent with a decline in effective sizes (Table 2). The various tests of heterozygosity deficit proposed in BOTTLENECK (Sign test and Wilcoxon test, TPM model) suggested a significant departure from constant population size in all populations.
Given an assumed generation time of 3–4 years, the estimated time since decline (Tf) ranged from 8000 to 20 000 years ago according to MSVAR 1·3, and around 13 000 years ago according to DIYABC. This was consistent for the four populations and for both sampling times (Table 2). The ancestral Ne values calculated by MSVAR 1·2 were approximately the same for all populations (approximately 50 000 individuals) suggesting a common ancestor. DIYABC revealed a smaller ancestral population size of 13 000–30 000 (Table 2) and two equally likely scenarios in which the populations were separated from a common ancestor by a star or cascade process from the southern (Scorff) to the northern population (Shin).
(ii). Estimators of effective population sizes
(a). Markers' number and heterozygosity variation
An increase in the markers' number and polymorphism led to a decrease in the variance of posterior distribution of current Ne estimates in almost all methods. While clearly shown for all populations with the MSVAR method, this phenomenon was less visible using TM3 and DIYABC methods for larger populations, and was absent for the Spey, the largest (c and d, Figs 4 and 6).
Using MSVAR (Fig. 3) the same trend was observed for the four rivers, showing an increase of mean values and a decrease of variances when the number and polymorphism of markers was increased. The TM3 (Fig. 4) method provided unexpected profiles with the subsets of five and ten loci (5 mH−, 5 mH+, 10 mH−, 10 mH+) for Oir and the subsets of five (5 mH−) for Scorff, generated large estimates of population size (Fig. 4 A and B). DIY ABC also generated large estimates for the Oir population size with the subset of low number of markers (Fig. 5 A).
Fig. 3. Posterior distribution of current effective population size estimators (log10 scale), for river Oir (2005 sampling), according to MSVAR. Curves refer to the eight subsets of markers (m) according to their lower heterozygosity (H−) in dotted line and their higher heterozygosity (H+) in solid line. The red curve refers to the full set of 37 markers. Similar results were obtained with other samples and rivers.
Fig. 4. Posterior distribution of current effective population size estimators (log10 scale) according to TM3, for the populations Oir (A), Scorff (B), Shin (C) and Spey (D). Curves refer to the eight subsets of markers (m) according to their lower heterozygosity (H−) in dotted line and their higher heterozygosity (H+) in solid line. The red curve refers to the full set of 37 markers. Similar results were obtained with other samples and rivers. See the legend in the corner of (A).
Fig. 5. Posterior distributions of current effective population size estimators (log10 scale) according to the DIYABC method for the populations Oir (A), Scorff (B), Shin (C) and Spey (D) in 2005. Similar results have been observed in the past samples (1988/1992). Curves refer to the eight subsets of markers (m) according to their lower heterozygosity (H−) in dotted line and their higher heterozygosity (H+) in solid line. The red curve refers to the full set of 37 markers. Similar results were obtained with other samples and rivers. See the legend in the left corner.
(b). Estimation when setting priors on the mean current Ne
Regarding the CIs and variance, short-term current Ne estimates using TM3 were more accurate for smaller populations (⩽1000) (Fig. 6 A and B) than for larger populations (⩾3000) (Fig. 6 C and D). Long-term current Ne estimates using both MSVAR and DIYABC were very consistent between sampling years, but their variances and 95% confidence intervals remained high for all samples. For example, the Spey estimate derived by MSVAR was 1076–50 364. Overall, current Ne for the four populations in 2005 were estimated to be 383 (Oir), 689 (Scorff), 304 (Shin) and 7344 (Spey) using MVAR, and 100 (Oir), 1174 (Scorff), 1842 (Shin) and 9417 (Spey) using DIYABC (Table 2). These values are nearer to the census population sizes (N) for the smaller populations (Oir and Scorff) than for the larger populations (Shin and Spey). However, if the variance around the median of Ne was considered the estimation of the effective size was near census size for all populations.
Fig. 6. Summary statistics of effective population size estimators (log10 scale) from the different methods (TM3, MSVAR and DIYABC) for the four populations: Oir (A), Scorff (B), Shin (C) and Spey (D). Census size is mentioned on the right y-axis. Median (horizontal black line), variance (box) and confidence interval at 95% (dotted black line on both sides of boxes) are given for each method. Priors are plotted as follows: Blue dotted lines: highest current Ne value for TM3 and DIYABC, starting value for MSVAR. Green dotted lines: lowest current Ne value for DIYABC.
Considering the differences in effective size values between methods, and taking into account higher ratios of Ne/N for smaller salmonid fish (Palstra & Ruzzante, Reference Palstra and Ruzzante2008), TM3 underestimated the Scorff population size (Fig. 6 B) and MSVAR underestimated the Shin population size (Fig. 6 C).
(c). Estimation when setting priors on the maximum current Ne
Changing the prior maximum size to 50 000 had an effect on point estimates of current Ne so that for all populations when running DIYABC, mean values were increased and the upper 95% limits of confidence intervals approached this maximum. On the contrary no effect was seen with TM3 for the smaller populations (Oir and Scorff). Setting the starting value at 50 000 using MSVAR had no appreciable effect on point estimates of current Ne in any population.
4. Discussion
Measurements of effective population size (Ne) are of importance in conservation and management because they give an overview on the evolution of genetic diversity. Effective size determines the rate at which genetic diversity is lost in the population by genetic drift (Franklin, Reference Franklin, Soule and Wilcox1980). The most genetically diverse populations are assumed to be ‘fitter’ (Ligoxygakis, Reference Ligoxygakis2001). A genetically viable population possesses the evolutionary legacy of the species and the genetic variability on which future evolutionary potential depends (Dodson et al., Reference Dodson, Gibson, Cunjak, Friedland, Garcia de Leaniz, Gross, Newbury, Nielsen, Power and Roy1998). Overfishing, blockage of migratory routes by hydroelectric dams and destruction of spawning habitat have severely depleted many wild salmon stocks (Waples, Reference Waples1990). Humans have altered natural ecosystems for many thousands of years, but the magnitude and rate of these changes have increased dramatically since the industrial revolution. Over recent decades, Atlantic salmon populations have declined or have been extirpated in many parts of its ancestral range (Parrish et al., Reference Parrish, Behnke, Gephard, McCormick and Reeves1998; Knockaert, Reference Knockaert2006) and measurements of Ne are few. Genetic estimates of effective size have never been conducted on the Oir, Scorff or Spey populations. Consuegra et al. (Reference Consuegra, Verspoor, Knox and De Leaniz2005) investigated Ne in the Shin, but based their estimate on four markers and did not provide estimates. These populations have management regimes and legislation in place to minimize the risks of further declines, but so far effective population size has not been implemented as a tool in their management strategies. Furthermore, the diversity of the rivers' habitats makes them interesting case studies of the effects of contrasting environments on their genetic structure (Whitlock & McCauley, Reference Whitlock and McCauley1999; Leberg, Reference Leberg2005; Wang, Reference Wang2005).
Before discussing the results on estimations of Ne obtained about the four salmon populations and their implication for conservation and management, it is useful to recall some features of these populations and to present the results on historical processes obtained by genetic analysis. Atlantic salmon populations are highly structured with strong genetic differentiation between regions and among rivers (Ståhl, Reference Ståhl, Ryman and Utter1987; Nielsen et al., Reference Nielsen, Hansen and Loeschcke1999; King et al., Reference King, Kalinowski, Schill, Spidle and Lubinski2001), with weaker differentiation within river systems (Nielsen et al., Reference Nielsen, Hansen and Loeschcke1999; Garant et al., Reference Garant, Dodson and Bernatchez2000; Primmer et al., Reference Primmer, Veselov, Zubchenko, Poututkin, Bakhmet and Koskinen2006; Vaha et al., Reference Vaha, Erkinaro, Niemela and Primmer2007, Reference Vaha, Erkinaro, Niemela and Primmer2008). This differentiation is partly explained by high fidelity to natal rivers (Stabell, Reference Stabell1984). Geographic genetic differences observed reflect the effects of underlying evolutionary forces, such as drift. Temporal (between samples from the same population) genetic differentiation (R ST, F ST and Nei distances) was significantly higher for the Scorff and lower for the Shin. These differences are explained in the Scorff by the introduction of non-native juveniles of unknown Scottish origin in the 1980s. If the Spey was the source, this last point could explain the highest values of numbers of migrants (Nm) between the Spey and Scorff in 1988. The evidence of reduced genetic variability in the Shin and of low temporal difference is probably from the result of long-term artificial stock enhancement program using fish of native (Shin) origin employed to mitigate the blocking of freshwater habitat by hydroelectricity dams in the 1950s. This native stock enhancement may explain the low temporal genetic distance. An effect of stock enhancement in the Oir is suggested by the STRUCTURE analysis. Significant differentiation among the Scorff samples, but not among the Oir samples can be explained by a brief and genetically distant gene flow into the Scorff from a Scottish source, whereas gene flow into the Oir has been more continuous and from a genetically similar French source, the Sélune. For the Spey, gene flow can be negligible given its large Ne. The clear differentiation between rivers and the high rate of assignment of individuals to their river of origin confirm that the populations maintained their originality.
The four populations also seem subject to a global decrease in wild salmon stocks. The genetic analysis revealed a negative growth rate and a bottleneck for all populations. Whereas allelic richness seemed to decrease in the Shin, increase in the Oir and Scorff and be stable in the Spey, genetic analysis on heterozygosity demonstrated a high and stable variability over time with no inbreeding for all samples. These results are concordant with previous studies that revealed high genetic diversity on wild populations of Atlantic salmon from Ireland and Norway with Nall and Hobs of 17·8 and 0·70, respectively, with 17 microsatellites (Norris et al., Reference Norris, Bradley and Cunningham1999) and from Scotland for the Shin population (sampled in 1987–1990) with a Hobs of 0·71 with 12 microsatellites (King et al., Reference King, Kalinowski, Schill, Spidle and Lubinski2001) and a bit lower with 4 microsatellites (Consuegra et al., Reference Consuegra, Verspoor, Knox and De Leaniz2005). The rate of loss of genetic diversity via genetic drift is higher in small populations and, in the absence of migration, this rate is expected to rise as the effective size decreases (Frankham et al., Reference Frankham, Ballou and Briscoe2002). Taking account of the low gene flow detected by IM for the smaller populations (0·003–0·004 for Oir and 0·004–0·005 for Scorff), we tested the bottleneck hypothesis by our new model DemoDivMs. Regarding the high genetic diversity (75–81%) detected in smaller populations (Oir and Scorff), we used it to evaluate the period when the bottleneck occurred. Our results suggest a common ancestor with a large size (10 000–50 000 effective individuals) dating back to the last glaciations, representing 2000 to 4000 generations. We assumed a scenario with an ancestral large population at great genetic diversity in equilibrium that splits into several populations. Trials with our method (using a mutation rate of 3×10−4 to 9×10−4 and the Single Step Mutation model for microsatellites) suggested that observed small populations (200–1000 effective size) could show the high observed heterozygosity if they derived from larger populations of size 2500–5000 that experienced a recent bottleneck 25–100 generations ago (i.e. 100–400 years), the bottleneck being all the more recent as gene flow occurs at a lower rate. Noteworthy, neither gene flow nor bottleneck alone could be considered as a single explanation, both phenomena seem necessary to explain the current status of these populations. However, the hypothesis of precocious parr contributing to the genetic variability in salmonid (Garcia-Vazquez et al., Reference Garcia-Vazquez, Moran, Martinez, Perez, de Gaudemar and Beall2001; Jones & Hutchings, Reference Jones and Hutchings2001, Reference Jones and Hutchings2002) cannot be neglected in view of the high proportions of fertilization of eggs by parr in Oir and Scorff (Baglinière & Maisse, Reference Baglinière and Maisse1985; Baglinière et al., Reference Baglinière, Maisse, Nihouarn, Gibson and Cuting1993).
Given the importance of identifying Ne for populations in decline, we compared recent coalescent methods to assess their sensitivity to the number of polymorphic markers, to the priors used in Bayesian approaches and to the structure of populations. Using different subsets of markers, we showed that methods providing estimates of current Ne may be more or less efficient. The efficiency of estimators is improved when the number of markers and their variability are increased. We chose markers which were independent, but violation of the independence among markers' loci should not affect the accuracy of estimation of current Ne (Wang & Whitlock, Reference Wang and Whitlock2003). The CIs and variance decrease by increasing the number and variability of markers with the three methods. However, this behaviour depended on populations, except for MSVAR. MSVAR results were improved on all populations by increasing the number of polymorphic markers (Fig. 4). The difference of accuracy between sets of markers (from low number and variability markers to high number and variability markers) was less accentuated for larger populations with TM3 and DIYABC (Figs 5 and 6). These methods presented more complex behaviour. They overestimate current Ne in smaller populations (Oir and Scorff) when using few markers. It appears that the used marker sets corresponded to uneven distributions of allelic frequencies. The sets corresponded to five and ten markers at low (H−) and high (H+) heterozygosity. The sets at low heterozygosity possess alleles at high frequencies around 60–80%. The sets at high heterozygosity (H+) did not show alleles with high frequency (>50%), but show uneven distributions between samples from the same population. The absence of alleles at intermediate frequency and high frequency of a single allele may lead to an overestimation of current Ne (Waples, Reference Waples1989). MSVAR did not seem to be impacted by the uneven distribution of allele frequencies in the smaller populations.
The incidence of Bayesian priors was weak using MSVAR, but quite strong using DIYABC for all populations, and TM3 for the largest populations (Spey and Shin). Hence, at least for these practical methods, Bayesian priors should be applied using reliable field data to avoid the derivation of inaccurate predictions.
Results from the set of 37 microsatellite markers demonstrated different sensitivity to the size of populations. The CIs and variance were greater for the larger populations (Shin and Spey) using TM3 (short-term Ne), large and homogeneous on all populations using MSVAR (long-term Ne) and heterogeneous using DIYABC (long-term Ne). The long-term current Ne estimates were strongly correlated between samples from the same population but CIs were larger for MSVAR than DIYABC except for the smallest population (Oir). Additional analysis using the classical temporal approach based on F-statistics (see equation (9) in Waples, Reference Waples1989) from NeEstimator v 1.3 (Peel et al., Reference Peel, Ovenden and Peel2004) demonstrated a high uncertainty in larger populations. The imprecision of temporal methods in large populations seems to occur with classical and coalescent models due to weak genetic drift. Berthier et al. (Reference Berthier, Beaumont, Cornuet and Luikart2002) compared TM3 with the classical F-statistical-based Ne estimator and showed narrower CIs and greater accuracy with TM3 when genetic drift was strong, but there was no improvement with weak genetic drift. Large CIs and variance observed in larger populations with TM3 can be explained by the ratio of sample size (S) to Ne. Too small a sample size may lead to insufficient extraction of genetic information in large populations. As pointed out by Nei & Tajima (Reference Nei and Tajima1981), precision increases with the ratio S/Ne. Waples (Reference Waples1989) demonstrated that a low ratio would translate into larger CIs, which means that populations with small Ne are most effectively studied. Furthermore, the effects of genetic drift in large populations might be swamped by sampling error using the temporal method, and if tS/Ne>√2, with t the sampling interval, increasing the number of loci (or alleles) has a greater effect on precision than increasing t or S (Waples, Reference Waples1989). This means that given the effective size to be estimated, minimum numbers of generations and sample sizes are required to achieve reasonable precision, unless a large number of loci can be surveyed.
The long-term methods (MSVAR and DIYABC) assumed that selection and migration were unimportant in changing population allelic frequencies relative to genetic drift. For the short-term current Ne estimates (TM3), all systematic forces (mutation, selection and migration) are assumed to be absent. The assumption of no mutation is reasonable because our sampling periods of 3–4 generations are sufficiently short to discount the effects of mutations. It may be reasonable to neglect the effects of selection because selection on most markers is unlikely to be strong enough to cause substantial changes in their frequencies (Wang & Whitlock, Reference Wang and Whitlock2003). Moreover, the loci used in our study have been tested to detect potential selection (Nikolic et al., Reference Nikolic, Fève, Chevalet, Høyheim and Riquet2009) and can be considered a reliable panel. In the context of this study, even if stock enhancement programs (based on non-native or native source stocks in Scorff and Shin, respectively) were associated with low values of Nm, they seem affected the temporal genetic differentiation, with stronger R ST, F ST and Nei distances in the Scorff and lower in the Shin. Since temporal estimators (TM3) interpreted allele frequency changes as due to genetic drift, they led to decreased estimated effective size in the Scorff, yielding a lower value than the long-term Ne. The effective size estimates derived from short-term methods on the Shin was higher than with the long-term methods because native reintroductions helped reduce the temporal genetic differences. Using model simulations of infinite source populations sending a constant supply of migrants with constant allele frequency Wang & Whitlock (Reference Wang and Whitlock2003) revealed that migration changes allele frequency more in the short term and less in the long term, leading to under and overestimation of current Ne, respectively. Regarding our results, introduction seems to have modified effective size estimators in different ways depending on the source. By comparing the effective size values between methods and by taking into account the results of higher ratio Ne/N on smaller salmonid fish by Palstra & Ruzzante (Reference Palstra and Ruzzante2008), TM3 seemed to underestimate the Scorff population size and MSVAR to underestimate the Shin population size. A native introduction, or from a genetically close source, could lead to underestimations of current Ne with long-term methods. A pulse of genetically different migrants coming between samples will underestimates current Ne with short-term methods. Even with small values of Nm, the temporal methods seemed more sensitive than expected to weak gene flow. In cases when artificial or natural gene flow is suspected, we suggest that both short-term and long-term estimates of population size should be evaluated.
In spite of these challenges, the current Ne estimates presented here are the first to be proposed for the French and the larger Scottish (Spey) studied populations, and provide potentially useful tool for their management. Efforts to establish lower population size thresholds have suggested current Ne=500 to maintain sufficient evolutionary potential (Franklin, Reference Franklin, Soule and Wilcox1980; Lande, Reference Lande1988; Franklin & Frankham, Reference Franklin and Frankham1998; Lynch & Lande, Reference Lynch and Lande1998). Estimating current Ne is thus an important tool in the assessment of genetic variability (Mace & Lande, Reference Mace and Lande1991). For the Oir, current Ne estimates of 383 and 100, compared with a 2005 estimated adult population of 130, suggesting that the river's IUCN (1994) conservation status of ‘vulnerable’ is still justified. The same is true for the Scorff where current Ne estimates are 689 and 1174, compared to a 2005 adult population of 1000. For the Shin, the 2005 adult population of 3000 exceeds the current Ne estimates of 304 and 1842. Similarly on the Spey current Ne estimates of 7344 and 9417 are probably exceeded by the 2005 population. In both cases, the IUCN (1994) status of ‘minor concern’ remains valid.
The discrepancies between current Ne estimators across the salmon populations contrasted here raise the question of which methods provide best estimates of Ne. All coalescent methods are potentially useful, but they may be biased because of their assumptions, particularly regarding migration. A short-term coalescent-based estimate, such as TM3, seems better suited to anadromous salmon populations than long-term estimates if the population is not large, if it does not undergo identified gene flow and if a high number of highly polymorphic markers are available. CIs become wider as current Ne values become larger. In the case of the long-term coalescence estimator, with a minimum sample of 50 individuals as recommended by Waples (Reference Waples1989), increasing the number of sampled individuals (S) may marginally improve the results because the variance of the estimator of θ (4Neμ) decays slowly, at a rate proportional to 1/log(S) (Deonier et al., Reference Deonier, Waterman and Tavaré2005). Hence, for larger populations, such as the Shin (up to 4000 individuals) and Spey (up to 60 000 individuals), the 37 markers may not be sufficient to accurately estimate Ne, and it would make little difference to increase sample sizes. For the short-term estimator, increasing sample sizes may be useful for larger populations. Using simulations, Ovenden et al. (Reference Ovenden, Peel, Street, Courtney, Hoyle, Peel and Podlich2007) estimated that 2000 sampled individuals are required to get a reliable estimate of an effective size of 8000 (Palstra & Ruzzante, Reference Palstra and Ruzzante2008). Another way to decrease CIs for short-term estimates (TM3) in the Shin and Spey populations would be to increase the number of generations between samples. Nevertheless, the amount of information is no longer simply proportional to sampling interval t, and too large a t may decrease the power of the method (Wang & Whitlock, Reference Wang and Whitlock2003). Usually, a short t is recommended (two to four generations) but if the migration rate is low, t could be larger to increase the estimation precision (Wang & Whitlock, Reference Wang and Whitlock2003). Finally, the temporal method (TM3) analysis is not equally efficient for all population sizes and other methods of estimating effective population size are necessary if current Ne is very large such as in the Shin and Spey, and if the influence of migration cannot be ignored. Regardless, the importance of selecting the appropriate tool for estimating current Ne is important for the conservation of wild salmon populations.
5. Conclusion
Despite their small size and declining status, French populations still show high levels of genetic diversity, similar to those found in the larger populations. Our results and coalescent simulations, where populations are assumed to derive from a common ancestor, suggest that a recent bottleneck rather than gene flow explains the high genetic diversity found in all studied populations. However, we cannot exclude the possibility that precocious parr also contribute to genetic variability. Concerning the sensitivity of methods, the results raise the importance of number and variability of neutral markers, of the Bayesian priors and of the structure of populations. Large populations require higher numbers of markers for long-term coalescence estimators and larger sample sizes for short-term coalescence estimators. Even in the case of low migration, establishing a conservation program should rely on both short-term and long-term estimates of population size.
6. Perspectives
Because variability of effective population size is a main factor determining the risk of extinction (Waples, Reference Waples2002), fluctuating population size is an important consideration for evolution and conservation. Some models have been recently developed (Drummond et al., Reference Drummond, Rambaut, Shapiro and Pybus2005; Heled & Drummond, Reference Heled and Drummond2008) to calculate the fluctuation of current Ne from the most recent ancestor. Although restricted for the moment to sequence data (mitochondrial or viral DNA) and large intervals of time, extension to microsatellite markers could help provide information on more recent history.
The authors wish to thank N. Jeannot, F. Marchand, G. Evanno and D. Azam (INRA Rennes, France); E. Prévost (INRA St Pée sur Nivelle, France); S. Burns, J. Woods, J. Reid, D. Ferguson and S. Grant (Spey Fishery Board Research Office, Scotland); D. Knox and E. Verspoor (Fisheries Research Services Perthshire, U.K) and J. Stoddart (Kyle of Sutherland District Salmon Fishery Board, U.K) for providing support collection of scale and fin samples. Genotyping were done at the Toulouse Genopole Platform (http://www.genotoul.fr/) with the help of K. Fève and J. Riquet that we fully thank. This work was financial supported by INRA, Laboratoire de Génétique Cellulaire (UMR 444), INRA-ENVT, BP 52627, 31326 Castanet Tolosan Cedex, France.
