Introduction 

Norway spruce (Picea abies [L.] Karst.) is among the most important conifer tree species in Europe and in Slovenia. It can grow on diverse sites, and often forms continuous stands, either alone or in mixtures ([38]). In general, it exhibits high levels of genetic variability within populations and low levels of differentiation ([18], [23], [41], [39], [37], [20]).

Still, the central-peripheral population hypothesis (CPH) of Mayr ([19]) and the extension of the abundant center model to population size and gene flow ([8]) predict that peripheral populations (i.e., populations at the range limits) would have reduced genetic diversity and increased genetic differentiation, as compared to core populations (i.e., those from the inner part of the species’ range), due to supposedly lower effective population size, higher genetic drift, founder effect and isolation, or a combination of these factors.

This hypothesis has been formerly investigated in many plant and animal species. As for plants, a reduced within-population genetic diversity towards periphery of the species’ range has been found in 64.2% and increased differentiation in 70.3% of the studies, as reviewed by Eckert et al. ([8]), thus supporting the CPH. Also for forest trees, empirical support of CPH varies across studies. For example, range-wide diversity and differentiation are consistent with CPH for Thuja occidentalis ([25]) and Pinus contorta ([11]), but not for Pinus strobus ([29]) and Liriodendron chinense ([44]). On the other hand, a simulation study considering the effect of the interaction between migration, genetic drift, and time of population establishment on genetic diversity of peripheral vs. core populations, found that the pattern of genetic diversity of passive dispersers such as trees is not expected to conform to CPH ([7]).

Looking at vertical gradients, genetic diversity varies following four distinct patterns, one of them being concordant with CPH ([24]). However, only 11 out of the 42 studies on forest trees reviewed by Ohsawa & Ide ([24]) were consistent with the expectation of highest diversity in the core (i.e., the vertical center of the gradient) and lower diversity at periphery (i.e., the distal ends of the gradient), while 53% of studies exhibit higher or equal genetic diversity at the upper range periphery than in the gradient center (core). This is important as populations at the upper range periphery are expected to colonize higher elevation sites in response to climate change ([9]).

Here, we contribute empirical knowledge to the debate on patterns of genetic diversity and differentiation between core and peripheral populations of Norway spruce. The studied populations are located in a small area at the upper tree line or in Karst sinkholes with vegetation inversion, which includes two distinct cases of vertical periphery. We employed richness and equitability measures based on isoenzyme and nuclear SSR markers to test CPH along an elevation gradient.

Genetic richness is a measure of the total amount of diversity which is more sensitive to the presence of rare alleles, and it is likely reduced by the stochastic processes occurring at the upper range limits (e.g., bottlenecks, fluctuations in population size). Contrastingly, equitability is more sensitive to high-frequency alleles and reflects how the diversity is partitioned among sites ([6], [8]).

Periphery may or may not represent ecological marginality or vice versa ([9]). Different habitats of peripheral Norway spruce populations at the upper tree line in the Alps may have favored specific genetic characteristics ([2], [14], [4]). Therefore, we also explored the correlation of genetic diversity estimates to environmental factors using the isoenzyme dataset.

  Materials and methods 

Study sites and field sampling

Dormant buds of vital trees > 30 m apart and belonging to dominant and codominant crown classes were sampled from 15 distinct populations of Norway spruce in 1999, 2000 and 2005 (Tab. 1). The number of sampled trees in each population ranged from 30 to 70. In no case, all dominant trees in a population were sampled. All sampled populations were naturally regenerated, putatively autochthonous and large enough, harboring at least 100 reproducing trees, for population size not to be a confounding factor. Natural regeneration was present at all study sites. The stand characteristics, ecological conditions and number of trees sampled are presented in Tab. 1. The location of the studied populations is reported in Fig. S1 (Supplementary material).

The sampled populations were classified as: (i) core populations, i.e., located in the core of the distribution range; or (ii) peripheral populations, growing either at the upper tree line or in a vegetation inversion. The latter is a common phenomenon in Dinaric Karst where the landform depresses below the surrounding area and the temperature decreases towards the depression bottom. The temperature inversion determines an inversion of the vegetation belts, whereby beech forests grow on the hilltop, followed by Norway spruce stands, and dwarf pine and other small shrubs at the bottom (see Fig. S2 in Supplementary material).

Isoenzyme analysis

Enzyme extraction from dormant buds was performed according to Rhodes ([31]). Electrophoresis conditions, staining, and genotyping followed Konnert & Maurer ([16]). The laboratory part of the analysis was carried out at the Bayerische Amt für forstliche Saat und Pflanzenzucht in Teisendorf, Germany. Overall, 9 enzyme systems yielding 15 coding gene loci were analyzed (E.C. Code, loci: 2.6.1.1, Aat-A, Aat-B, Aat-C; 1.1.1.42, Idh-A, Idh-B; 1.1.1.37, Mdh- A, Mdh-B, Mdh-C; 1.4.1.2, Gdh-A; 3.4.11.1, Lap-B; 5.3.1.9, Pgi-B; 2.7.5.1, Pgm-A; 1.1.1.25, Skdh-A; 1.1.1.44, 6-Pgdh-B, 6-Pgdh-C).

SSR analysis

Total DNA was isolated from dormant buds using the DNeasy^® plant kit (QIAGEN, Germany), as per the manufacturer’s specifications. All samples were genotyped at 6 highly variable microsatellite loci: SpAGG3 ([26]); EAC1F04 ([36]); PAAC19, PAAC3 ([35]); paGB8 ([3]); and PGL15 ([30]). Polymerase chain reactions (PCRs) were performed as described by Westergren ([42]). The sizing of the PCR products was performed on an ABI PRISM 310^® automatic sequencer with the accompanying software Gene Mapper^® v. 3.7 (Thermo Fisher Scientific Inc., Waltham, MS, USA).

Data analysis

Genetic diversity analysis was carried out separately for the isoenzyme and SSR datasets. Estimates of genetic diversity (allelic richness (A_R), expected (H_E) and observed (H_O) heterozygosity, and inbreeding coefficient (F_IS) were calculated using the software package SpaGeDi v. 1.4 ([13]). To test for significant difference in mean values of genetic diversities between peripheral and core populations, a t-test with Welch’s modification for unequal variances between groups was calculated in R ([28]). Departure of F_IS from zero was tested using bootstrapping (10 000 permutations) in SpaGeDi. Deviations from the Hardy-Weinberg equilibrium and linkage disequilibrium between loci were tested running 10 000 permutations with the software Genepop v. 4.0 ([32]). Presence of null alleles was checked using the EM algorithm implemented in FreeNA ([5]) followed by correction of the dataset to meet the estimated frequencies of each allele in each population according to the EM algorithm. For the selectively neutral SSRs, the number of migrants (Nm) was estimated based on the private allele method implemented in Genepop.

Pairwise and global F_ST values were calculated and their significance tested using 10 000 permutations with SpaGeDi. Genetic differences between populations were also investigated using a model-based clustering algorithm implemented in Structure v. 2.3.4 ([27], [10], [15]). To estimate the best number of distinct clusters and to average the results of the replicated runs, the CLUMPAK software ([17]) was used. The default model parameters using populations’ priors were used for simulations, allowing number of populations K to vary from 1 to 7 for both datasets. Each run, replicated 8 times, consisted of 400 000 burn-in iterations and 600 000 data collection iterations.

Bonferroni’s corrections were applied to adjust critical values in case of multiple comparisons.

Correlations of genetic diversity estimates to environmental factors

In Alpine conifers, precipitation, temperature and soil play a significant role in shaping genetic variation and local adaptation ([21], [22]). To understand whether genetic diversity (A_R, H_E, H_O,) and F_IS correlate with environmental factors and elevation, Spearman’s correlation coefficient (ρ) was used. Cambium activity in temperate regions is regulated by temperature, rainfall and photoperiod, and reflects the environmental adaptivity of trees ([1]). Radial growth, a result of cambial activity, of Norway spruce adult trees in the Alps is affected positively by high precipitation during May-June period, while higher mean temperature under the same period had a negative effect ([34]). Therefore, we focused our environmental association analyses on mean May to June precipitation and temperature sums. Additionally, a t-test with Welch’s modification was used to test for significant difference in mean values of genetic estimates between groups denoted by nominal variables: bedrock (calcareous and igneous), as a determinant of soil formation, and leaf flushing (21-31 May, 1-10 June). Environmental and flushing data were obtained from the Slovenian Environment Agency data layers; precipitation and temperature for the period 1981-2010 and phenology averages from the 1971-2000 data series.

  Results and discussion 

Quality of the selected markers

Null alleles with frequency ranging between 0.05 and 0.36 (average = 0.21) were detected in all populations at the SSR locus EAC1F04. The average frequency of null alleles at the other SSR loci was much smaller. Still, seven population/locus combinations had frequency of null alleles higher than 0.15, and 11 higher than 0.10 (but lower then 0.15), accounting for half of population/locus combinations. Therefore, the whole dataset was corrected for null alleles to fit allele frequencies estimated by the EM algorithm implemented in FreeNA. All values reported for SSRs are based on this corrected dataset.

The frequencies of null alleles for isoenzymes varied between zero and 0.14, on average 0.01 per population/locus, exceeding 0.10 in five out of 143 population/locus combinations. Therefore, they were not corrected to fit allele frequencies estimated by the EM algorithm.

For none of the pairwise loci combinations within each population the null hypothesis of independent assortment of genotypes could be rejected.

Allele frequency distribution

Isoenzymes

Out of 15 examined loci, 13 were polymorphic in at least one of the examined populations; the two non-polymorphic loci (Mdh-A, Aat-A) were excluded from further analysis. Two to seven alleles were detected for each locus, with a total of 49 alleles across all populations and loci. Twenty-four of these were rare (frequency < 0.05) or absent and 24 common (frequency ≥ 0.05) across all populations. Five alleles were rare in some and common in other populations, but there was no pattern in regards to peripheral and core populations. However, nine of the rare alleles (frequencies between 0.008 and 0.014) were found either in only peripheral (six alleles) or core (three alleles) populations, in five cases in more than two populations (Fig. 1). On the other hand, only two private alleles were found in peripheral populations (both in populations with temperature inversion) and three in core populations.

Fig. 1 - Distribution of allele frequencies (y-axis) at 5 isoenzyme loci exhibiting rare alleles restricted to peripheral or core populations (marked with an arrow). Each of the studied populations is represented by a column. Numbers on the x-axis represent different observed alleles.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

SSR

All loci were highly polymorphic, with 17 to 53 alleles for each locus, totaling 227 alleles across all populations and loci. Of these, 114 were rare (frequency < 0.05) or absent and six common (frequency ≥ 0.05) across all seven studied populations. Except for one rare allele at locus PAAC19 (present only in core populations) and one allele at locus PGL15 (rare in all four peripheral populations and common in all three core populations), no pattern of allele distribution between peripheral and core populations was detected.

Relationships among populations

Previous studies reported that Norway spruce in Slovenia forms two distinct demes: (i) the Alpine and the Dinaric deme, concordant with recolonization from a glacial refugia identified by Tollefsrud et al. ([40]) in the Eastern Alps; and (ii) a second deme presumably located in the north Dinaric mountains. Such evidence was not confirmed in this study (Fig. 2). The analysis of population structure revealed that all populations belong to the same deme, both using isoenzyme and SSR markers, thus excluding the effect of hidden population structure on the findings of this study.

Fig. 2 - Lack of genetic structure among the analyzed populations as identified by Structure for SSR (left panel) and isoenzyme (right panel) markers. Each sampled tree is represented by a vertical line where a color denotes the part of the genome belonging to each of the K clusters (2 ≤ K ≤ 7).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Isoenzymes

Using our dataset based on 13 polymorphic isozyme loci, only one genetic group was identified by Bayesian clustering (Fig. 2). The lack of genetic differentiation among populations was confirmed by pairwise F_ST estimates between populations, which were small, ranging from 0.00 to 0.03 (Tab. S1 in Supplementary material), and were significantly different from zero (with α = 0.05 after Bonferroni’s correction) only for population pairs SD-TR (both inversion) and PP-GS (both upper tree line). Out of six pairs of populations with F_ST estimates higher than 0.02, five were pairs of peripheral populations. Population genotypic differentiation (as calculated by Genepop after the application of Bonferroni’s corrections) was significant for nine population pairs; in six cases these population pairs included a peripheral population (Tab. S1). Differentiation within peripheral populations, measured as global F_ST, was twice as high as in core populations (F_ST = 0.007). Inversion populations had the highest differentiation among each other (F_ST = 0.021 - Tab. 2). Nevertheless, the F_ST did not differ significantly between core and peripheral populations (p = 0.068).

SSRs

All analyzed populations belonged to the same genetic group based on the results of Bayesian clustering implemented in Structure (Fig. 2). The number of detected migrants between all populations after correction for size was moderately high, 6.6. Number of migrants estimated was slightly smaller for core (Nm = 4.8) than for peripheral populations (Nm = 6.3). With Nm being moderately high in both peripheral and core populations, ample gene flow likely counterbalances the contemporary effects of drift in the sampled populations, whereby all sampled populations seem to belong to one large panmictic population. Indeed, no differentiation was observed between core and peripheral populations using Bayesian clustering and differentiation was low for both isoenzymes and SSRs (F_ST-ISO = 0.002, F_ST-SSR = 0.002), although significant in the latter case (p = 0.013 - Tab. 2). High gene flow among and between core and peripheral populations may also be a consequence of the overlap in leaf phenology (Tab. 1), taken as a proxy of flowering phenology. As estimation of Nm based on the private alleles method assumes neutrality, and any kind of selection might lead to bias ([43]), Nm was only estimated for the neutral SSRs.

Consistent with high Nm, pairwise F_ST values were mostly low and ranged from zero to 0.21 (Tab. S2 in Supplementary material). No systematic difference in differentiation between peripheral and core populations was observed. Rather, population VP (inversion) differed significantly from all other populations except LIP (upper tree line). Of the seven significant cases of genotypic differentiation, three were among peripheral populations, three among peripheral and core populations and one among core populations (Tab. S2). Here, inversion population VP differed significantly only from the core population CV. Population VP also exhibited departure from Hardy-Weinberg equilibrium (see next section).

Despite the extensive gene flow, global F_ST estimates significantly differed from zero in peripheral and core populations, being higher in peripheral than core populations (p = 0.013 - Tab. 2). Similar to the isoenzyme dataset, global F_ST was almost twice as high in peripheral populations than in the core ones, and three times as high in inversion then in core populations.

Neither Bayesian nor classical methods revealed structuring of the studied populations using the two datasets, and gene flow was high among populations. Still, genetic differentiation was twice as high in peripheral then in core populations, in inversion populations three times as high as in core populations in both datasets. This is in concordance with the CPH and consistent with more than 70% of the 37 studies on plant species reviewed in Eckert et al. ([8]). For forest trees, the same pattern was observed for Thuja occidentalis ([25]) and Pinus contorta ([11]) but not in Liriodendron chinense ([44]) or Pinus strobus ([29]).

Genetic diversity estimates

Isoenzymes

Multilocus estimates of genetic diversity, i.e., allelic richness (A_R), observed heterozygosity (H_O), expected heterozygosity (H_E) and inbreeding coefficient (F_IS) averaged 2.204 ± 0.051, 0.175 ± 0.004, 0.183 ± 0.003 and 0.043 ± 0.024, respectively, for peripheral populations, and 2.143 ± 0.101, 0.168 ± 0.016, 0.185 ± 0.009 and 0.092 ± 0.055 for core populations, respectively. The number of private alleles and F_IS were lower by 62% and 74% in peripheral than core populations. Nevertheless, the null hypothesis of equal means of peripheral and core populations could not be rejected at α = 0.05 for neither these two estimates nor any of the other estimates. Hardy-Weinberg’s equilibrium was rejected for two out of 11 populations: the peripheral population VP and the core population TB. Isoenzyme diversity estimates in each population are presented in Tab. 3.

SSR

Multilocus estimates of genetic diversity, i.e., allelic richness (A_R), observed heterozygosity (H_O), expected heterozygosity (H_E) and inbreeding coefficient (F_IS) averaged 17.158 ± 0.219, 0.908 ± 0.0213, 0.926 ± 0.007 and 0.020 ± 0.021, respectively, for peripheral populations and 17.680 ± 0.624, 0.924 ± 0.003, 0.935 ± 0.006 and 0.014 ± 0.009 for core populations, respectively. Similar to isoenzymes, the null hypothesis of equal means of peripheral and core populations could not be rejected at α = 0.05 for none of the SSR diversity estimates nor for F_IS, number of rare and number of private alleles. Nevertheless, average F_IS was by 157% higher in peripheral than in core populations (Tab. 4), opposite to the pattern observed with isoenzymes. Higher F_IS in peripheral populations was also reported for Picea sitchensis ([12]). Despite corrections for null alleles, the Hardy-Weinberg’s equilibrium was rejected for population VP (p = 0.04). SSR diversity estimates for each population are presented in Tab. 4.

None of the datasets exhibited significant differences in genetic diversity estimates between core and peripheral populations for neither richness nor equitability measures. Moreover, neither the presence and quantity of rare nor the higher frequency alleles differed between core and peripheral populations, thereby rejecting the expectations based on the CPH hypothesis. Similar to Picea sitchensis ([12]), Pinus strobus ([29]) and Liriodendron chinense ([44]), genetic diversity in the studied peripheral Norway spruce populations was not lower than in core populations. In contrast, core populations of Thuja occidentalis ([25]) had higher diversity of richness measures. Still, the same study revealed similar expected heterozygosity in both core and peripheral populations, as found in our study.

Genetic diversity in relation to environmental factors

Correlation analysis was carried out on the isoenzyme dataset. The only significant correlation was found between F_IS and elevation (ρ = 0.637, p = 0.035 - Fig. 3); the populations from lower elevations showed higher F_IS than populations from higher elevations. Temperature and A_R tended towards a negative correlation (ρ = -0.516, p = 0.104), while no pattern was observed between A_R and precipitation (ρ = 0.345, p = 0.299) nor among other genetic and environmental indices. Mean values of genetic diversities A_R, H_E, H_O and F_IS did not differ significantly neither between populations growing on different bedrock nor for populations with different start of leaf flushing.

Fig. 3 - Plot of F_IS estimates vs. elevation.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Conclusions 

Our study was based on a limited number of populations from a geographically restricted area in the center of distribution of the southern lineage of Norway spruce in Slovenia. The studied populations were classified as core and peripheral, i.e., growing at the range margins (upper tree line, inversion). The added value of the study is: (i) the comparison of populations growing in a geographic area much smaller then those of similar studies carried out so far, therefore giving insights on the conformity to CPH over small distances; (ii) addressing periphery along an altitudinal gradient; and (iii) comparing the results based on two types of markers.

Our results showed that differentiation is higher in peripheral populations according to predictions from CPH, while genetic diversity is similar in peripheral and core populations. Moreover, the number of private and rare alleles were similar among core and peripheral populations (the number of private alleles in the isoenzyme dataset was much lower in peripheral populations, though not significantly different). This supports the evidence that peripheral populations along an altitudinal gradient may harbor alleles valuable for future adaptation ([9]).

The studied populations are geographically close to each other or close to other Norway spruce populations. Therefore, a high gene flow among populations is not surprising given the biology of the species. Higher differentiation of peripheral populations at the upper tree line and in inversions might be the result of lower effective population sizes and shorter lifespan at higher altitudes ([33]) or presumably at the bottom of vegetation inversions, leading to faster adaptation and consecutively, differentiation.

  Acknowledgements 

This study was done within the framework of research tasks of projects: V4-1438, V4-1614, COST Action FP1202 Strengthening conservation: a key issue for adaptation of marginal/peripheral populations of forest trees to climate change in Europe (MaP-FGR), and the Research Programme P4-0107 financed by the Slovenian Research Agency, V4-1438 and V4-1614 co-financed by the Ministry for Agriculture, Forestry and Food of the Republic of Slovenia. The authors are grateful to Dr. Monika Konnert (ASP, Teisendorf, Germany) for her supervision in the isoenzyme analysis. Data visualization was aided by Daniel’s XL Toolbox Add-in for Excel, version 6.52, by Daniel Kraus, Würzburg, Germany.

  References

  Supplementary Material

Fig. S1 - Geographic location of the sampled populations in Slovenia.
Fig. S2 - The Norway spruce population of Smrekova Draga (SD), Slovenia, located along an inverted vegetation gradient occurring in a Karst sinkhole.
Tab. S1 - Pairwise F_ST values and p-values for genotypic differentiation among Norway spruce populations based on 13 polymorphic isoenzyme loci.
Tab. S2 - Pairwise F_ST values and p-values for genotypic differentiation among Norway spruce populations based on six SSR loci.