Local Adaptation in European Firs Assessed through Extensive Sampling across Altitudinal Gradients in Southern Europe

Background Local adaptation is a key driver of phenotypic and genetic divergence at loci responsible for adaptive traits variations in forest tree populations. Its experimental assessment requires rigorous sampling strategies such as those involving population pairs replicated across broad spatial scales. Methods A hierarchical Bayesian model of selection (HBM) that explicitly considers both the replication of the environmental contrast and the hierarchical genetic structure among replicated study sites is introduced. Its power was assessed through simulations and compared to classical ‘within-site’ approaches (FDIST, BAYESCAN) and a simplified, within-site, version of the model introduced here (SBM). Results HBM demonstrates that hierarchical approaches are very powerful to detect replicated patterns of adaptive divergence with low false-discovery (FDR) and false-non-discovery (FNR) rates compared to the analysis of different sites separately through within-site approaches. The hypothesis of local adaptation to altitude was further addressed by analyzing replicated Abies alba population pairs (low and high elevations) across the species’ southern distribution range, where the effects of climatic selection are expected to be the strongest. For comparison, a single population pair from the closely related species A. cephalonica was also analyzed. The hierarchical model did not detect any pattern of adaptive divergence to altitude replicated in the different study sites. Instead, idiosyncratic patterns of local adaptation among sites were detected by within-site approaches. Conclusion Hierarchical approaches may miss idiosyncratic patterns of adaptation among sites, and we strongly recommend the use of both hierarchical (multi-site) and classical (within-site) approaches when addressing the question of adaptation across broad spatial scales.


Introduction
Local adaptation is the evolutionary process by which populations diverge toward different phenotypic and genetic optima in response to their local ecological conditions.Forest trees provide numerous examples of adaptive divergence across a variety of spatial scales, and local adaptation is supposedly a key process of tree populations' evolution and species' diversification [1].Indeed, forest tree species are often widely distributed across sharply contrasted environmental conditions, and local adaptation is favored in trees as a result of large population sizes and high levels of genetic variation for fitness-related traits [2][3][4][5][6][7][8][9][10][11][12].In particular, climate is one of the most important drivers of adaptation in forest tree populations [1,10,11,[13][14][15][16][17][18] and understanding the molecular bases of adaptation to climatic conditions is essential to accurately predict trees' responses to global climate change (GCC, [19]).European forests are facing enormous threats from rapid GCC with increasing frequency and intensity of summer drought, while considerable uncertainties exist about plants potential to respond to future warming and declining moisture availability [20].Although recent evidence suggests that pollen flow can connect populations more than 10 2 −10 3 km apart [21], escaping GCC through migration will also require substantial seed dispersal, which can be a limiting factor for many tree species.Those populations will thus have to cope locally with environmental changes, through individual physiological tolerance (i.e.phenotypic plasticity) in a proximate time (a time period corresponding to the individual lifespan), and through evolutionary change (i.e.genetic adaptation over generations, [22,23]).Notwithstanding large genetic variation and potentially fast adaptation, there is only anecdotal evidence that forest trees can genetically adapt to contemporary environmental change over a limited number of generations.In this context, studying species at the limits of their distribution range is particularly important for predicting the future evolution of species and their peripheral populations [10].
Outlier methods (i.e.F ST -based outlier detection tests) are frequently used to study local adaptation by identifying loci showing strong differentiation across populations [13,24].Their basic rationale is that loci influenced by natural selection toward different genetic optima across populations (i.e.divergent selection) are expected to be more differentiated than neutral loci, while loci subject to selection toward the same optimum across populations (i.e.homogenizing selection) are expected to be less differentiated than neutral ones [25][26][27][28][29].Despite its conceptual simplicity, the F ST -outlier approach suffers from several limitations, e.g. its sensitivity to heterogeneity in demographic histories among populations and lineages [30,31], subsequent complex spatial structuring [27,32,33] and variations in recombination and linked selection across the genome [34].To overcome these problems, it is recommended to apply different methods to reduce the risk of false positives [35][36][37] and to use replicated population comparisons [31,32].Two approaches are widely used to identify outliers under selection: the coalescent approach (FDIST) by Beaumont and Nichols [26] and the Bayesian method (BAYESCAN) by Foll and Gaggiotti [28].In particular, Bayesian modelling is very powerful to empirically calibrate complex models without a priori assumptions about the magnitude of parameters to infer [38].This approach has revolutionized the field of population genetics [39,40] and is now commonly used to analyze the genetic structure of populations [41][42][43] and to identify outliers for selection [28].Moreover, the Bayesian F ST -outlier selection test developed by Foll & Gaggiotti [28] has been shown to generally result in lower type I and II error rates than the original coalescent method by Beaumont and Nichols [26].
Although candidate gene approach has proved to be effective in detecting adaptive genetic differentiation in species with large genomes [18,[44][45][46][47], only few studies attempted to detect divergent selection to altitude using replicated population pairs spread across large spatial scales.Indeed, developing ad hoc models able to analyze replicated population pairs [31] is necessary to address the question of local adaptation in natural settings using replicated sampling designs.This is important since failing to account for demographic history [30,48] and, more generally, for hierarchical genetic structure [27,36,49,50] can lead to a lack of detection power or to the detection of false positives when analyzing replicated populations pairs simultaneously.It is thus necessary to develop methods adapted to large datasets and able to reduce the detection of false positives [27,32,36,[51][52][53][54][55][56].
The main objective of this work was to identify molecular evidence of adaptation to altitude in a conifer tree, Abies alba Mill., by sampling replicated populations pairs at study sites spread across its southern distribution range, where populations are expected to face strong selective pressures [10,57].To this end, a Bayesian model that explicitly considers both the replication of the environmental contrast 'low versus high' elevation and the hierarchical structure among replicated sites was developed.Pairs of populations were sampled at nine A. alba study sites plus one additional site of the congeneric species A. cephalonica Loudon.About 1600 individuals were genotyped using 273 SNPs within expressed sequences (ESTs), with the aim of deciphering the genetic outcomes of local selection in a keystone European forest tree species.

Ethics statement
Adult trees of A. alba were sampled in nine study sites in France, Italy, Bulgaria, and Romania (sites n°1 to 9).In France (sites n°1 to 5), no specific permission was required.Leaf sampling was nondestructive and carried out on public forest lands.In Italy (sites n°6 and 7), the sampling was authorized by the 'Ente Parco Nazionale del Gran Sasso e Monti della Laga'.In Bulgaria (site n°8), the sampling was carried out with permission and under the supervision of the 'Pirin National Park Administration'.In Romania (site n°9), the sampling was carried out with permission from the private owner of the forest site.
In addition, one natural population of A. cephalonica was sampled in the Peloponnese, Greece (site n°10).The permission was given by the General Directorate for the Protection and development of Forests and Rural Environment of the Greek Ministry for Environment and Energy.
Field studies did not involve endangered or protected species.

Biological model
Silver fir (A.alba Mill., Pinaceae) is a widespread European conifer occurring in mountainous regions of central and southern Europe, usually between 500 and 1500 meters [58].It is a keystone species of mountain forests with high ecological and economic value.A. alba populations belong to at least three main lineages (Pyrenees, Apennines and Balkan) that probably diverged during the Pleistocene as they remained at least partly isolated in multiple glacial refugia during Quaternary climatic cycles, including the last glacial maximum [59].In addition, adaptive variations have already been reported in this species for quantitative traits, isozymes, and candidate genes [18,[60][61][62][63]].
An additional study site of its congeneric species, Greek fir (A.cephalonica Loudon, [64]) was included in this study and treated separately from A. alba study sites to test whether local adaptation may be shared by two congeneric species.Greek fir is an endemic fir species widespread in the mountains of central and southern Greece from 400 to 1800 meters.Despite slight genetic divergence between the two fir species, they are closely related and able to produce fertile hybrids [59,65,66].

Sampling strategy and DNA extraction
Adult trees of A. alba were sampled in nine study sites located along the southern edge of the species distribution range, from the French Pyrenees to Romania (Fig 1, sites n°1 to 9), at two discrete elevations within each site ('low' versus 'high').One natural population of A. cephalonica located in the Peloponnese (Greece) was sampled according to the same strategy (site n°1 0).The location of the study sites and the sampling scheme are presented in Table 1.Material from the four French Alps sites is the same as in Roschanski et al. [62].Genomic DNA was extracted by LGC Genomics (Middlesex, United Kingdom).

SNP design, genotyping, blastX and functional annotation
A total of 763 SNPs and surrounding sequences derived from a transcriptome assembly [67] were sent to LGC Genomics (Middlesex, UK) for KASP assay design.Candidate genes were selected based on a specific annotation procedure described in details by Roschanski et al.
(2013) [67].This set of SNPs has further been enriched by including transcripts related to   [62] (for more details see these two related articles).Among the 763 SNPs tested for genotyping in the present study, 273 SNPs located within 177 transcripts were found to be amplifiable and polymorphic.These 273 SNPs were thus selected to genotype samples from the nine A. alba sites (n°1 to 9) plus the unique A. cephalonica site (n°10).Due to technical reasons, however, only 243 SNPs were successfully genotyped for the silver fir site of the Pirin Mountains (site 8, Bulgaria) and for the Greek fir site of Menalo Mt (site 10, Greece).After genotyping, a standardized (with uniform parameters) and deeper (including KEGG pathways and enzyme codes) annotation of the targeted transcripts was realized.For this, the targeted transcripts were blasted (blastX) and annotated using Blast2Go software [68] with minimum e-value = 10 −6 , annotation cut-off = 55 and GO Weight = 5.This extended functional annotation was consistent with the hypothesis of climatic adaptation for a majority of them, confirming thus that the targeted transcripts can be reasonably considered as candidates for selection.This dataset is thus a non-random genome-scan enriched in candidate genes which is more relevant than classical genome-scans [69] because SNPs located in expressed candidates are more likely to experience selection than SNPs randomly sampled from the entire genome.

Genetic structure analysis
Genetic structure was analyzed using STRUCTURE version 2.3 [42,70].Five independent runs for each K value ranging from 1 to 20 were performed under the admixture model, with a burn-in period of 15,000 followed by 25,000 iterations.The rate of change in L(K) across successive K values (ΔK) was calculated following Evanno et al. [71] using the web application 'STRUCTURE Harvester' [72].

Bayesian modelling, step 1: Bayesian allelic frequencies inference and G ST estimation
The Bayesian modelling approach is composed of two independent steps.The first step aims at inferring allelic frequencies and estimating pairwise G ST , while the second step aims at partitioning the G ST inferred through the first step into genome-wide and locus-specific parameters.Allelic frequencies (p i,m and q i,m ) at each bi-allelic marker (m) within each population i (i.e. each site×elevation combination) were inferred independently for each marker and population from allelic counts (n i,m and N i,m ) through the equations: The number of observations of a given allele (n i,m ) within a population i at a bi-allelic marker m is sampled from a binomial distribution with parameters p i,m (which is the allelic frequency in population i at marker m) and N i,m (which is the total number of allelic counts in population i at marker m).A Beta distribution was chosen as prior for allelic frequencies p i,m : p i,m ~β (1,1).

Bayesian modelling, step 2: Bayesian hierarchical outlier detection (HBM = Hierarchical Bayesian Model)
Because our sampling design consisted of 9 replicates of local elevations spread across a broad geographical scale, we designed a hierarchical Bayesian model (hereafter 'HBM') to partition pairwise G ST into genome-wide and locus-specific components under a two-level hierarchical model.
The median of G ST inferred through the first step was used as input in the second step.A logit transformation was applied to G ST values and a classical linear model was used to partition transformed G ST into genome-wide effects corresponding to the different levels of neutral genetic structuring (clusters, sub-clusters within clusters and elevation), and into locus-specific effects related to elevation: where k Clus(i,j) , k SubClus(i,j) and k Elev(i,j) are binary matrices describing whether the two populations (i,j) belong to different clusters, to different sub-clusters within a same cluster, or to different elevations: ) belong to different sub-clusters within the same cluster, k SubClus(i,j) = 0 otherwise; μ G , μ Clus , μ Subclus and μ Elev are genome-wide parameters describing the extent of genomewide differentiation between clusters, sub-clusters within clusters and between elevations.μ G is a global mean of differentiation that is inferred from all population pairs and captures the global extent of differentiation among all population pairs.μ Clus captures the genome-wide effect of belonging to different clusters, μ SubClus the genome-wide effect of belonging to different sub-clusters within a same cluster, and μ Elev the genome-wide effect of belonging to different elevations (whatever the cluster or sub-cluster).
θ Elev(m) are locus-specific parameters describing the extent of locus-specific genetic differentiation among populations belonging to different elevations.These parameters capture locusspecific effects of differentiation caused by elevation, assuming a common effect of elevation in all study sites whatever the cluster and sub-cluster.
Genome-wide and locus-specific parameters had non-informative priors: τ R, is the residual precision (i.e.inverse of the residual variance 1/σ² R ).
To make the model identifiable, the constraint Sθ Elev(m) = 0 was applied to the estimation of θ Elev(m) parameters.θ Elev(m) were sorted in decreasing order.Outliers were detected by fitting a posteriori a normal law from the inferred distribution of ordered(θ Elev(m) ) and by attributing a probability of being under selection to each marker.Outliers on the upper-tail of the distribution are significantly more differentiated than expected under neutrality (i.e.above the neutral background), while outliers on the lower tail are significantly less differentiated than expected under neutrality (i.e.below the neutral background).Outliers above the neutral background are thus good candidates for divergent selection between environmental contrasts repeated among the different study sites.Outliers below the neutral background may have different origins: they may alternatively reflect equilibrium allelic frequencies between elevations resulting from homogenizing selection, or rare alleles in one or few study sites caused by purifying selection or other neutral evolutionary processes.
6.1.Running HBM on simulated data and empirical assessment of false-discovery rate (FDR) and false non-discovery rate (FNR).To assess the statistical power and to explore the limits of our Bayesian modelling approach, HBM was applied to data simulated under different selection scenarios.Simulations were performed using the software SIMUPOP [75], an individual-based forward simulation of population evolution in Python language (http://simupop.sourceforge.net/Main/HomePage),S1 Method and Figshare doi: 10.6084/m9.figshare.1385219(https://figshare.com/articles/Populations_simulation/1385219).
The statistical power of HBM to detect divergent selection in the replicated population pairs was assessed through simulated datasets composed of 16 populations of N = 200 diploid individuals each (totaling N = 3200 diploid individuals) under a hierarchical model of migration, by varying both the number of selected loci and the selection strength (s = 1-ω with ω the fitness, S1 Method).
Two thousand iterations with a burn-in of 1000 were largely sufficient to reach chain convergence and to obtain informative posteriors.The model was thus applied to simulated populations with 1, 5 or 10 loci (out of 100) under divergent selection of different strengths (uniform among loci): s = 0.05 (weak selection), s = 0.075 (moderate selection), s = 0.1 (strong selection), s = 0.5 and s = 0.99 (very strong selection leading to differential allele fixation at different elevations).The model was also applied to simulated data composed of 10 loci under divergent selection with variable selection strengths (s) among loci: 0.06 to 0.15, 0.15 to 0.25 and 0.06 to 0.25.Finally, the power of the model was investigated in the cases of homogenizing selection and of a combination of homogenizing and divergent selection.In the former case, populations with 5 loci under divergent selection with selection strength of 0.1 were simulated.In the latter case, populations with 2 loci under divergent selection and 2 loci under homogenizing selection with selection strength of 0.1 were simulated.
False-discovery rate (FDR) and false non-discovery rate (FNR) were empirically assessed [76,77]: FDR = V/R where V is the number of false positives and R the number of rejected null-hypotheses, and FNR = T/(m-R) where T is the number of false negatives and (m-R) the number of accepted null hypotheses.
6.2.Running HBM on A. alba data.Each A. alba population was assigned to a particular cluster and sub-cluster according to the structuring detected using STRUCTURE (for K = 2 and K = 4, see results section, S1B Fig) , and the model was applied to the whole A. alba dataset using 10000 iterations with a burn-in period of 5000 iterations.The SNPs that were not genotyped in the nine A. alba study sites were discarded, resulting in a dataset composed of 243 SNPs successfully genotyped in all study sites (as described in section 3 'SNP design & genotyping').

Within-site outlier detection tests (FDIST, BAYESCAN and SBM)
Because the Bayesian model was developed to find replicated evidence of divergent selection in replicated population pairs, it was inherently unable to detect evidence of divergent selection particular to one or a few study sites.Thus, we tested for divergent selection between elevations within each site through classical coalescent (FDIST) and Bayesian approaches implemented in Arlequin and BAYESCAN respectively, plus a simplified (single-site) version of the Bayesian model introduced in this study (hereafter 'SBM').
The coalescent FDIST approach [26], implemented in Arlequin [78], was run under 20000 iterations assuming an island model.BAYEsCAN v.2.1 [28] was run with a burn-in of 50000 iterations, a sample size of 5000, a thinning interval of 10 (resulting in a total of 100000 iterations), and a prior odd of 10000 [32] indicating that for every locus the neutral model is much more likely than the model with selection (i.e.10000 neutral loci for every one under selection).
Loci displaying a significant effect of the α-component under a FDR threshold of 0.1 were retained.Lastly, SBM (a simplified, single-site, Bayesian model excluding genome-wide effects of differentiation among clusters and sub-clusters) was also applied to each site using 10000 iterations with a burn-in period of 5000 iterations.
Within-site outlier-detection tests (FDIST, BAYESCAN and SBM) were applied to both A. alba and A. cephalonica study sites.The power of within-site Bayesian methods (BAYESCAN and SBM) was also assessed through simulated data composed of two populations of N = 1600 diploid individuals each (totaling N = 3200 diploid individuals).Monomorphic loci within each site were discarded before within-site analyses.

Within-site Genotype-Environment Association tests (GEAs)
We assessed whether the genotypic frequencies were structured between elevations within each study site through genotype-environment association tests (see the methodology flowchart in S2 Fig) .For each study site: 1.An empirical distribution of the fixation index (F IS ) across loci in the entire population (low and high elevations confounded) was drawn, with: 2. n1 = 100 F IS -values were sorted in the empirical distribution of F IS , and n1 populations of N = 1000 diploids individuals were simulated under Hardy-Weinberg equilibrium corrected by the different F IS -values.
3. n2 = 100 sub-populations pairs of size 500 diploid individuals (corresponding to each elevation in our real dataset) were randomly sorted from the simulated populations of N = 1000 individuals.
4. The distribution of expected genotypic frequencies across the n1×n2 = 10000 iterations was then empirically drawn and corresponds to the expected distribution of genotypic frequencies at each elevation of each study site under the hypothesis of no structuring.
5. The cumulative probability of that distribution was drawn, and a p-value was attributed to observed genotypic frequencies at the different loci in the different sub-populations (elevations).It expresses the probability of a given genotype to be more or less present than expected.A negative association means that the genotype is less abundant at this elevation than expected under the hypothesis of balanced genotype frequency between elevations, while a positive association means that the genotype is more abundant at this elevation than expected under the hypothesis of balanced genotype frequency between elevations.
The test was applied to each study site separately.Monomorphic loci were discarded before analyses.The test is written in R and provided in appendix (S2 Method).The genetic structure remained concordant with the geographical distribution of the studied sites until K = 4, at which a secondary ΔK peak was detected, suggesting complex hierarchical structuring.At K = 4, the clusters detected at K = 2 were both split into two sub-clusters: the western Europe cluster into two sub-clusters corresponding to western (Pyrenees) and eastern (Alps) France, and the central-eastern Europe cluster into Italy and Balkan (Romania and Bulgaria) sub-clusters.K = 2 and K = 4 were thus retained for defining the assignment of each population to a given cluster and sub-cluster in the following hierarchical Bayesian approach (HBM), as described in the Materials and Methods.

Multi-site outlier detection under the HBM
2.1.Results using simulated data.Running HBM on simulated data revealed high power to detect divergent selection in the replicated population pairs under a hierarchical island model compared to within-site approaches.Table 2 shows a summary of the empirical assessment of FDR and FNR under the different selection scenarios at 5% and 1% p-value thresholds (see also S4-S7 Figs).FDR was always very low, even in the case of a combination of homogenizing and divergent selection and whatever the selection strength (s).No outliers were detected at 1% threshold when no selection was applied in spite of strong complex neutral genetic structure among populations, except for one false positive detected at a selection strength (s) of 0.5.FNR was low under scenarios of strong homogenizing and divergent selection (FNR = 0%), but was high under weak to moderate selection, even more when the number of loci was high and the selection strength was variable among loci: at a 1% threshold, the model missed all outliers when a selection strength of 0.05 was applied, and FNR varied between 0% and 8% at a selection strength of 0.075 depending on the number of loci under selection.While a 1% threshold may lack true-positives, a 5% threshold may produce false-positives.Thus, it is more conservative to retain only outliers detected at 1% threshold, as it is less damageable to miss outliers than to produce false-positives.For this reason, a threshold of 1% was retained.
2.2.Results obtained using A. alba data.The distributions of inferred G ST between populations are displayed in S8 and S9 Figs.As expected, the mean G ST decreased across levels of genetic divergence: the differentiation between populations belonging to different clusters and sub-clusters (G ST = 0.063, S8A and S8B Fig) was slightly higher than the differentiation Table 2. Power of the Bayesian approaches (HBM, SBM and BAYESCAN) to detect outliers for selection from simulated data.False-discovery rates (FDR) and false non-discovery rates (FNR) were empirically assessed under different scenarios of divergent and/or homogenizing selection, by varying the proportion of selected loci, and selection strength.

Multi-site analysis (Hierarchical model HBM)
Within   ).However, these two outliers were caused by rare variants rather than by homogenizing selection.SNP 26 was monomorphic except in sites 2 and 5 (Ventoux and Vesubie, East France): the minor allele was brought by three and one heterozygotes at low elevation in each site, respectively.SNP 111 was monomorphic except in site 5 at low elevation and was represented by one homozygote.These two SNPs were removed from within-site outlier detection tests in sites in which they were monomorphic.

Within-site outlier detection
3.1.Results using simulated data.SBM was applied to simulated data corresponding to individual study sites, and was less powerful than HBM as it failed to detect true positives in  2).
The detection power of SBM (with 1% threshold) and BAYESCAN (with prior odd 10000) was quite comparable under scenarios of uniform selection among loci, with minor differences only.However, BAYESCAN was less powerful than SBM under scenarios of variable selection among loci, as it produced more false-positives and more false-negatives.
The hierarchical approach implemented in HBM (under 1% threshold) was generally more powerful than within-site analyses, as the FDR was equal to 0% in all cases except one (one loci under very strong selection of strength s = 0.5).The FNR of HBM was also lower than the FNR of intra-site approaches, except in one scenario (10 loci under variable selection among loci, with selection strength ranging between s = 0.16 and s = 0.25) at which HBM detected only 6 outliers (i.e.four false-negatives) while SBM detected 8 outliers (i.e. two false-negatives) and BAYESCAN detected 11 outliers (i.e. one false-positive).
3.2.Results obtained using Abies data.A summary of the outliers detected using the different within-site approaches is provided in Fig 3, Table 3 and S1 Table.The coalescent approach (FDIST) detected a total of 97 outliers, all sites confounded: 64 outliers were below the neutral background and 33 outliers were above the neutral background (S1 Table ).The number of outliers in each study site detected by FDIST varied between five in the site 7 (Italy 'Colle dell'Abete') and twenty-two in site 9 (Romania).In site 8 (Bulgaria), the differentiation between elevations was too low to infer meaningful migration rates and to identify outliers through FDIST.BAYESCAN detected no outliers with a prior odd of 10000.

Within-site Genotype-Environment Association tests (GEAs)
GEAs were conducted in order to address the question of altitudinal adaptation at the genotype level.GEAs aimed at testing whether a genotype was more or less abundant at a given elevation than expected under the hypothesis of balanced genotypic frequencies between elevations.Among the 13374 associations tested (all sites, elevations, SNPs and genotypes confounded), 10306 (77.06%) were non-significant, while 3068 (22.94%) were significant, suggesting either a positive or a negative association of a given genotype (homozygote or heterozygote) at a given (low or high) elevation.Because a large proportion of outliers under divergent selection previously detected are expected to be false-positives, we further focused on the 12 consistent outliers for divergent selection detected twice by the two independent approaches FDIST and SBM respectively (Table 3).For all of them, significant GEAs were observed in the study sites where they were detected, with contrasted patterns between low and high elevations (S2 Table ).For example, at SNP 157, homozygotes (GG) were more abundant at low elevation than expected under the hypothesis of balanced genotypic frequencies between elevations (resulting in a positive association) in the study site 5 (where this SNP was detected as outlier for divergent selection).Homozygotes (GG) were also less abundant at high elevation than expected (resulting in a negative association) and heterozygotes (AG) were more abundant at high elevation than expected (resulting in a positive association).

Power of Bayesian modelling approaches (HBM, SBM and BAYESCAN) to detect homogenizing and divergent selection in replicated population pairs (simulations)
Bayesian approaches are known to be more stringent than methods based on coalescent simulations for identifying F ST -outliers [36,50]: even if Bayesian approaches may miss some true-positives, they often produce less false-positives than the coalescent approach by Beaumont and Nichols [26].This observation was confirmed in our study, where significantly more outliers were detected using FDIST than using BAYESCAN or SBM when applied to A. alba data (S1 Table ).
Table 3.Consistent outliers detected twice above the neutral background (by two different approaches).The first column describes the SNP number, the second column the SNP ID.The third column describes the study site in which the outliers were detected and the method used: FDIST (within-site coalescent method) and SBM under a 1% threshold (within-site Bayesian method).Study sites IDs are described in Table 1.The complete list of outliers detected, all methods confounded, is provided in S1 The power of BAYESCAN (with prior odd 10000) and SBM was also explored in more details through simulated data, revealing minor differences in the detection power of the two approaches (Table 2).These differences are probably caused by differences between the two methods, including mathematical formulation, programming, Markov Chain Monte Carlo sampling algorithm, and decision strategy.Especially, BAYESCAN and the Bayesian models developed here (HBM and SBM) use different differentiation indices: BAYESCAN estimates a locus-specific-population F ST close to Wright's F ST , while the models introduced here use pairwise Nei's G ST .Locus-specific G ST values were globally lower than F ST values (μ = 0.010 vs 0.017), but also more variable (σ = 0.013 vs 0.010) which probably facilitates the inference of locus-specific parameters (S11 Fig) .In addition, HBM and SBM partition pairwise G ST between populations into genome-wide and locus-specific parameters, while BAYESCAN partitions locus-population-specific F ST coefficients that measure the shared ancestry within each of the subpopulations.Moreover, the models encoded here do not use a composite likelihood to identify outliers and thus do not require locus-specific prior odds.Instead, there is no a priori assumption and a p-value is directly attributed to each marker from the inferred normal distribution of locus-specific parameters, leading to faster computation and allowing the detection of outliers below and above the neutral background.
Simulations revealed the high power of the hierarchical approach to detect both homogenizing and divergent selection in replicated population pairs: the false discovery rate was always very low, and outliers detected under a 1% threshold may be reasonably considered as truly under selection (Table 2 and S5-S7 Figs).However, HBM may fail to detect some loci under selection and the false-non-discovery rate may be high when the selection strength is weak (because the extent of differentiation is not strong enough to be outside the neutral background), particularly when many loci are selected and the selection strength is variable among loci (because the distribution of locus-specific parameters is empirically drawn and becomes wider when the number of selected loci increases).Nevertheless, HBM was more powerful than within-site Bayesian approaches (BAYESCAN and SBM) as it detected no (or few) false-positives with a FNR lower than intra-site approaches.In particular, it was able to detect replicated patterns of locus-specific divergence caused by moderate to strong selection strength that have not been detected when analyzing the different study sites independently through SBM (Table 2).This reinforces the idea that testing for local adaptation not only requires replicated population pairs at different study sites, but also the use of ad hoc models able to take into account replicated sites simultaneously with respect to their genome-wide genetic structure.On the other hand, HBM was designed to detect patterns of local adaptation replicated in the different study sites, but it was intrinsically unable to detect divergent selection particular to one or a few sites.Therefore, such hierarchical approach remains highly sensitive to heterogeneous selection among replicated population pairs, and testing for adaptation within each site remains essential to finely explore the process of adaptation within sites spread across broad geographical scales.

Altitudinal adaptation across the southern distribution range of Abies alba
Altitudinal gradients have attracted much attention because they display strong ecological variations (e.g.temperature, solar irradiance and precipitation) over short geographical scales [10,37,79,80], and many studies have provided evidence of local adaptation to elevation in a variety of plant species at both phenotypic [81][82][83] and molecular levels [18,19,61,[84][85][86][87][88].This study addresses the question of altitudinal adaptation across large spatial scales, from Western to Eastern Europe, using a dataset composed of 273 expressed SNPs located in candidate genes.
The analysis of neutral genetic variation showed a clear hierarchical genetic structure into clusters (K = 2) and sub-clusters (K = 4), as further confirmed by decreasing estimates of genome-wide parameters across geographical scales (μ Cluster > μ SubCluster > μ Elev ).This observation was concordant with a classical 'isolation-by-distance' model, in which the extent of neutral genetic differentiation increases with geographic distance because of limitations in gene flow caused by the geographic distance itself.
The hierarchical Bayesian model (HBM) detected only two SNPs below the neutral background (SNPs 26 and 111) but no common trend of divergent selection between elevations in sites spread across the whole southern distribution of silver fir.On the contrary, within-site analyses allowed the detection of many outliers (S1 Table ) but most of them are probably falsepositives, as the majority was detected in only one study site and by the coalescent (FDIST) approach only, and not validated by any of the Bayesian approaches (SBM or BAYESCAN).These genes may either be neutral or under too weak selection strength to be detected by the Bayesian approaches, thus making the distinction between true and false-positives extremely complicated.Especially, the demographic histories experienced by populations inhabiting the different sites (for example bottleneck or fast expansion) would also affect allele frequency spectra, resulting in the discovery of false-positives or false-negatives by selection tests [27,89,90].Even if Bayesian approaches [25,28] are known to be less sensitive to neutral processes and less prone to detect false-positives than the coalescent approach [36,50,91], they may remain sensitive to demographic processes [32,36,50,55,56,92].However, the confounding effects of demography and selection have not been fully explored yet, and it is impossible to conclude whether the genes showing a departure from neutral expectations in one site and only detected by a single approach are true-positives or not, without knowing if the different populations have experienced strong demographic changes.In addition, the relative measures of divergence used by these methods (F ST , G ST ) are sensitive to heterogeneous recombination rates in the genome.For example, centromeric or rearranged regions, where recombination is reduced, are expected to be more differentiated than the rest of the genome [34], and SNPs located in these regions are likely to be false-positives.However, the location of the targeted contigs and SNPs on the chromosomes of the species is unfortunately unknown.At last, the SNPs detected below the neutral background cannot be viewed as being involved in local adaptation caused by elevation and will not be considered thereafter.
Twelve SNPs (Table 3) may nonetheless be reasonably considered under divergent selection in one study site (SNPs 29, 58, 61, 65, 84, 99, 113, 157, 161, 203, 255 and 258), as they were detected above the neutral background by both FDIST and SBM and they were corroborated by contrasted GEAs between elevations at study sites in which they were detected.It is however unclear why these 12 SNPs were not detected by BAYESCAN, and several reasons may be invoked: (i) the selection strength may have been too weak to allow their detection by BAYES-CAN (as suggested by the results on simulated data), (ii) the dataset size may be too small to allow a reliable outlier detection, and/or (iii) a prior odd of 10000 may be too stringent to detect outliers in scans of expressed candidate SNPs [88].In fact, prior odds up to 10000 (that means that each SNP is ten thousand times more likely to be neutral than to be under selection) are commonly used to identify candidate loci in the context of genome wide association studies with millions of SNPs.However, each SNP analyzed through a genome scan of candidate genes is expected to be as much or more likely selected than neutral, and a prior odd of 10000 would probably lack true positives in this context.In addition, local demographic events may also produce false-positives, but the sensitivity of the different methods used (FDIST, BAYESCAN, SBM) has not been explored in this study.
These 12 consistent outliers (out of 273 analyzed) were particular to one study site, and the proportion of consistent outliers within each study site varied between 0% and 0.95% (S3 Table ).These estimates may however be biased by the size of the dataset.Indeed, this study focused on a set of 273 SNPs within 177 transcripts, which is far from being representative of the entire genome, and local adaptation may have targeted other genes not examined in this study.Moreover, the SNPs analyzed here are located in expressed candidate genes and therefore this set of SNPs is probably enriched in selected SNPs, which may result in an over-estimation of the true proportion of selected SNPs in the entire genome.Nevertheless, these estimates are comparable to those obtained by other experimental surveys of altitudinal adaptation in Fagus sylvatica.In the Mediterranean area, Jump et al. [93] detected one outlier locus out of 241 scored (0.4%) in the Mountains of Catalonia in north-eastern Spain, while Csilléry et al. [88] detected only three marginally significant outliers out of 546 scored SNPs (0.5%) in the Mont Ventoux study site in south-eastern France.Also, Pluess and Weber [94] detected 13 outlier loci out of 517 scored (2.5%) in the lowland forests of Switzerland.In Abies alba, adaptive variations have recurrently been reported for quantitative traits, isozymes, and candidate genes (sequences and SNPs data) in European wild populations [18,[60][61][62][63].In particular, Roschanski et al. [62] recently provided evidence of local adaptation to drought and cold tolerance in the French Alps at both molecular and phenotypic levels.They identified 16 outlier out of 267 SNPs (~6%) that showed patterns of divergent selection, and 8 genes containing SNPs at which allelic frequencies were correlated with bioclimatic variables.Among the SNPs they analyzed, one (0.37%) showed evidence of adaptation to altitude.In addition, two consistent SNPs detected here (SNP 65 in the Pyrenees and SNP 157 in the French Alps) were also detected as being involved in adaptation by Roschanski et al., providing thus additional evidence of local adaptation in A. alba.
Taken together, the results obtained through hierarchical and within-site approaches revealed only weak evidence of local adaptation to altitude, as idiosyncratic patterns among study sites were detected by the within-site FDIST and SBM approaches only (see for example variations in genotypic frequencies among sites at SNP 255, Fig 4).It is however not surprising to find site-specific evidence of selection not detected by the hierarchical (multi-site) approach, as the process of local adaptation to elevation is probably confounded with the processes of adaptation caused by large-scale factors not explored in this study (i.e.variations in environmental conditions between clusters, sub-clusters or sites).Indeed, the analyzed study sites may not be exactly 'true' replicates of the altitudinal contrasts as the range of ecological conditions might vary among them.In particular, the sampled populations occur at different elevations (Table 1), and the ecological conditions they experience may differ in many ways, including temperature, solar irradiance, length of the growing season, rainfall, soil type, water and nutrients availability and/or pathogens.In addition, the populations inhabiting the different study sites belonged to different lineages/gene pools as revealed by the complex biogeographic structure from West to East Europe: these populations may thus have adapted differently because different ancestral gene pools were exposed to selection.It is thus not surprising that different loci would have been under selection in the different study sites.
The 12 consistent outliers for divergent selection detected twice by within-site approaches and validated by significant GEAs were the strongest possible targets of natural selection, and are thus good candidates for exploring in detail the molecular bases of altitudinal adaptation in Abies.These SNPs are probably independent and not significantly linked to each other, or to other SNPs not analyzed in this study.They are located in different expressed genes (i.e.contigs, with size varying between 418 and 2340 bp, S4 Table) and we do not expect the presence of other -not genotyped-SNPs in the same transcripts or in other transcripts adjacent to those targeted in this study.Moreover, linkage disequilibrium is supposed to decay fast within genes in natural settings (within 100 bp), as suggested by Heuertz et al. [95] and Pavy et al. [96] in other coniferous species.For example, SNPs 157 and 158 were located in the same contig (428 bp apart) but only one of them (SNP 157) was detected by SBM and considered as a consistent outlier for selection.It confirms that linkage disequilibrium decreases sufficiently fast to prevent hitchhiking between close SNPs, and it proves the robustness of the Bayesian method SBM to potential linkage between a selected SNP and neutral SNPs located in the same or in adjacent contigs.
In this respect, functional annotation may help to identify the metabolic pathways as well as the biological process involved in local adaptation (S4 Table ).Especially, it needs to be stressed that six consistent outliers are located in transcripts involved in primary metabolic processes (SNPs 29, 58, 84, 113, 157 and 161) and three in a gene involved in nitrogen compound metabolic process (SNPs 65, 99 and 258).Primary metabolic processes -among which photosynthesis, carbon assimilation and sugar metabolism-are crucial for plants survival and development, and our results suggest that genes involved in primary metabolism pathways would be targeted by natural selection in European fir species.In particular, solar irradiance is a major constraint for plant development, and the processes related to photosynthesis and carbon assimilation have repeatedly been identified as being targets of natural selection in plant populations [10,[97][98][99][100][101].In addition, the primary products of photosynthesis and the activity of enzymes involved in sugar metabolism may vary with altitude as demonstrated by Kumar et al. [102] in grass species.Our results open the question of how metabolic variation in photosynthetic metabolism, sugar biosynthesis and their regulatory pathways may be genetically driven and structured by natural selection along altitudinal gradients.However, fleshing out our understanding on local adaptation to altitude would require a more detailed investigation of the adaptation process at the phenotypic level through quantitative genetic approaches (such as reciprocal transplants), and a fine exploration of the genetic architecture behind metabolic, physiological and morphological variations in wild populations (via association genetics and QTL analyses).

Fig 1 .
Fig 1. Genetic structure.Location of the study sites (A.alba, sites 1 to 9 and A. cephalonica, site 10) and genetic structure of A. alba populations revealed by STRUCTURE for K = 2 (top) and K = 4 (bottom).The map was created in ArcMap v.9.3 (ESRI.Redlands, CA).The European basemap is copyrighted by EUROSTATS (EuroGeographics for the administrative boundaries) and is available at: http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units.The black area shows the distribution range of A. alba (according to EUFORGEN 2009, http://www.euforgen.org).Study sites IDs are described in Table1.
STRUCTURE revealed a clear spatial structuring up to K = 4, with the highest ΔK at K = 2 followed by a secondary peak at K = 4 (Fig 1 and S3 Fig), suggesting that the most likely number of genetic clusters is 2: one group corresponding to western Europe (French sites) and one group corresponding to central-eastern Europe (Italian, Bulgarian and Romanian sites).
between populations belonging to the same cluster but different sub-clusters (G ST = 0.056 S8C and S8D Fig), which was subsequently higher than the differentiation between populations belonging to the same cluster and the same sub-cluster but different elevations (G ST = 0.013, S8E Fig).This resulted in decreasing genome-wide parameters inferred by HBM from broad to local scales: μ Cluster = 1.3 (with 95%CI = [1.28;1.35])> μ SubCluster = 1.14 [1.1;1.19]> μ Elev = -0.03[-0.06;-0.001].The extent of differentiation between elevations within sub-clusters was quite low compared to the extent of differentiation among populations belonging to different clusters (G ST = 0.063) and sub-clusters (G ST = 0.056): the G ST varied between 0.003 in site 9 (Romania) and 0.015 in site 4 (Issole, France) (S9 Fig).HBM did not detect outliers above the neutral background suggesting no tendency of divergent selection between high and low elevations replicated in the different sites.Two SNPs were detected below the neutral background: SNP 26 and SNP 111 (Fig 2 and S1 Table

Fig 3 .
Fig 3. Within-site outlier detection.Results of FDIST and SBM within each A. alba and A. cephalonica study sites.It is noteworthy that some outliers were detected in two study sites.doi:10.1371/journal.pone.0158216.g003

Fig 4 .
Fig 4. Idiosyncrasy between sites.Observed genotypic frequencies in the different sites at SNP 255 (genotypes CC, CT and TT).SNP 255 was detected by two within-site outlier detection methods (FDIST and SBM) in site 9 (purple line), but in no other site.In addition, significant GEAs were detected in site 9 (S2 Table).doi:10.1371/journal.pone.0158216.g004

Table 1 .
A. alba and A. cephalonica study sites and sampling design.