Abstract

DNA polymorphism at 22 loci was studied in an average of 47 Norway spruce [Picea abies (L.) Karst.] haplotypes sampled in seven populations representative of the natural range. The overall nucleotide variation was limited, being lower than that observed in most plant species so far studied. Linkage disequilibrium was also restricted and did not extend beyond a few hundred base pairs. All populations, with the exception of the Romanian population, could be divided into two main domains, a Baltico–Nordic and an Alpine one. Mean Tajima's D and Fay and Wu's H across loci were both negative, indicating the presence of an excess of both rare and high-frequency-derived variants compared to the expected frequency spectrum in a standard neutral model. Multilocus neutrality tests based on D and H led to the rejection of the standard neutral model and exponential growth in the whole population as well as in the two main domains. On the other hand, in all three cases the data are compatible with a severe bottleneck occurring some hundreds of thousands of years ago. Hence, demographic departures from equilibrium expectations and population structure will have to be accounted for when detecting selection at candidate genes and in association mapping studies, respectively.

LEVEL of nucleotide polymorphism, extent and pattern of linkage disequilibrium (LD), and degree of population differentiation are fundamental population genetics parameters that are strongly influenced by evolutionary forces that acted in the past. Their analysis can therefore be used to infer past demographic history and selection events. Solid reconstructions of past demographic events based on a large number of loci are needed to detect genomic areas that are under selection since, if the population departs from the standard neutral model, current neutrality tests that compare the observed polymorphism pattern to that expected under the standard neutral model cannot be used (see, for example, Thornton and Andolfatto 2005). In a few intensively studied species, the availability of extensive genomic data and powerful coalescent-based estimation methods are enabling such reconstructions, thereby greatly facilitating the detection of loci under selection in genome scans (e.g., Akey et al. 2002; Schaffner et al. 2005; Wright et al. 2005). In other organisms, while such fine-tuned reconstructions are still out of reach, more limited surveys of nucleotide variation, coupled to coalescent simulations still do allow the evaluation of different demographic models. For example, Haddrill et al. (2005) used multilocus neutrality tests, measures of linkage disequilibrium, and coalescent simulations to show that simple bottleneck models were sufficient to account for most, if not all, polymorphism features of Drosophila melanogaster. Such approaches have not yet been applied to conifer species, although they may be the key to the understanding of some of the intriguing patterns of nucleotide polymorphism that have emerged from initial surveys. Estimates of nucleotide diversity reported so far in conifers have been much lower than expected on the basis of their life-history traits and the high heterozygosity levels observed at allozyme loci for these species (Hamrick and Godt 1996). The average πsilent was 0.0064 in Pinus taeda (Brown et al. 2004) and ∼0.0041 in P. sylvestris (Dvornyk et al. 2002; García-Gil et al. 2003). In Norway spruce, nucleotide diversity seems also low (πs = 0.0041 for 21 EST loci sequenced across 12 individuals; S. degli Ivanissevich and M. Morgante, unpublished data). In P. taeda, Brown et al. (2004) concluded that the low nucleotide diversity could be the result of a particularly low mutation rate (on the order of 1.7 × 10−10/bp/year, i.e., an order of magnitude lower than in angiosperms) combined with a low effective population size (5.6 × 105) due to population fluctuations during the late Pleistocene and the Holocene. This low effective population size was derived from the relationship θ = 4Neμ, using the mutation rate per generation, and hence a standard neutral model was assumed. An alternative explanation of the low nucleotide diversity could be the presence of repeated selective sweeps but this seems unlikely in conifers as current estimates suggest that LD does not extend beyond a few hundred or thousand base pairs (Neale and Savolainen 2004). However, LD is known to vary extensively along the genome, and at different scales (e.g., Myers et al. 2005), and current estimates in conifers are based only on a handful of loci in a few species, so this picture might change drastically as data accumulate.

In this study, we surveyed DNA polymorphism at 22 loci in an average of 47 haplotypes from seven populations representative of the Norway spruce natural range. Eleven loci were candidate genes for seasonal growth cessation and the remaining ones were randomly chosen from an EST database. The latter are a priori not related to seasonal growth cessation, a trait showing strong clinal variation (Ekberg et al. 1979). The aim of this study was to assess nucleotide diversity, population structure, and LD and address the following questions:

  1. Are nucleotide polymorphism and LD in spruce as limited as in other conifer species and do the patterns indicate departure from the standard neutral model?

  2. Do some of the candidate genes depart from the average pattern?

  3. Do nucleotide polymorphisms display patterns of population structure similar to allozymes and cytoplasmic markers? Those markers distinguished two main domains, one covering northeastern Russia and Scandinavia (Baltico–Nordic domain) and the other centered on the Alps and extending into Poland (Alpine and Central European domain, hereafter called the Alpine domain) (Lagercrantz and Ryman 1990; Bucci and Vendramin 2000; Vendramin et al. 2000; Sperisen et al. 2001). The domains mirror the natural distribution of the species into two main geographical areas with smaller, more isolated pockets in the Carpathians and the Balkans.

  4. Is it indeed so that the Alpine domain has a lower level of diversity than the Baltico–Nordic domain and went through a bottleneck as suggested by Lagercrantz and Ryman (1990) while the Baltico–Nordic domain is closer to an equilibrium neutral model?

  5. If populations went through a bottleneck, what were its characteristics (time of occurrence and overall severity) and could a bottleneck help explain the particularly low level of nucleotide polymorphism?

To address these questions various population growth and bottleneck models were evaluated through coalescent simulations.

MATERIALS AND METHODS

Plant material:

Picea abies seeds were collected from nonadjacent maternal trees in seven natural populations or artificial populations representing local gene resources (Table 1). Seedlots were partitioned between the Uppsala and the Udine laboratories. Seeds were soaked in water overnight and haploid DNA was extracted from megagametophytes using a CTAB method (Doyle and Doyle 1990). In each population, both laboratories mostly used megagametophytes from the same individuals for sequencing but in a few cases additional megagametophytes from different individuals were included.

View this table:
TABLE 1

Geographical coordinates of the Picea abies populations analyzed in this study

Sequencing:

To identify functional candidate genes, we performed BLAST, BLASTX, and TBLASTX searches (Altschul et al. 1997) in the NCBI and the loblolly pine EST (http://pinetree.ccgb.umn.edu/) databases, using published sequences of genes from the photoperiod and vernalization pathways in model organisms, mainly Arabidopsis thaliana (Simpson and Dean 2002; Yanovsky and Kay 2003; Hayama and Coupland 2004). A total of 11 growth cessation candidate genes were chosen in P. abies (Table 2), showing similarity with the A. thaliana genes co (constans, Putterill et al. 1997), cry1 (cryptochrome1, Lin et al. 1996), ebs (early bolting in short days, Piñero et al. 2003), gi (gigantea, Fowler et al. 1999), pat1 (phytochrome A signal transduction1, Bolle et al. 2000), phyA and phyB (phytochrome A and phytochrome B, Sharrock and Quail 1989), and vip3 (vernalization independence 3, Zhang et al. 2003). Consistent with the three phytochrome gene lineages reported in gymnosperms (phyN, phyO, and phyP; Schmidt and Schneider-Poetsch 2002), we identified multiple gene copies and pseudogenes within the P. abies phytochrome gene family by cloning (data not shown). Therefore, nonoverlapping regions of phyN and phyP genes were treated as different loci. PCR primers were designed from loblolly pine or Scots pine EST sequences or from P. abies specific sequences, obtained through RT–PCR and RACE reactions, using the Primer3 software (Rozen and Skaletsky 2000). Control loci a priori not involved in the photoperiod or vernalization pathways (se121, se129, se1100; se1151, se1358, se1364, se1368, se1390, se1391, xy225; xy1420) were selected from a pilot resequencing survey of 21 EST-based loci across 12 P. abies individuals (S. degli Ivanissevich and M. Morgante, unpublished data). Selection criteria included ease of amplification and sequencing with the amplification primers and the presence of at least two polymorphic sites detected in the pilot survey. All genes were amplified from haploid megagametophyte DNA with the Phusion DNA Polymerase (Finnzymes, Espoo, Finland) or AmpliTaq Gold DNA Polymerase (Applied Biosystems, Foster City, CA) and directly sequenced from the PCR product either with BigDye v3.1 and run on a ABI 3730 (Applied Biosystems) or with Dyenamic ET terminators and run on a MegaBace 1000 (GE Healthcare, Piscataway, NJ). Most gene regions were covered by two or more reads. Sequences were base called and assembled with PHRED and PHRAP (Ewing and Green 1998; Ewing et al. 1998) and visualized and edited with CONSED (Gordon et al. 1998). A putative SNP was considered true when PHRED quality scores of the different variants exceeded 25.

View this table:
TABLE 2

Nucleotide variation, haplotypic diversity, and neutrality tests in 22 Picea abies loci sequenced across seven populations

Nucleotide diversity analysis:

Estimates of standard population genetics parameters and neutrality test statistics were calculated for each locus with the DnaSP v. 4.0 software (Rozas et al. 2003). Insertions or deletions are reported, but were excluded from further analyses. The level of polymorphism for each locus was estimated as both haplotype and nucleotide diversities.

Population structure:

Population differentiation was first estimated with Wright's fixation index FST (Wright 1951). F-statistics of each gene were computed from allele frequencies as variance component ratios with the locus-by-locus AMOVA approach (Excoffier et al. 1992) implemented in the Arlequin software (Schneider et al. 2000). Single-gene FST's and overall FST were obtained by summing variance components over nucleotide loci Math according to Weir and Cockerham (1984). The significance of FST was tested by comparing the observed value with the distribution of FST after 10,000 permutations of sequences among populations.

The genetic structure of the Norway spruce sample was also investigated with the model-based clustering algorithm implemented in STRUCTURE v. 2.1 (Pritchard et al. 2000; Falush et al. 2003). We used the admixture model on a subset of the data represented by 105 parsimony informative unlinked loci, that is, loci between which Fisher's exact test with Bonferroni correction (see Linkage disequilibrium) was not significant. Ten runs with a burn-in of 100,000 and a run length of 500,000 iterations were performed for a number of clusters from K = 1 to K = 7, allowing for correlation of allele frequencies between clusters. As individuals are assigned to clusters to achieve linkage and Hardy–Weinberg equilibria, the fact that different loci were sequenced from different megagametophytes from the same individual should not affect the results. When, as occurred in a few cases, loci were sequenced in gametophytes from different mother trees, sequences were not pooled to create haplotypes. Instead other loci were coded as missing values. In any case, it should be emphasized that we use Structure here primarily as a tool to explore the data rather than a first step in an association study, which would clearly require a stricter control of the population structure. More generally, and as pointed out by Setakis et al. (2006), the notion of subpopulation is a theoretical construct that will only imperfectly reflect reality and therefore the resulting clusters should not be interpreted too literally.

Linkage disequilibrium:

The level of linkage disequilibrium between parsimony-informative sites within genes was estimated as r2, the mean squared correlation in allelic state between pairs of SNPs, using DnaSP. Significance of the associations between SNPs was determined with Fisher's exact test with Bonferroni correction. The overall decay of LD with physical distance within genes was evaluated by nonlinear regression of r2 on distance between sites in base pairs (Remington et al. 2001). We used the Hill and Weir (1988) expectation of r2 between adjacent sites,Mathwhere C is the population recombination parameter (ρ = 4Ner) and n the sample size, and replaced C by C × distance in base pairs when fitting the formula to our data using the nonlinear regression (nls) function in the R software (R Development Core Team 2005).

Statistical test of neutrality and evaluation of alternative models:

To test for departure from the standard neutral model the mean value of Tajima's D and Fay and Wu's H over loci was compared with the distribution of mean values from coalescent simulations using code kindly provided by P. Andolfatto and described in Haddrill et al. (2005). Both test statistics compare two estimates of θ: Tajima's D measures the standardized difference between π and θW (Tajima 1989) while Fay and Wu's H (Fay and Wu 2000) measures the difference between π and θH. The former is most sensitive to an excess of rare variants whereas the latter is most sensitive to an excess of high-frequency-derived variants. Both D and H are expected to be close to zero under the standard neutral model. All tests were carried out with recombination, as the lack of recombination makes the tests overly conservative. An estimate of the population recombination parameter ρ = 4Ner was obtained with the composite-likelihood method of Hudson (2001) adapted to finite-sites models as implemented in the software LDHat v. 2.0 (McVean et al. 2002), for col1 (this study) and for two other genes, ft1 and toc1 (our unpublished data). Estimates of ρ were 8.5 × 10−3, 4.7 × 10−3, and 2.6 × 10−3/bp, respectively. For col1 we used the estimated ρ-value for that gene, while for all other genes we used the average of the three estimates, 5.3 × 10−3. The ancestral state of nucleotides, required for Fay and Wu's H-test, was inferred by using a single sequence of P. mariana, P. glauca, P. sitchensis, or, in the case of col1, P. sylvestris as an outgroup.

Coalescent simulations were also used to evaluate two types of alternative models: exponential growth models and bottleneck models followed by exponential growth. Briefly, it is now well established that the ranges of tree taxa went through cycles of contraction and expansion in response to climate changes during the late Quaternary (Bennett 1997). In Norway spruce, as in most species, however, the severity of the contractions, the size and location of the refugia, and the rate of the ensuing growth are still poorly characterized and, consequently, we modeled bottlenecks of various severity and ages and considered models with different growth rates. We also assessed the effects of repeated bottlenecks (data not shown). Importantly, because all times are in units of effective population size for which we do not have any independent estimate, the age and severity of the bottleneck cannot be defined exactly. So our primary aim in this study was to test whether the data could be better explained by a bottleneck than by the standard neutral model in the first place rather than to obtain a fine characterization of that bottleneck. The difficulty in characterizing a bottleneck is compounded by the fact that the effect of a bottleneck on the frequency distribution of mutations segregating in a population depends on the time at which the bottleneck ended and its strength, which is approximately a function of the ratio of its severity (the magnitude of the reduction in population size) to its duration. Hence different combinations of the three parameters can lead to the same nucleotide frequency spectrum (Fay and Wu 1999; Voight et al. 2005; Wright et al. 2005). The bottleneck models were tested over a grid of parameter values: the severity varied between 0.0004 and 0.001 and the time at which the bottleneck ended (t_end) varied between 0.001 and 0.0095. The length of the bottleneck was fixed to 0.0015 in all bottleneck models. Time measures are in units of 4N0 generations from the present and the severity of the bottleneck is in units of the current population size. The coalescent simulations were run with recombination estimated as above. Because we would actually need to run the simulation using the recombination rate in the ancestral population, for which we have no estimates, to assess the robustness of the results, we also ran a subset of demographic scenarios with twice that value and without recombination. Details on the methods can be found in supplemental material at http://www.genetics.org/supplemental/.

RESULTS

Nucleotide variation:

Sequence variation was obtained for all 22 loci (supplemental Table 1 at http://www.genetics.org/supplemental/) in an average of 47 megagametophytes, ∼7 from each of seven P. abies populations. A total of 16,161 bp were aligned over the 22 genes, of which more than half was coding sequence, resulting in a total of ∼760 kb of sequence information across individuals. Insertions/deletions (indels) covered 130 bp. They comprised seven microsatellites with a motif length of 1, 3, 4, or 9 bp in, respectively, ebs, phyo, and se1390 (1 bp); se1364 and xy225 (3 bp); ebs (4 bp); and se1100 (9 bp). The microsatellites were located in noncoding regions, except in the case of se1364, where a 3-bp repeat in the coding region produced no shift in the reading frame. The remaining indels were located in noncoding regions, namely an 11- and a 54-bp stretch in col1, a 27-bp stretch in se1100, and an 8-bp stretch in each pat1 and se1368. Indels were excluded from further analyses.

We identified a total of 230 segregating sites, of which 89 were singletons and 141 were parsimony-informative sites (Table 2). This corresponds to 1 SNP every 69 bp. One parsimony-informative site with three variants was found in se1420; it was excluded from further analyses. Forty (17.4%) SNPs were amino acid replacement substitutions. Statistics of sequence variation are summarized in Table 2. Total nucleotide diversity πt was between 0.0002 and 0.0068 (average πt = 0.0021) and silent nucleotide diversity, including synonymous and noncoding positions, varied between 0 and 0.0176 (average πs = 0.0039; loci for which no similarity was found in Blast searches were considered only in the calculation of the total nucleotide variation). Nonsynonymous nucleotide diversity was on average 4.2 times smaller than silent nucleotide diversity and varied from 0 to 0.0029 (average πa = 0.0009). Control loci were more polymorphic than candidate loci [πt(controls) = 0.0027 ± 0.0019 (SD) vs. πt(candidates) = 0.0014 ± 0.0008 (SD)]. It is difficult to speculate on the cause of this difference since (i) it might result simply from differences in sampling approach for the two groups of genes and (ii) different species were used as outgroups for the different genes, making estimates of the average mutation rate in the two groups complicated. With these caveats in mind, considering only the genes for which the same outgroup, P. taeda, was used, the average divergence at silent sites among control loci was larger than the average divergence among the candidate loci but the standard deviations were very large [0.1777, SD = 0.115 (n = 6) vs. 0.077, SD = 0.065 (n = 6)].

Population structure:

FST-values varied between 0 and 0.289 among loci and revealed substantial differentiation among the seven populations, with an overall value of FST = 0.117 over all 229 SNPs (supplemental Table 2 at http://www.genetics.org/supplemental/). Romania was the most differentiated population in pairwise comparisons (0.194 ≤ FST ≤ 0.262, Table 3).The two-level variance partitioning revealed a high differentiation (FST = 0.147) between the Baltico–Nordic domain (Northern Sweden, Southern Sweden, Russia), the Alpine domain (Italy, Switzerland, Germany), and the Carpathian domain (Romania) (data not shown). Populations within the Baltico–Nordic domain were significantly, though weakly differentiated (FST = 0.025, P ≤ 0.05), whereas populations from the Alpine domain were not significantly differentiated (FST = 0.015, data not shown).

View this table:
TABLE 3

Pairwise measures of population differentiation

The STRUCTURE program revealed the highest likelihood for K = 4 clusters (average log probability of data Ln P(D) = −1452.01 ± 10.39, SD); however, biologically meaningful genetic structure was already detected at K = 3 with Ln P(D) = −1518.26 ± 34.16 (SD). With K = 3, the Baltico–Nordic, the Alpine, and the Carpathian domains were essentially distinguished (Figure 1). The Romanian population was the most distinct with 92.7 ± 0.32% (SD) of ancestry in cluster 1. All other populations were fairly admixed. Populations from the Baltico–Nordic domain had their largest proportion of ancestry, 59.0 ± 0.69% (SD) in cluster 2 while populations from the Alpine domain had theirs in cluster 3 (68.0 ± 0.66%). The populations from southern Sweden (41.0% in cluster 2, 36.8% in cluster 3) and Italy (33.9% in cluster 1, 42.8% in cluster 2) were even more admixed. With K = 4, the structure of K = 3 was confirmed and a fourth cluster accounted for 20–27% of the ancestry of populations from southern Sweden, Germany, Switzerland, and Italy.

Figure 1.—

Structure analysis of the seven populations when K = 3 clusters are assumed. The Baltico–Nordic domain includes southern Sweden, northern Sweden, and Russia and the Alpine domain includes Switzerland, Germany, and Italy.

The results on among-population differentiation suggest different evolutionary histories for the Baltico–Nordic vs. the Alpine part of the Norway spruce range and a particular situation for Romania. Diversity estimates were lowest in Romania, with πT = 0.0012 lower than πT ≥ 0.0016 in other populations (unilateral paired t-tests: P ≤ 0.05 except for Germany where P = 0.053 and Switzerland where P = 0.091). Mean genetic diversities in the Baltico–Nordic and in the Alpine domains were very close (πT = 0.0022 vs. πT = 0.0017).

Linkage disequilibrium:

A low level of linkage disequilibrium was observed within genes, with an average r2 = 0.115 and 75 significant exact tests after Bonferroni correction among 1411 pairwise comparisons between informative SNPs. LD decayed fast within genes, with r2 dropping below 0.2 within ∼100 bp (Figure 2).

Figure 2.—

Plot of the squared correlation of allele frequencies (r2) vs. distance in base pairs between polymorphic sites across 22 loci for different subsets of populations.

Statistical tests of neutrality and evaluation of alternative models:

Mean values of both Tajima's D- and Fay and Wu's H-statistics were negative, with values of −0.88 and −0.74, respectively (excluding ebs for which no outgroup was available; Table 4). Coalescent simulations using both statistics led to the rejection of the standard neutral model and population growth, but simulations that assumed a severe and rather ancient bottleneck followed by moderate population growth were consistent with the data. Table 4 gives the values of mean Tajima's D and mean Fay and Wu's H for a subset of the models that were tested. Tajima's D obtained from simulations under the standard neutral model is significantly larger than the observed value, whereas under the growth model Tajima's D no longer differs from the observed value but Fay and Wu's H is now significantly larger than the observed value. Various growth models were tested (supplemental Table 3 at http://www.genetics.org/supplemental/): none led to negative values for both mean D and mean H. Unless it was extremely severe a recent bottleneck was inconsistent with the data as it led to an excess of common variants and a positive Tajima's D or required an extremely large ancestral effective population size (θ as large as 16; Figure 3). On the other hand, as long as the bottleneck is ancient enough, the data can be explained by different combinations of time at which the bottleneck ended and severity without requiring unrealistically large θ-values (Figure 3). The same analysis was also carried out within the Baltico–Nordic and Alpine domains and led to similar conclusions, the acceptance regions being somewhat larger in the Baltico–Nordic domain than in the Alpine domain (supplemental Figures 1 and 2 at http://www.genetics.org/supplemental/).

Figure 3.—

Evaluation of different bottleneck models in the total data. (Top) Significance level of the multilocus neutrality test for different combinations of severity and time at which a bottleneck ended (t_end). The duration of the bottleneck was 0.0015. The P-value reported was in all cases that for Tajima's D. The P-value for Fay and Wu's H was always >0.05. From darker to lighter shading: P > 0.05, 0.01 < P ≤ 0.05, 0.001 < P ≤ 0.01, and P ≤ 0.001. (Bottom) Corresponding average θ-values used in coalescent simulations. Lightest shading, θ > 16; darkest shading, θ < 5.

View this table:
TABLE 4

An evaluation of alternative demographic models for the total population and the Baltico–Nordic and the Alpine domains

Only five genes showed significant Tajima's D-values, namely col1, ebs, phynrII, vip3, and se121 when Tajima's D and Fay and Wu's H were calculated for individual loci (Table 2). To assess whether demography alone could explain those departures or whether the frequency spectrum at those loci would still depart from the rest of the genome when demography is taken into account we tested them against a bottleneck model that could not be rejected globally. Ebs was not considered as we did not have an outgroup and phynrII was discarded because it had only two segregating sites. The bottleneck model was accepted for col1 and se121 but was rejected for vip3 (Table 5). Additional factors, such as selection or a more complex demographic model, might therefore need to be invoked to account for the polymorphism at vip3. However, we note that no significant signal of selection could be detected on any SNP using Beaumont and Balding's (2004) method (data not shown).

View this table:
TABLE 5

Evaluation of a bottleneck model at individual genes that depart from the standard neutral model

DISCUSSION

Norway spruce has a low to moderate level of nucleotide diversity (πs = 0.0039, θWs = 0.0058), low levels of LD, which decayed by 50% within <100 bp, and a moderate level of population structure (FST = 0.12). Using multilocus tests based on summary statistics of the allele frequency spectrum we showed that the standard neutral model can be rejected and that a severe bottleneck predating the Last Glacial Maximum is sufficient to explain the data. This is true when all populations are considered but also within both the Baltico–Nordic and Alpine domains when those are analyzed separately. Hence, although nucleotide diversity was slightly higher in the Baltico–Nordic than in the Alpine domain the two domains seem to have experienced rather similar demographic histories.

Nucleotide diversity:

The average level of silent nucleotide diversity in P. abies, πs = 0.0039, confirmed earlier results from S. degli Ivanissevich and M. Morgante (unpublished data), who found πs = 0.0041 for 21 EST loci sequenced across 12 individuals (note that the 11 control loci we analyzed were selected on the basis of that study), and supports the contention that conifers are characterized by a low level of nucleotide diversity. Compared to other conifers, the level of silent polymorphism was of the same order of magnitude as that in Cryptomeria japonicas = 0.0038 across 7 genes, Kado et al. 2003) and P. sylvestriss ≈ 0.0041 across 14 genes, Dvornyk et al. 2002; García-Gil et al. 2003), but lower than that in P. taedas = 0.0064 across 19 wood-production candidate genes, Brown et al. 2004; πs = 0.0079 across 18 drought stress candidate genes, González-Martínez et al. 2006). These estimates of silent nucleotide diversity are higher than that in soybean (πs = 0.0015, Zhu et al. 2003) but twofold lower than that in A. thalianas = 0.0083, Schmid et al. 2005) and an order of magnitude lower than estimates in aspen (πs = 0.0160, Ingvarsson 2005) and wild relatives of maize (πs = 0.012–0.013, Tiffin and Gaut 2001). Hence, our results indicate that nucleotide diversity in the Norway spruce gene pool is indeed low.

Variation in average nucleotide diversity estimates across species can be caused by a combination of factors such as differences in individual sampling strategies, parts of the genome considered, selection, demographic history, and differences in mutation rate. The studies cited above were based on samples covering the entire species distribution ranges, or wild and cultivated genotypes of different origins, so differences in individual sampling are unlikely to account for the magnitude of the variation in estimates among species. Silent variation varied 30-fold across genes in our study and 50-fold in P. taeda (Brown et al. 2004), so the studied genes do not appear to be biased toward a particular group. Because the EST genes were selected to be variable, our estimate might even be biased upward. Estimates of average polymorphism could also be reduced by selective sweeps that diminish variation at and around particular genes or by purifying selection. However, the limited amount of LD makes selective sweeps an unlikely explanation for low nucleotide diversity and recurrent hitchhiking events would not account for the negative values of Fay and Wu's H (Przeworski 2002; Haddrill et al. 2005). Similarly, models of weak negative selection predict an excess of low-frequency-derived mutations and hence a positive value of Fay and Wu's H, which is not consistent with our data. In brief, while selection may partly explain the low level of nucleotide variation at individual loci it does not seem to be sufficient to explain the low level across loci.

Two possible explanations therefore remain for the relatively low level of polymorphism in conifers, namely demographic history and/or mutation rates. On the basis of sequence variation at 19 nuclear genes (amounting to a total of almost 18,000 bp) Brown et al. (2004) estimated the substitution rate per year to be 1.17 × 10−10 in P. taeda, a value similar to that reported for P. sylvestris (Dvornyk et al. 2002) and an order of magnitude lower than angiosperm mutation rates. This estimate was based on a divergence time between P. pinaster and P. taeda of 120 million years. There are grounds, however, to question this mutation rate estimate. First, the divergence time retained to calculate it corresponds to the early diversification of the Pinaceae in the early Cretaceous (∼120–140 MYA), not to the divergence time of two species within the genus. Wang et al. (2000) inferred that Pinus species diverged from one another in the early to mid-Cretaceous (∼70 MYA), which is consistent with the first appearance of Pinus in the fossil record in the early Cretaceous. Picea species appeared later on, in the middle Pliocene (∼45 MYA), and apparently diversified ∼20 MYA, if not later (Bouillé and Bousquet 2005). Hence 120 MYA is likely to be a gross overestimate of the divergence between P. pinaster and P. taeda and 1.17 × 10−10 an underestimate of the mutation rate. Bouillé and Bousquet (2005), considering three nuclear genes (amounting to a total of ∼2000 bp), obtained a mutation rate of 2.23 × 10−10 to 3.32 × 10−10 in Northern American Picea species. As nucleotide diversity varies a lot across loci this estimate cannot be taken at face value but is, in any case, two- to threefold higher than the one reported by Brown et al. (2004). Second, estimates of molecular divergence between Pinus and Picea based on an extensive EST database lead to an estimate of the mutation rate of ∼1 × 10−9/year if the divergence time between pine and spruce is that of the diversification of the Pinaceae in the early Cretaceous (∼120–140 MYA) (Savolainen and Wright 2004). Finally, Willyard et al. (2006) used divergence at multiple nuclear and chloroplast loci, exemplar taxa, and two calibration points to show that divergence times among pine lineages have often been overestimated and, consequently, absolute mutation rates have been underestimated. They obtain a nuclear silent mutation rate in Pinus of 0.70–1.31 × 10−9 sites/year. Hence, the particularly low levels of nucleotide diversity in Norway spruce are probably not exclusively due to low mutation rates and we may have to turn to population demographic history for additional explanations.

Population history:

The main outline of Norway spruce population history inferred from this DNA polymorphism survey is the following. As shown previously with other molecular markers [allozymes (Lagercrantz and Ryman 1990), AFLP (Acheré et al. 2005), and cytoplasmic DNA (Vendramin et al. 2000; Sperisen et al. 2001)], the Norway spruce population is today genetically and geographically divided into two main domains, namely the Baltico–Nordic domain and the Alpine Central European domain, and a more limited one, the Carpathian domain, represented in the present survey by a single population, Romania. This population had a very limited polymorphism and may not be representative of the Carpathian domain. The estimate of overall population differentiation that we obtained using FST (0.117) is substantially higher than that previously reported by Lagercrantz and Ryman (1990) using isozymes (0.052) on a larger set of 70 populations covering a similar geographic range and by Acheré et al. (2005) using AFLPs and SSRs (0.02) on a more limited set of populations. This could, in part, be due to differences in sampling and differences in levels of variation among different types of markers (Charlesworth 1998). The choice of candidate genes putatively involved in controlling phenology-related traits for which ample variation exists among these populations (QST = 0.729, R. Liesch and M. Lascoux, unpublished data) could also explain this difference, if they were under selection. However, no significant difference between candidates and control genes is visible in population differentiation levels as estimated by FST (data not shown).

The split among the two main geographic domains has been dated to a maximum of 40,000 years (Lagercrantz and Ryman 1990), coinciding with the time estimated from pollen analysis. Previous to that, our analysis suggests that the whole population went through a rather severe bottleneck. We did not attempt to date the bottleneck precisely but recent bottlenecks failed to generate negative values for both H and D and too severe ones would require unrealistically high values of θ in the ancestral population to explain the data. A more recent bottleneck may also be incompatible with the low level of linkage disequilibrium. Because of the fairly large set of bottleneck parameters that are compatible with the data it is difficult to associate the bottleneck with a particular climatic event. However, climate reconstructions extending back 400,000 years (e.g., Petit et al. 1999) show that the average temperature fluctuated with an amplitude of ∼10° and a periodicity of ∼100,000 years. The bottleneck(s) suggested by the genetic data could then correspond to one of the abrupt changes in temperature that took place during the quaternary. More complex demographic models, metapopulation models, or glacial cycles models, for instance (Wakeley and Aliacar 2001; Jesus et al. 2006), may provide an even better fit to our data, but would be more difficult to justify and model at that stage. Finally, our inference was based on both coding and noncoding DNA. It would certainly have been better to use only noncoding DNA as was done by Haddrill et al. (2005). However, as there was no strong evidence of selection at loci considered individually, and given the low level of linkage disequilibrium ruling out strong hitchhiking effects, we feel that this may not have altered our conclusion. In summary, we therefore conclude that, even if demography alone is unlikely to explain the low nucleotide variation in all coniferous species, it provides a simple explanation, at least in Norway spruce.

Linkage disequilbrium:

Linkage disequilibrium was limited within genic regions; LD decayed by half within <100 bp, confirming earlier results of Rafalski and Morgante (2004). However, the present analysis of LD is biased by the pattern in col1, the only gene for which a fragment >3000 bp was sequenced, and indeed estimates of LD at two other long fragments (our unpublished data) were somewhat higher. The pattern of LD was only weakly influenced by population structure, since similar results were obtained when sequences from Romania, the most divergent population, were not included. The rapid decay of LD, consistent with the prevailing outcrossing mating system and the high level of heterozygosity of this species, was similar to the one observed in another outcrossing plant, maize (Tenaillon et al. 2001). The very low level of LD could also provide an explanation for the differences in variability at allozyme and nucleotide levels: a limited number of segregating sites per locus recombining freely can lead to a high haplotype diversity.

Conclusions:

The level of population structure detected in this study and the overall departure from the standard neutral model of spruce populations imply that these factors will have to be taken into account when carrying out association-mapping studies (Marchini et al. 2004; Campbell et al. 2005; Helgason et al. 2005) and when interrogating SNP databases for signatures of natural selection (Akey et al. 2002), respectively. The rapid decay of LD in spruce will allow high-resolution mapping in association studies, given that the right candidate genes are chosen, but will also require a high-density marker screening due to the limited predictive power of single SNPs over neighbor sequence diversity. Finally, if confirmed by more extensive studies, the rapid decay of LD also implies that hitchhiking is likely to have played a limited role in the species evolution.

Acknowledgments

We thank Skogforsk, Randolph Schirmer, Felix Gugerli, Magdalena Palada, Vladimir Semerikov, and Giovanni Vendramin for providing Picea abies seeds. M.L. thanks Peter Andolfatto for kindly providing the code used to carry out multilocus tests of neutrality and answering questions about it and Erik Lagercrantz for writing a Ruby script that greatly facilitated these analyses. We are grateful to two anonymous referees for constructive comments on the manuscript. This research was funded by the European Commission, project TREESNIPS (QLRT -2001-01973), the Carl Tryggers Foundation, the Philip-Sörensen Foundation, and the Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning. M.H. acknowledges a “Mobility of Researchers” grant from the National Research Fund of Luxembourg and is currently a postdoctoral researcher at the National Fund for Scientific Research of Belgium.

Footnotes

  • Sequence data from this article have been deposited with the EMBL/GenBank Data Libraries under accession nos. AM267499AM268023 and AM268895AM269360.

  • 1 These authors contributed equally to this study.

  • 2 Present address: Laboratoire d'Eco-Ethologie Evolutive - CP 160/12, Université Libre de Bruxelles. 50, Av. F.D. Roosevelt, B-1050 Bruxelles, Belgium.

  • Communicating editor: M. Nordborg

  • Received August 20, 2006.
  • Accepted October 2, 2006.

References

View Abstract