Research

Environmental versus geographical determinants of genetic structure in two subalpine conifers Elena Mosca1, Santiago C. Gonz
alez-Mart
ınez2 and David B. Neale1,3 1

Research and Innovation Centre, Fondazione Edmund Mach (FEM), Via E. Mach 1, S. Michele all’Adige, 38010, Italy; 2INIA, Forest Research Centre, Carretera de La Coru~ na km 7.5,

Madrid 28040, Spain; 3Department of Plant Sciences, University of California at Davis, Davis, CA 95616, USA

Summary Author for correspondence: Elena Mosca Tel: +39 0461 650956 Email: [email protected] Received: 22 March 2013 Accepted: 30 July 2013

New Phytologist (2014) 201: 180–192 doi: 10.1111/nph.12476

Key words: environmental gradient, isolation by adaptation, isolation by distance, landscape genetics, outlier locus detection.

 Alpine ecosystems are facing rapid human-induced environmental changes, and so more knowledge about tree adaptive potential is needed. This study investigated the relative role of isolation by distance (IBD) versus isolation by adaptation (IBA) in explaining population genetic structure in Abies alba and Larix decidua, based on 231 and 233 single nucleotide polymorphisms (SNPs) sampled across 36 and 22 natural populations, respectively, in the Alps and Apennines.  Genetic structure was investigated for both geographical and environmental groups, using analysis of molecular variance (AMOVA). For each species, nine environmental groups were defined using climate variables selected from a multiple factor analysis. Complementary methods were applied to identify outliers based on these groups, and to test for IBD versus IBA.  AMOVA showed weak but significant genetic structure for both species, with higher values in L. decidua. Among the potential outliers detected, up to two loci were found for geographical groups and up to seven for environmental groups. A stronger effect of IBD than IBA was found in both species; nevertheless, once spatial effects had been removed, temperature and soil in A. alba, and precipitation in both species, were relevant factors explaining genetic structure.  Based on our findings, in the Alpine region, genetic structure seems to be affected by both geographical isolation and environmental gradients, creating opportunities for local adaptation.

Introduction Retrospective studies of population genetic structure and adaptation in forest trees provide insights into how forests will respond to future environmental changes (Petit et al., 2008). Range expansions since the end of the last glaciations and the presence of specific adaptations in recently colonized areas have demonstrated that tree species can adapt rapidly in response to geographically variable selection (Holliday et al., 2010a; Keller et al., 2011). Gene flow has a pivotal role in this process, as it can both promote and prevent local adaptation (Kremer et al., 2012). Conifers are widespread in northern temperate forests. A long generation time and long dispersal distance of pollen and seeds make conifers valuable case studies for investigating genetic differentiation in the presence of gene flow (Savolainen et al., 2007), and for testing hypotheses about the relative importance of population isolation and environment-dependent selection in creating population genetic structure in long-lived organisms (Nosil et al., 2008; Andrew et al., 2012 for short-lived species). Species that grow across a wide range are subjected to distinct evolutionary forces, which may lead to local adaptation for ecologically important traits (Eckert et al., 2010a; Holliday et al., 2010a). A long history of common garden experiments has 180 New Phytologist (2014) 201: 180–192 www.newphytologist.com

repeatedly demonstrated that the interaction of different types of selection, large amounts of standing genetic variation and gene flow, and environmental heterogeneity promotes local adaptation in forest trees (Kremer et al., 2010; Alberto et al., 2013). Genomic regions affected by selection show specific signatures, such as high differentiation among populations (corresponding to an elevated estimate of genetic distance, FST) and a decrease in polymorphism within populations, and increased linkage disequilibrium (Schl€otterer & Harr, 2002; Beaumont, 2005; Keller et al., 2012). Genetic differences between populations are expected to increase with increasing distance between them, as a result of higher population isolation (Rousset, 2004). Moreover, environmental differences between sites are also expected to increase with distance. Therefore, genetic structure could also be affected by the association of geographical distance with a specific environmental factor, as demonstrated in both herbaceous species (Leimu & Fischer, 2008) and some forest trees (De Carvalho et al., 2010; Chen et al., 2012). One of the key goals of ecological genomic studies is the identification of loci that underlie local adaptation; this often involves identifying loci that show a different polymorphism pattern compared with the whole genome, that is, ‘outlier’ loci (Coyer et al., 2011). There are several strategies to detect such outliers, Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

New Phytologist depending on summary statistics and test assumptions (De Mita et al., 2013). Early methods were based on the allele frequency distributions (Lewontin & Krakauer, 1973; Akey et al., 2002), or on the FST distribution of gene diversity (Beaumont & Nichols, 1996). Both approaches focused on the identification of loci that show significantly higher or lower FST than the neutral expectation. A more recent development implemented a Bayesian method that considers a multinomial Dirichlet distribution (Balding, 2003). An alternative is to consider measures of genetic diversity, such as the expected heterozygosity (He), instead of FST. Schl€otterer & Harr (2002) developed such an approach for microsatellite markers that strictly follows the stepwise mutation model. Finally, under a model selection framework where selection effects can be included in or excluded from the model, Foll & Gaggiotti (2008) developed a Bayesian method that tests for specific population and locus effects. Recently these approaches have been increasingly applied to detect candidate genes potentially involved in adaptation in forest trees (Namroud et al., 2008; Eckert et al., 2010a). Several other studies focused on the association between genetic polymorphism and climatic variables (Eckert et al., 2010a,b; Holliday et al., 2010b; Prunier et al., 2011; Tsumura et al., 2012). For example, Prunier et al. (2011) found 26 single nucleotide polymorphisms (SNPs) associated with differences in temperature and precipitation in Picea mariana. To understand the complexity of adaptive processes and the main forces creating population genetic structure in natural populations, it is crucial to investigate the main selective forces acting upon the target species and environments (Joost et al., 2007). Alpine landscapes are characterized by heterogeneous topography as a consequence of the presence of physical barriers, such as valleys and high mountains, that can limit gene flow and influence genetic diversity within and among subalpine plant populations (Theurillat & Guisan, 2001; K€orner, 2003). As a consequence of this topographic heterogeneity, steep environmental gradients are also created in the Alps. Alpine plant communities are especially sensitive to change in climate (Theurillat & Guisan, 2001); therefore, temperature, precipitation and altitudinal gradients are expected to be the key forces shaping forest communities at high elevation. Nevertheless, genetic differentiation among populations can also be increased by population isolation and both topographic and climatic factors could be important (Nosil et al., 2008). In this context, it is relevant to determine whether isolation by geographical distance (IBD) or ecologically dependent reproductive isolation (i.e. ‘environmental isolation’ or ‘isolation by adaptation’ (IBA) – the restriction of gene flow as adaptive divergence increases; Nosil et al., 2008) is the main driver of genetic differentiation in alpine ecosystems. Under IBD, it is expected that species’ range fragmentation will cause an increase in genetic differentiation with the distance between populations, while under IBA the differentiation between populations is correlated to the relative influence of landscape and environmental variables on gene flow (Nosil et al., 2008; Andrew et al., 2012). In this study, we focus on the question as to whether geographical distance or habitat difference is more important for determining patterns of genetic variation in two subalpine forest trees, silver fir (Abies alba) and European larch (Larix decidua), using Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

Research 181

SNP genotyping data. The genotyped SNPs were located across a set of 150 genes (see details in Mosca et al., 2012a), involved in broadly different cellular mechanisms, which might underlie traits potentially adaptive to climate change, and are thus suitable to identify both general patterns of population structure and specific signatures of selection acting on some loci. In a previous study of four subalpine conifers (Mosca et al., 2012b), weak overall population genetic structure was found in A. alba and L. decidua, whereas some significant correlations were detected between genetic variation and environmental factors presumably as a result of natural selection. This study pointed to environmental selection as a more important force than historical isolation to explain population genetic structure in subalpine conifers. In the current study, we formally tested this hypothesis in different ways. First, we constructed geographical and environmental groups based on known topographical features and main environmental variables, respectively, in the Alps. Then, we tested for genetic differentiation among geographical and environmental groups using hierarchical analysis of molecular variance (AMOVA), and used different approaches to test for outlier loci based on these groups. Secondly, to test for overall geographical versus environmentally driven isolation, population genetic differentiation was tested for association with physical distance (IBD) and for correlation (after removing the geographical distance effect) with altitude, temperature and precipitation (IBA), three major environmental factors in Alpine ecosystems.

Materials and Methods Focal species, sampling, and SNP data Abies alba Mill. and Larix decidua Mill. were chosen for the present study because they are primary components of European alpine landscapes and have great ecological importance in these environments. Abies alba is a mountain species, broadly distributed throughout Europe. Its distribution is patchy, as a result of its demographic history, which is characterized by the presence of two main lineages, in the Alps and the Balkans, and by several migratory pathways from southern refugia following the end of the last ice age (Liepelt et al., 2009). Larix decidua is naturally distributed at high elevation (1000–2200 m) in the mountains of Central Europe; it may occur at the timberline in the Central Alps (Farjon, 1990). Its demographic history is characterized by a population expansion after a bottleneck, probably occurring during the last glaciations (Semerikov & Lascoux, 1999). Larix decidua was at that time a key pioneer species that colonized virgin soil after the glaciers retreated (Pleuss, 2011). For A. alba, 36 natural populations were sampled along its natural distribution in the Italian peninsula, from southernmost Serra San Bruno in Calabria to the Alps, to cover different ecological and climatic conditions (Fig. 1, Supporting Information Table S1). For L. decidua, 24 natural populations were sampled across its range in the Italian Alps, also covering a wide range of environmental conditions (Fig. 1; Table S1); two populations having < 10 successfully genotyped samples were removed. The populations used in this study are a subsample of those included New Phytologist (2014) 201: 180–192 www.newphytologist.com

New Phytologist

182 Research 13

2

6 7

32

33

24,23

29-31 9 34 20-22 28 3-5 10-12 35 8 18

19

16 36 1

15

26

25

17

14 27

0

300 km

(a)

0

200 km

1 4

12-14 9-11 18 7 6 16, 5

19-21

2

3 8 24

18

15 22 23

(b)

in Mosca et al. (2012b). In particular, populations were sampled across an altitudinal gradient ranging from 650 to 2197 m in A. alba and from 1123 to 2218 m in L. decidua. The average New Phytologist (2014) 201: 180–192 www.newphytologist.com

Fig. 1 Sampling locations across species distribution for Abies alba (a) and Larix decidua (b) across the Italian peninsula. The map was created using QGIS 7.1 (Quantum GIS Development Team, 2011) software.

annual temperature varied greatly across sampling locations in A. alba, from 3.71 to 12.08°C, while it ranged from 1.86 to 6.97°C in L. decidua. The cumulative annual precipitation Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

New Phytologist ranged from 606 to 1692 mm in A. alba and from 758 to 1646 mm in L. decidua (Table S1). Each sampled tree was georeferenced using Trimble GPS technology (http://www.trimble.com/) and marked; needles were collected from 25 to 65 individuals per population (n = 1108 in A. alba and n = 824 in L. decidua) for SNP genotyping. Generation of SNP data is explained in detail in Mosca et al. (2012b) and briefly outlined as follows. SNP discovery was carried out using a Sanger re-sequencing approach, using 800 PCR primer pairs from Pinus taeda in 12 individuals belonging to the two studied species (Mosca et al., 2012a). To obtain gene functional annotation, BLASTx and BLASTn analyses were performed on the P. taeda expressed sequence tag (EST) sequences in the National Center for Biotechnology Information, NCBI (http://www.ncbi.nlm.nih.gov/) and the Arabidopsis Information Resource, TAIR (http://www.arabidopsis.org/) databases, using published sequences of Arabidopsis (Notes S1 and Supplementary File 2 in Mosca et al., 2012a). A subset of the discovered SNPs (384 and 528 SNPs in A. alba and L. decidua, respectively) were then chosen for genotyping based on their Illumina (San Diego, CA, USA) design score, their coverage across the amplicon, the putative gene function, and the SNP allele frequency. The SNP genotyping was carried out at the University of California, Davis Genome Center using the Golden Gate platform (Illumina). SNP arrays were displayed on a Bead Array reader (Illumina) and analysed using GENOMESTUDIO V2009.1 software (Illumina). The successfully genotyped SNPs were selected using a minimum threshold of 0.25 for the GenCall50 score, a value of minor allele frequency > 1% and a Wright’s inbreeding coefficient (FIS) in the range of 0.35 to 0.35. In addition to these criteria, SNPs with a percentage of missing data higher than 20% were removed from the analysis. A suite of summary statistics was calculated to quantify the information content across loci in each sampling location and to further control genotyping quality. For each locus, the observed (Ho) and expected (He) heterozygosities and Wright’s inbreeding coefficient (FIS = 1 Ho/He) were calculated by population (Table S1). After this quality check, a total of 231 and 233 high-quality and polymorphic SNPs (conversion rates of 60.16% and 44.13%) remained for further analyses in A. alba and L. decidua, respectively. Environmental data In addition to geographical coordinates (see previous section for the sampling), we recoreded elevation, slope and aspect for each sampled tree. The average georeferenced point of each sampling location was used to identify the geographic position of the sampled population (Table S1). Climatic variables were obtained as described into Mosca et al. (2012b). Before the analysis, aspect was first transformed into ‘folded aspect’ about the north–south line, using the equation (folded aspect = 180 – |aspect – 180|), as suggested in McCune & Dylan (2002). Average monthly data and seasonal quarter averages were calculated for minimum, maximum and mean temperatures for each sampling site. Monthly cumulative precipitations and seasonal cumulative precipitations were calculated for each annual quarter, defined with 1 January Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

Research 183

to 31 March as the first quarter. The ‘growing degree days’ parameter with base temperature 5°C (GDD5), which is the threshold temperature for growth suggested for boreal conifers (Prentice et al., 1992, 2011; Sork et al., 2010), was computed for each species and site as the difference between average temperature and base temperature. Each sample site was also assigned to either the carbonate or silicate soil type category according to the Ecopedological Soil Map of Italy (European Communities, 2001). For the more densely sampled Trentino Province, a local soil map was used to assign the soil type category (Sartori & Mancabelli, 2009; Panagos et al., 2011). Soil type was used as a categorical variable by assigning a value of +1 to the carbonate soil and 1 to the silicate soil. To summarize highly correlated environmental data, a multivariate exploratory data analysis was performed using the FACTOMINER package (Husson et al., 2012) in R (R Development Team, 2011). A Multiple Factor Analysis for Mixed Data (AFDM) with both quantitative (climate) and categorical (soils) data was carried out for each species. The contribution of each variable to the axes is listed in Table S2A, where five axes (called ‘dimensions’) were taken into account. The correlation coefficient of each variable and the individuals’ coordinates on the axes were calculated (data not shown). The variable contributions were determined using the square cosine parameter and only variables with square cosine > 0.90 were considered to be variables with major contributions and used to define environmental groups (Table S2B; see next section). Geographical and environmental groups Based on current knowledge of the demographic history of subalpine trees (Semerikov & Lascoux, 1999; Liepelt et al., 2009; Koz
akov
a et al., 2011) and the location of physiographical barriers (e.g. the Po plain), four geographical groups were defined across the sampling range: one in the Apennines and three in the Alps, using the Dora Baltea and Adige rivers as physical borders (Table S1). In A. alba, two geographical groupings were tested: one considering two regions that corresponded to Alpine and Apennine populations and another considering four groups, one per geographical group described above. In L. decidua, which is absent from the Apennines, we only considered the three groups across the Alps (see Table S1). The main explanatory variables detected in the AFDM analysis (see previous section) were used to construct environmental groups, sometimes including populations that are geographically distant from each other in the same group. Nine environmental groups (Fig. S1) were defined in A. alba using March mean temperature (tmean_03) and cumulative precipitation of the first seasonal quarter (Q1_prec). Nine environmental groups were also defined in L. decidua, using the GDD5 measured in June (gdd_06) and the monthly cumulative precipitation of December (prec_dec). Genetic differentiation (AMOVA) Genetic differentiation among geographical groups and populations was examined for each species using hierarchical New Phytologist (2014) 201: 180–192 www.newphytologist.com

New Phytologist

184 Research

AMOVAs (Excoffier et al., 1992) in ARLEQUIN v. 3.5.1.3 (Excoffier & Lischer, 2010). The total variance was divided into three components: among geographical groups, among populations within groups, and within populations. The significance for each variance component was determined using a nonparametric permutation procedure (Excoffier et al., 1992). In a similar way as for geographical groups, hierarchical AMOVAs for each species using ARLEQUIN v. 3.5.1.3 (Excoffier & Lischer, 2010) were used to test for genetic differences among environmental groups. Finally, we tested the effect of environmental groups nested within the geographical groups using the function test.within of the HIERFSTAT package (Goudet, 2005) in the R environment. Detection of outlier loci Outlier locus detection analyses (i.e. detection of loci showing unusually low or high FST values) were conducted for both geographical and environmental groups. In this way, we were able to investigate whether geographical isolation or environment is driving adaptation in subalpine conifers at the molecular level. We used three complementary approaches to identify outliers. First, outliers were detected considering finite island (method 1) and hierarchical island (method 2) models using ARLEQUIN v. 3.5.1.3 (Excoffier & Lischer, 2010). The finite island model assumes an equal probability of migration between populations and mutation-drift equilibrium (Beaumont & Nichols, 1996), while the hierarchical island model assumes structured populations (Slatkin & Voelm, 1991). To correct the P-value obtained with ARLEQUIN for multiple testing, the function p.adjust was calculated in R (R Development Team, 2011), using Bonferroni’s correction. Secondly, a Bayesian method considering specific population and locus effects was used (method 3), as implemented in BAYESCAN v. 2.0 (Foll & Gaggiotti, 2008). The method is based on the estimation of two alternative models (including/excluding the effect of selection on single loci); their respective posterior probabilities are estimated using a Monte Carlo Markov chain. Model choice decision is performed using the Bayes factor. A prior odds equal to 10 was applied in both species with a false discovery rate (FDR) equal to 0.001. The gene functional annotation for outlier loci was obtained from the BLASTx of the P. taeda EST sequence versus Arabidopsis, while the SNP functional annotation (i.e. whether SNPs are noncoding, synonymous or nonsynonymous) refers to the SNPs submitted to the GenBank database (JQ440374– JQ445205) and EMBL website (HE663538–HE663608 and HE681087–HE681096) in the re-sequencing project (Mosca et al., 2012a). Isolation by distance (IBD) versus isolation by adaptation (IBA) Mantel tests and multiple regression matrix analyses were both applied to determine whether patterns of overall genetic differentiation in A. alba and L. decidua are attributable to IBD or to environment-driven selection. New Phytologist (2014) 201: 180–192 www.newphytologist.com

First, IBD was investigated using a Mantel test to correlate pairwise physical spatial distances and pairwise genetic differentiation (FST), calculated according to Weir & Cockerham (1984). Both matrices were generated with SPAGEDI 1.3a (Hardy et al. 2009) and scaled prior to the Mantel test. Mantel tests and partial Mantel correlations were also used to investigate the association between genetic distance and environmental variables (IBA), that is, elevation (E), annual mean temperature (T), annual cumulative precipitation (P) and soil type (S), with (to remove purely geographical effects) and without introducing spatial distance as a covariate. IBA analyses were performed considering both all populations and populations within geographical groups (see ‘Geographical and environmental groups’ section). All distance matrices were scaled prior to the analysis. All tests were performed with the mantel function in the ECODIST package (Goslee & Urban, 2007) in R (R Development Team, 2011). Secondly, to further investigate the effects of environment on genetic distance and test for IBA, environmental distance matrices were constructed (Table S1). Sampling site elevation (E), temperature (T), and precipitation (P) were used directly, while two categories (silicate or carbonate) were used for soil type (S), assigning 0 to silicate soils and 1 to carbonate soils. Then, for each variable, environmental distances among sites were computed as the Euclidean distance between populations with 10 000 simulations. All matrices were scaled before the analysis. Environmental distance matrices were used as a predictor in Multiple Regression Distance Matrix (MRDM) analyses, in which genetic distance was the response variable, both using and not using physical spatial distances as a factor. All matrix correlations were performed with the ECODIST package (Goslee & Urban, 2007) using the MRM function in R (R Development Team, 2011). The MRDM analyses were performed both using all sampled populations in each species and for each of the geographical groups defined in the previous section.

Results The number of SNPs successfully genotyped in each species was 231 SNPs across 150 genes in A. alba and 233 SNPs across 151 genes in L. decidua. The expected heterozygosity (He) calculated for each sampling site (population) ranged between 0.247 and 0.305 for A. alba and between 0.170 and 0.301 for L. decidua, while the observed heterozygosity (Ho) ranged from 0.237 to 0.315 in A. alba and from 0.180 to 0.406 in L. decidua (Table S3). Levels of linkage disequilibrium (LD) were low (Fig. S2), as observed in other conifers (Pavy et al., 2012). Genetic differentiation among geographical and environmental groups AMOVA was used to test for significant genetic differentiation among geographical and environmental groups in each species. For geographical groups (Table 1a) in A. alba, both genetic differentiation among groups (FCT) and genetic differentiation among populations within groups (FSC) were highly significant (P < 0.0001), with a percentage of variation explained of 2.18% Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

New Phytologist

Research 185

Table 1 (a) Analysis of molecular variance (AMOVA) based on single nucleotide polymorphism (SNP) data assuming a hierarchical geographical population structure with two regions and four geographical groups in Abies alba and with three geographical groups in Larix decidua; (b) AMOVA based on single nucleotide poloymorphism (SNP) data assuming population structure with nine environmental groups in A. alba and in L. decidua

Species

Code

Source of variation

df

Sum of square

Variance components

% of variation

P

Fixation index

(a) A. alba

Regions

A. alba

Geo-groups

L. decidua

Geo-groups

Among groups Among populations within groups Within populations Among groups Among populations within groups Within populations Among groups Among populations within groups Within populations

1 34 21 808 3 32 21 802 2 19 1626

347.818 2580.713 57831.32 732.914 2195.617 57831.318 549.221 973.052 29702.47

0.6096 0.8036 26.528 0.36306 0.6895 26.528 0.63047 0.45297 18.2672

2.18 2.88 94.94 1.32 2.50 96.18 3.26 2.34 94.40

< 0.0001 < 0.0001 < 0.0001 < 0.0001 < 0.0001 < 0.0001 < 0.0001 < 0.0001 < 0.0001

0.0218 0.0294 0.0506 0.0132 0.0253 0.0132 0.0559 0.0242 0.0326

(b) A. alba

Env-groups

L. decidua

Env-groups

Among groups Among populations within groups Within populations Among groups Among populations within groups Within populations

8 27 2180 8 13 1626

872.46 2056.07 57831.32 912.343 609.930 29702.47

0.1676 0.7976 26.528 0.3826 0.3965 18.2672

0.61 2.90 96.49 2.01 2.08 95.91

0.0020 < 0.0001 < 0.0001 0.0020 < 0.0001 < 0.0001

0.0061 0.0292 0.0351 0.0201 0.0212 0.0409

The population assignations to the groups are reported in Supporting Information Table S1 and described in the Materials and Methods section.

and 2.88% for the two regions and 1.32% and 2.50% for the four geographical groups, respectively. In L. decidua, genetic differentiation among groups (FCT) and genetic differentiation among populations within groups (FSC) were also highly significant (P < 0.0001), with percentages of variation explained of 3.26% and 2.34%, respectively. With respect to environmental groups (Table 1b), AMOVAs in both A. alba and L. decidua were highly significant (P < 0.0001) for genetic differentiation among populations within groups (FSC) and less significant (P = 0.002) for differentiation among groups (FCT), explaining substantially lower percentages of variance in A. alba (0.61%) than in L. decidua (2.01%). Finally, genetic differentiation of environmental groups nested within the geographical groups was not significant. Detection of outlier loci Three methods were applied to detect outlier SNP loci in each species. Using geographical groups, the number of outliers detected assuming a finite island model (method 1) was 15 in A. alba and 34 in L. decidua, whereas the hierarchical island model (method 2) detected 12 and 23 outliers in A. alba and L. decidua, respectively (Table S4). These numbers were greatly reduced after adjusting the P-value with Bonferroni’s correction (Tables 2, S4). In A. alba, after multiple test corrections, method 1 detected four outlier SNPs, while only one and two outlier loci were found with method 2, considering two regions and four geographical groups, respectively. While two loci were detected as outliers by both methods (2_6313_01_Abal_160 and CL3116Contig1_03_Abal_118), two other loci (0_7009_01_Abal_212 and 0_5361_01_Abal_287) were only found with the island model. In L. decidua, after multiple test corrections, more outliers were found with method 1 (six loci) than with method 2 (two loci). In Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

this species, the outlier loci were characterized by a higher value of FST than those in A. alba. To identify outliers for environmental groups, only hierarchical island models (method 2) can be applied. In A. alba, 17 outlier loci were detected among nine environmental groups (Table S4) and five loci were still significant after Bonferroni’s multiple testing corrections (Table 2). These loci were also considered outliers in analyses based on geographical groups or the island model, with the exception of locus CL1455Contig1_06.Abal.152. In L. decidua, 39 outlier loci were found and seven remained significant after applying Bonferroni’s correction (Table 2). Four of these loci were specific for analyses based on environmental groups (0_7810_01.Lade.525, 0_9284_02.Lade.470, CL1045Contig1_01.Lade.380 and CL1077Contig1_02.Lade.225). Finally, the rate of outliers detected with Foll and Gaggiotti’s Bayesian approach (method 3) ranged from 1.73% in A. alba to 6.86% in L. decidua (Table 3; Fig. S3). In both species, the majority of these outliers (three loci out of four in A. alba and 12 loci out of 16 in L. decidua) were only detected with the Bayesian approach (Tables 3, S5) and not with other methods. Functional annotation of outlier loci The total number of outlier loci detected with the different methods ranged from eight loci in A. alba to 24 loci in L. decidua (Table S5), corresponding to 3.46% and 10.30% of the total number of SNPs analysed in each species, respectively. The functional annotation of the outlier genes describes the putative protein they code for; for each SNP, annotation type is also reported (NC, noncoding; SY, synonymous; NS, nonsynonymous) (Table S5). Even though most SNPs were found in noncoding regions or were synonymous, one SNP in A. alba was nonsynonymous New Phytologist (2014) 201: 180–192 www.newphytologist.com

New Phytologist (2014) 201: 180–192 www.newphytologist.com

0.355

FST

1.0E-07

P

0.3650

0.3748

0.1308

FST

0.0836

He

1.0E-07

7.3E-05

P

Three geographical groups

0.359

He

Two regions

***

**

Pp.adj

***

Pp.adj

0.183

0.137

FST

1.0E-07

1.8E-04

P

0.157 0.073 0.090 0.228 0.143 0.232 0.117

He

0.248 0.270 0.172 0.194 0.000 0.163 0.301

FST

1.0E-07 1.0E-07 1.0E-07 1.0E-07 2.0E-04 1.0E-07 1.0E-07

P

Nine environmental groups

0.283

0.197

He

Four geographical groups

*** *** *** *** * *** ***

Pp.adj

***

*

Pp.adj

0.146 0.115 0.160

0.116 0.104

FST

NM_120865 NM_111095 NM_106106 NM_202043 NM_001035973 NM_121136 NM_112228

P. taeda BLASTx

0.200 0.377 0.275

0.137 0.504

He

*** *** ***

*** ***

Pp.adj

NM_118832 NM_106473 NM_122008

NA NM_104402

P. taeda BLASTx

NA Kelch repeat-containing F-box family protein Unknown protein Protein binding Member of RAN GTPase gene family

Protein description

Beta-glucuronidase Flavodoxin family protein GDSL-motif lipase/hydrolase family protein Phosphate translocator-related Protein binding/protein homodimerization/transcription repressor Histone H3 Endomembrane protein 70

Protein description

1.0E-07 1.0E-07 1.0E-07

9.8E-05 1.0E-07

P

Nine environmental groups

NA1 NA NA NA1 NA SY NA

SNP code

NA NA NS

NA SY1

SNP code

The P-value is adjusted (p.adj) using Bonferroni’s correction for multiple testing. Only significant analyses after Bonferroni’s correction are shown. The SNP code is as follows: NA, no annotation; SY, synonymous; NS, nonsynonymous. The protein description refers to the BLASTx of the Pinus taeda expressed sequence tag (EST) sequence versus Arabidopsis. 1 Loci moderately associated with environment in Mosca et al. (2012b). ***, Pp.adj-value < 0.0001; **P-value < 0.001; *, Pp.adj-value < 0.05.

(b) 0_17790_01_159 0_7001_01_260 0_7810_01_525 0_9284_02_470 CL1045Contig1_01_380 CL1077Contig1_02_225 CL1634Contig1_03_108

Larix decidua Locus

2_6313_01_160 CL1455Contig1_06_152 CL3116Contig1_03_118

(a) 0_5361_01_Abal_287 0_7009_01_212

Abies alba Locus

Table 2 Outlier loci found in Abies alba (a) and Larix decidua (b) with the hierarchical island model (method 2) with two regions or three/four geographical groups (in L. decidua and A. alba, respectively), and with nine environmental groups

186 Research

New Phytologist

Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

New Phytologist

Research 187

Table 3 Detection of outliers using the Bayesian approach (method 3) with false discovery rate (FDR) = 0.001 and prior odds equal to 10 Functional annotation

Species

SNP

Abies alba

0_7009_01_Abal_212

Larix decidua

ID 46

Prob

log10 (PO)

1

1000

1_6493_01_Abal_226 CL4354Contig1_01.Abal.147

119 60

0.998 0.997

CL4354Contig1_01.Abal.202

185

1

Alpha

2.795 2.584 1000

FST

Sign

P. taeda blastx

Protein description

1.158

0.125

***

NM_104402

1.724 1.616

0.01 0.011

*** ***

NA NM_115712

1.366

0.147

***

NM_115712

Kelch repeat-containing F-box family protein NA Protein serine/threonine phosphatase (PP2A-3) Protein serine/threonine phosphatase (PP2A-3) Glycoprotease M22 family protein Serine-tRNA ligase CpDNA Beta-glucuronidase Unknown protein CF9 mRNA Peroxidase Phosphate translocator-related UDP-glucuronate 4-epimerase/catalytic (GAE1) Heat shock protein 101 Putative LRR receptor-like serine/threonine protein kinase MRH1 GTP-binding protein GTP-binding protein Histone H3 Exostosin family protein Endonuclease exonuclease phosphatase family protein

0_11772_01.Lade.137

25

0.997

2.493

1.180

0.111

***

NM_118398

0_14221_01.Lade.224 0_14591_02.Lade.108 0_17790_01.Lade.159 0_18644_02.Lade.469 0_5038_01.Lade.226 0_6659_01.Lade.101 0_9284_02.Lade.470 2_10352_02.Lade.72

72 212 216 248 102 38 138 16

1 1 0.987 0.998 1 1 1 1

1000 1000 1.894 2.698 1000 1000 1000 3.398

1.474 1.325 1.050 1.204 1.617 1.628 2.369 1.470

0.140 0.124 0.100 0.113 0.155 0.156 0.258 0.140

*** *** *** *** *** *** *** ***

NM_179316 AP000423 NM_120865 NM_100301 NM_104644 NM_103577 NM_202043 NM_119190

2_6317_01.Lade.233 2_8011_02.Lade.468

259 255

1 1

1000 3.398

2.246 1.375

0.239 0.130

*** ***

NM_106091 NM_116232

2_9465_01.Lade.228 2_9465_01.Lade.430 CL1077Contig1_02.Lade.225 CL3832Contig1_05.Lade.89 CL4776Contig1_03.Lade.552

215 179 181 35 60

1 1 0.994 0.997 0.997

1000 1000 2.191 2.552 2.467

2.408 1.863 1.104 1.256 1.216

0.266 0.185 0.104 0.118 0.114

*** *** *** *** ***

NM_123975 NM_123975 NM_121136 NM_118384 NM_115718

SNP code SY1 SY SY NC SY NA NA NA1 NA SY NA NA1 NA NA SY

NA NA SY NA NC

The single nucleotide polymorphism (SNP) code is as follows: NA, no annotation; NC, noncoding; SY, synonymous; NS, nonsynonymous. The protein description refers to the BLASTx of the Pinus taeda expressed sequence tag (EST) sequence versus Arabidopsis. 1 Loci moderately associated with environment in Mosca et al. (2012b).

and located in a gene (CL3116Contig1_03) encoding the GTPbinding nuclear protein Ran-1 in Arabidopsis. In both species, the remainder of the outlier loci were silent and located in genes encoding unknown proteins or proteins broadly involved in several metabolic processes (Table S5; see also the Discussion section). Isolation by distance (IBD) versus isolation by adaptation (IBA) The Mantel test based on pairwise genetic and geographical distance matrices was positive and highly significant in A. alba and L. decidua (Table 4; P = 0.00001 in both species). This result suggested that IBD may account for most of the differentiation among populations found in both species. By geographical group, IBD was not significant for A. alba in geographical groups 1 (western Alps) and 4 (Apennines), whereas it was significant in the other groups (central and eastern Alps). Conversely, in L. decidua the genetic distance by geographical group was only correlated with the spatial distance in the central Alps. With respect to environmental distances (not corrected by geographical distance), overall correlations with genetic distance were positive in A. alba for temperature (P = 0.00140) and Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

precipitation (P = 0.04486), but not for elevation. The correlations calculated within geographical groups were positive for soils (P = 0.00026) in geographical group 2 (central Alps), and soils (P = 0.02422) and precipitation (P = 0.00007) in geographical group 3 (eastern Alps). In L. decidua, the (noncorrected) Mantel test including all populations was significant only between genetic distance and precipitation (P = 0.00014). The same result (P = 0.00007) was found within geographical group 2 (central Alps), as well as a weaker correlation with the same variable (P = 0.0170) in geographical group 1 (western Alps). Partial Mantel correlation was also used to examine the contribution of each environmental variable to the differentiation among populations when the geographical distance was taken into account. In A. alba, the overall correlation between genetic distance and both temperature and precipitation remained significant when taking into account geographical distance (P = 0.00001); moreover, in this case, the correlation was also significant with elevation (P = 0.00001). In L. decidua, after correcting for geographical distance, the correlation between genetic distance and precipitation was still significant considering both all populations (P = 0.0010) and geographical group 2 (central Alps; P = 0.00152). New Phytologist (2014) 201: 180–192 www.newphytologist.com

New Phytologist

188 Research

Table 4 Mantel and partial Mantel correlation coefficients used to test for association of genetic distance (F) considering all populations (‘All sites’) and only populations within geographical groups (‘Geo-groups’) with physical distance between pairwise populations (D), elevation (E), annual mean temperature (T), annual cumulative precipitation (P) and soil type (S) Geo-group 1 (western Alps)

All sites Abies alba Mantel test

r

F~D F~E F~T F~P F~S F~E|D F~T|D F~P|D F~S|D

0.5491 0.0339 0.2140 0.0902 0.0124 0.5492 0.5184 0.5294 0.0183

P1

r

P1

0.0000 0.6340 0.0014 0.0448 0.3570 0.0000 0.0000 0.0000 0.2952

0.0306 0.5532 0.3753 0.0488 0.1913 0.2341 0.4464 0.1608 0.2341

0.4504 0.9665 0.7348 0.5032 0.5654 0.6174 0.8346 0.6167 0.6175

Geo-group 1 (western Alps)

All sites Larix decidua Mantel test F~D F~E F~T F~P F~S F~E|D F~T|D F~P|D F~S|D

r 0.8519 0.0081 0.0546 0.5628 0.0147 0.0715 0.0107 0.4416 0.0424

Geo-group 2 (central Alps)

Geo-group 3 (eastern Alps) P1

r 0.9228 0.0771 0.0261 0.0813 0.4963 0.1488 0.0448 0.2163 0.2292

r

P1

0.0000 0.7957 0.2192 0.0001 0.4713 0.7427 0.4948 0.0010 0.2765

0.6227 0.1605 0.1594 0.4962 NA 0.2402 0.4479 0.2366 NA

0.0164 0.3488 0.3669 0.0171 NA 0.9669 0.2178 0.8509 NA

r 0.6843 0.1024 0.1976 0.6174 0.0611 0.0610 0.1311 0.2666 0.0567

0.8382 0.1088 0.3611 0.8078 0.3305 0.1076 0.3554 0.0364 0.0007

0.0000 0.7395 0.0165 0.0001 0.0242 0.7332 0.0167 0.3528 0.4061

P1

r 0.0186 0.0524 0.3143 0.0525 0.1767 0.0569 0.3398 0.0634 0.1851

0.4291 0.5486 0.8896 0.4174 0.6399 0.5504 0.8898 0.3655 0.6380

Geo-group 32 (eastern Alps)

Geo-group 2 (central Alps)

P1

P1

r

0.0000 0.2593 0.5584 0.7873 0.0003 0.1337 0.3595 0.9885 0.0765

Geo-group 4 (Appenines)

P1

r

P1

0.0000 0.1466 0.0360 0.0001 0.2527 0.2473 0.0944 0.0152 0.6536

NA NA NA NA NA NA NA NA NA

NA NA NA NA NA NA NA NA NA

All distance matrices were scaled before the analysis. Bold font indicates significant values, r is the Pearson coefficient of correlation and P is the associated P-value. 1 One-tailed P (null hypothesis: r ≤ 0). 2 Not calculated; only two populations belonged to Geo-group 3.

These results were further confirmed by MRDM analyses (Table S6), showing a significant overall association between genetic distance and both temperature and precipitation in A. alba, as well as for the eastern Alps (geographical group 3). The MRDM analysis also found a strong significant effect of soils for A. alba in the central Alps (geographical group 2) and marginally in the eastern Alps (geographical group 3). In L. decidua, the analysis confirmed the overall association between genetic distance and precipitation, and also for the western and central Alps (geographical groups 1 and 2). Moreover, a weak significant association with temperature was found in the central Alps (geographical group 2).

Discussion This study presents a comprehensive investigation of patterns of genetic variation across environmental gradients in A. alba and L. decidua in the Alps, with a focus on the causes of genetic divergence among populations. Weak but significant genetic structure was found in both species. Among the potential outliers detected, two loci (out of five) in A. alba and two (out of seven) in L. decidua were related to both geography and specific environmental variables (based on environmental groups); four and one additional loci were outliers in relation only to New Phytologist (2014) 201: 180–192 www.newphytologist.com

environment in L. decidua and A. alba, respectively. A stronger effect of IBD versus IBA was found for population genetic structure in both species. Finally, we demonstrated how the effect of environment can be partially separated from other confounding factors, such as geography. Specifically, together with historical isolation, both temperature and precipitation in A. alba and only precipitation in L. decidua were identified as relevant factors explaining population genetic structure in these keystone subalpine trees. The presence of a geographical effect on genetic variation was examined by assigning sampling sites to up to four geographical groups (eastern, central and western Alps, and Apennines). In A. alba, higher genetic differentiation was found in the two regions sampled (Alps and Apennines) than in the four geographical groups. Populations from the Apennines probably belong to a different gene pool from populations from the Alps (Liepelt et al., 2009), which would explain this result; moreover, within the Alps, some genetic structure was found but not clearly corresponding to geographical regions. In L. decidua, genetic differentiation was stronger (FCT = 0.0559), showing that genetic diversity was primarily distributed according to geographical regions. Despite the significant differences detected among geographical regions and among populations within geographical region, the main source of variation occurred within population Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

New Phytologist in both species, which is typical for forest trees (M€ uller-Starck et al., 1992; Savolainen et al., 2007). The case for environmental selection was formally tested by defining nine distinct environmental groups in each species, using the principal explanatory variables obtained from an AFDM analysis. Selected climate variables were related to winter precipitation (both species) and to March temperature (A. alba) or GDD5 measured in June (L. decidua). The AMOVA for environmental groups showed low but significant among-groups genetic differentiation, with higher values for L. decidua (2.01% of variation) than A. alba (0.61%). In addition, outlier loci exclusive to environmental groups included one locus in A. alba and four in L. decidua. These findings suggested the presence of environmental discontinuities that might have contributed to shaping population genetic structure, as well as fostering local adaptation. In A. alba, our results suggested differences in the response to climate across the Alps, as previously described based on tree ring growth (Carrer et al., 2010). In L. decidua, our findings suggested that environmental groups based on December precipitation and GDD5 are relevant for adaptive processes. Even when this species loses its needles, minimizing winter desiccation damage, these variables may still be important in regulating species growth phenology. Similar results were found for bud burst in relation to temperature in Picea sitchensis (Mimura & Aitken, 2010) and for water availability in a wild relative of Arabidopsis (Lee & Mitchell-Olds, 2011). The stronger genetic structure and higher number of outlier loci found in L. decidua for environmental groups suggest that adaptation in this species is more environment-dependent than in A. alba. However, these findings could also be related to the higher number of loci potentially under selection found in this species compared with A. alba and/or to sampling biases (i.e. more unbalanced sampling) in L. decidua. The overall rate of outlier loci was comparable to that found in other genome scans in boreal conifers (3.7% in Picea glauca (Namroud et al., 2008) and 2.4% and 2.7% for temperature and precipitation clustering, respectively, in Picea mariana (Prunier et al., 2011)), and slightly lower (6% and 4% for positive and negative outliers, respectively) than that found in Pinus pinaster (Eveno et al., 2008). The detection of outlier loci was similar across methods for A. alba, but varied for L. decidua from 2.6% using the island model method of Beaumont & Nichols (1996) to 6.86% using the Bayesian approach implemented in BAYESCAN (Foll & Gaggiotti, 2008), whereas approaches using a more complex demographic model (i.e. the hierarchical island model of Excoffier et al., 2009) produced intermediate numbers (3%). These differences are probably related to population genetic structure in L. decidua that makes the relatively simple model used for data simulation in Beaumont & Nichols (1996) unrealistic (Helyar et al., 2011). Moreover, the method of Excoffier et al. (2009) is known to provide a more conservative estimation of the number of outlier loci, with a higher false-negative rate (Excoffier et al., 2009). In both species, the majority of the outlier loci detected with BAYESCAN were not significant with the other methods, which reflects the complexity of identifying the precise causal loci for environmental adaptation. Nevertheless, our results indictate more widespread signatures of local adaptation in L. decidua, as Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

Research 189

well as confirming the advantage of combining several approaches for the identification of candidate functional loci associated with environment (De Mita et al., 2013). Although gene annotation in nonmodel species is still relatively poor (Ekblom & Galindo, 2011), functional annotation of the Arabidopsis homologues of outlier loci highlighted the presence of several interesting targets for further investigation of subalpine forest tree molecular adaptation. Among the loci detected under the geographical group clustering, the Arabidopsis homologues of locus CL3116Contig1_03 (NM_122008; AT5G20010) found in A. alba encode a protein located in the cell wall, which is produced in response to cadmium ion and salt stress (Meier & Brkljacic, 2010). A member of the same family protein (the RAN GTPase gene family) was found to be associated with phenology and cold tolerance in Pseudotsuga menziesii (Eckert et al., 2009). Another interesting locus detected under the finite island model in L. decidua, CL71Contig1_04 (NM_001036704), encodes a putative disease resistance protein, ADR1-like 1, which is involved in apoptosis and the defence response. Among the outliers detected with the Bayesian simulation, two loci were found in A. alba on a gene (CL4354Contig1_01) encoding a protein serine/threonine phosphatase (PP2A-3; NM_115712) and other two loci were identified in L. decidua on a gene (2_9465_01) encoding a GTP-binding protein (NM_123975). Interestingly, most outlier loci in this study were not among those involved in environmental associations in a previous study (Mosca et al., 2012b), which suggests complex interactions between environment and genetic differentiation, and the existence of other environmental variables (not measured in Mosca et al., 2012b) responsible, at least partially, for patterns of population genetic structure found in subalpine forest trees (e.g. GDD5 and soil type). A similar discordant pattern has been shown in other species (Nosil et al., 2008; Eckert et al., 2010b) but not in P. mariana, where allele frequency was also correlated with temperature and precipitation in 62% of the outlier loci (Prunier et al., 2011). The lack of association between allele frequency and climatic variables may be attributable to a more complex adaptation to local environment, which is only partially explained by the studied variables (Prunier et al., 2011). Finally, two outlier loci in L. decidua were detected in genes significantly associated with a phenotypic trait (autumn cold hardiness or budset timing) in Picea sitchensis: locus 0_6659_01 encoding a peroxidase and locus CL1077Contig1_02 encoding a histone (Holliday et al., 2008, 2010b). Further insights into geographical and environmental (temperature, precipitation, altitude and soil type) effects on genetic structure came from Mantel tests and MRDM analyses. In both species, geographical distance (i.e. significant IBD) accounts for most of the genetic differentiation among populations. A strong correlation of geographical distance and pairwise genetic distance has been found in other widespread forest trees, such as P. sitchensis (Gapare et al., 2005; Mimura & Aitken, 2007), Pinus mugo (Heuertz et al., 2010), and Picea abies (Tollefsrud et al., 2009). Interestingly, Mantel tests showed a significant effect of both temperature and precipitation in A. alba and of precipitation in L. decidua, after correction for the geographical distance. New Phytologist (2014) 201: 180–192 www.newphytologist.com

New Phytologist

190 Research

Alpine vegetation is strongly dependent on climate and plants must be adapted to rapid changes (Rebetez et al., 2004). Several studies found adaptive loci along steep environmental gradients, such as those created by altitude (Gonzalo-Turpin & Hazard, 2009), latitude (Hall et al., 2007; Chen et al., 2012) and precipitation (Tsumura et al., 2012). A significant altitudinal cline for growth in A. alba and other forest species was confirmed using a common garden experiment (Vitasse et al., 2009). A significant effect of soils in the central and eastern Alps was also apparent from MRDMs. To our knowledge, the potential effect of carbonate versus silicate soil type on genetic structure has not been investigated in other forest tree species. The greatest proportion of the forest root system is concentrated in the upper soil horizon (Vogt et al., 1983), where fine roots are often abundant and the effect of bedrock type is minimal. Soil physical and chemical properties are also affected by climate, as was demonstrated along an altitudinal gradient characterized by decreasing temperature and increasing precipitation (Bockheim et al., 2000). Although the present study did not include all factors involved in soil formation and affecting soil properties, it is one of the first to report an association between genetic structure and soil type. This result highlights the importance of further investigations on this aspect of forest ecosystems. Conclusions Subalpine forest trees have large amounts of gene flow compared with herbaceous and annual plants (Savolainen et al., 2007). Gene flow plays an important role in promoting or constraining adaptive variation (Hendry & Taylor, 2004; Garant et al., 2007; Kremer et al., 2012). This relationship is complex, as a geographical barrier can increase population isolation and still, if genetic variation is high, result in local adaptation. Conversely, strong environmental selective pressure can result in local adaptation even in the presence of high gene flow. In the present study, we have found evidence of both IBD and environment-driven adaptation (IBA) in A. alba and L. decidua. We also demonstrated how the effect of environment can be partially separated from other confounding factors, such as geographical distance. Moreover, we identified winter precipitation and early spring temperature (March) as key environmental factors that contribute to the genetic structure of two conifers of the Alps. This study found some potentially adaptive loci based on outlier detection for environmental groups constructed using these variables. However, further investigations are necessary to confirm their involvement in adaptation. For example, the role of these SNPs in physiological processes could be investigated in large-scale association studies. By integrating climatic analysis with a landscape genetic approach, we found that in both species precipitation is involved in creating population genetic structure, whereas in A. alba temperature and soil may also drive genetic differentiation across the landscape, at least in some regions. Our conclusion is that, in subalpine trees that have large amounts of gene flow and little population genetic structure, environment-driven selection is still important and, together with geographical isolation, may promote genetic adaptation. This process is favoured by the large New Phytologist (2014) 201: 180–192 www.newphytologist.com

amount of seeds produced during the life of a tree (i.e. its life-history traits) and the strong levels of selection at early stages of development (i.e. seedlings and saplings).

Acknowledgements The authors are grateful to Erica Di Pierro, Lorenzo Bonosi and Alessio Fortunati for their useful comments during preparation of the manuscript. We would like to thank Piero Belletti, Andrew J. Eckert, Alessandro Mancabelli and Nicola La Porta for their help with the sampling design, and the Italian State Forest Service of Alpine Regions, David Blanco, Yuri Gori, Stefano Maffei, Ambrogio Molinari, Marta Scalfi and Daniele Sebastiani for their support during the sampling. We acknowledge Katie Tsang and Randi Famula for their laboratory work, Ben Figueroa for managing data storage, and Jill L. Wegrzyn, John D. Liechty, Vanessa K. Rashbrook and Charles M. Nicolet for their support with the genotyping. We thank all members of the GIS and Remote Sensing Unit at Fondazione Edmund Mach for providing us with the environmental data and for their help. The ACE-SAP project was partially funded by the Autonomous Province of Trento (Italy), with the regulation No. 23, June 12, 2008, of the University and Scientific Research Service. Thanks are extended to the ERA-Net BiodivERsA (LinkTree project, EUI2008-03713), with the Spanish Ministry of Economy and Competitiveness as national funder, part of the 2008 BiodivERsA call for research proposals, which supported the work of S.C.G.-M.

References Akey JM, Zhang G, Zhang K, Jin L, Shriver MD. 2002. Interrogating a high-density SNP map for signatures of natural selection. Genome Research 12: 1805–1814. Alberto F, Aitken SN, Al
ıa R, Gonz
alez-Mart
ınez SC, H€anninen H, Kremer A, Lefevre F, Lenormand T, Yeaman S, Whetten R et al. 2013. Potential for evolutionary responses to climate change – evidence from tree populations. Global Change Biology 19: 1645–1661. Andrew RL, Ostevik KL, Ebert DP, Rieseberg LH. 2012. Adaptation with gene flow across the landscape in a dune sunflower. Molecular Ecology 21: 2078– 2091. Balding DJ. 2003. Likelihood-based inference for genetic correlation coefficients. Theoretical Population Biology 63: 221–230. Beaumont MA. 2005. Adaptation and speciation: what can FST tell us? Trends in Ecology and Evolution 20: 435–440. Beaumont MA, Nichols RA. 1996. Evaluating loci for use in the genetic analysis of population structure. Proceedings of the Royal Society of London, Series B 263: 1619–1626. Bockheim JG, Munroe JS, Douglass D, Koerner D. 2000. Soil development along an elevational gradient in the southeastern Uinta Mountains, Utah, USA. Catena 39: 169–185. Carrer M, Nola P, Motta R, Urbinati C. 2010. Contrasting tree-ring growth to climate responses of Abies alba toward the southern limit of its distribution area. Oikos 119: 1515–1525. Chen J, K€a llman T, Ma X, Gyllenstrand N, Zaina G, Morgante M, Bousquet J, Eckert A, Wegrzyn J, Neale D et al. 2012. Disentangling the roles of history and local selection in shaping clinal variation in allele frequency and gene expression in Norway spruce (Picea abies). Genetics 191: 865–881. Coyer JA, Hoarau G, Pearson G, Mota C, J€ uterbock A, Alpermann T, John U, Olsen JL. 2011. Genomic scans detect signatures of selection along a salinity gradient in populations of the intertidal seaweed Fucus serratus on a 12 km scale. Marine Genomics 4: 41–49. Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

New Phytologist De Carvalho D, Ingvarsson PK, Joseph J, Suter L, Sedivy C, Macaya-Sanz D, Cottrell J, Heinze B, Schanzer I, Lexer C. 2010. Admixture facilitates adaptation from standing variation in the European aspen (Populus tremula L.), a widespread forest tree. Molecular Ecology 19: 1638–1650. De Mita S, Thuillet A-C, Gay L, Ahmadi N, Manel S, Ronfort J, Vigouroux Y. 2013. Detecting selection along environmental gradients: analysis of eight methods and their effectiveness for outbreeding and selfing populations. Molecular Ecology 22: 1383–1399. Eckert AJ, Bower AD, Gonz
a lez-Martínez SC, Wegrzyn JL, Coop G, Neale DB. 2010b. Back to nature: ecological genomics of loblolly pine (Pinus taeda, Pinaceae). Molecular Ecology 19: 3789–3805. Eckert AJ, Bower AD, Pande B, Jermstad KD, Krutovsky KV, St. Clair JB, Neale DB. 2009. Association genetics of coastal Douglas fir (Pseudotsuga menziesii var. menziesii, Pinaceae). I. Cold-hardiness related traits. Genetics 182: 1289–1302. Eckert AJ, van Heerwaarden J, Wegrzyn JL, Nelson CD, Ross-Ibarra J, Gonz
a lez-Mart
ınez SC, Neale DB. 2010a. Patterns of population structure and environmental associations to aridity across the range of loblolly pine (Pinus taeda L., Pinaceae). Genetics 185: 969–982. Ekblom R, Galindo J. 2011. Applications of next generation sequencing in molecular ecology of non-model organisms. Heredity 107: 1–15. European Communities. 2001. Georeferenced Soil Database for Europe. Manual of procedures Ver 1.1. Eveno E, Collada C, Guevara MA, Leger V, Soto A, D
ıaz L, Le´ger P, Gonz
a lez-Mart
ınez SC, Cervera MT, Plomion C et al. 2008. Contrasting patterns of selection at Pinus pinaster Ait. drought stress candidate genes as revealed by genetic differentiation analyses. Molecular Biology and Evolution 25: 417–437. Excoffier L, Hofer T, Foll M. 2009. Detecting loci under selection in a hierarchically structured population. Heredity 103: 285–298. Excoffier L, Lischer H. 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources 10: 564–567. Excoffier L, Smouse P, Quattro J. 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131: 479–491. Farjon A. 1990. Pinaceae: drawings and descriptions of the genera Abies, Cedrus, Pseudolarix, Keteleeria, Nothotsuga, Tsuga, Cathaya, Pseudotsuga, Larix and Picea. K€onigstein, Germany: Koeltz Scientific Books. Foll M, Gaggiotti OE. 2008. A genome scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics 180: 977–993. Gapare WJ, Aitken SN, Ritland CE. 2005. Genetic diversity of core and peripheral Sitka spruce (Picea sitchensis (Bong.) Carr) populations: implications for conservation of widespread species. Biological Conservation 123: 113–123. Garant D, Forde SE, Hendry AP. 2007. The multifarious effects of dispersal and gene flow on contemporary adaptation. Functional Ecology 21: 434–443. Gonzalo-Turpin H, Hazard L. 2009. Local adaptation occurs along altitudinal gradient despite the existence of gene flow in the alpine plant species Festuca eskia. Journal of Ecology 97: 742–751. Goslee SC, Urban DL. 2007. The ecodist package for dissimilarity-based analysis of ecological data. Journal of Statistical Software 22: 1–19. Goudet J. 2005. HIERSTAT, a package for R to compute and test hierarchical F-statistics. Molecular Ecology Notes 5: 184–186. Hall D, Luquez V, Garcia VM, St Onge KR, Jansson S, Ingvarsson PK. 2007. Adaptive population differentiation in phenology across a latitudinal gradient in European Aspen (Populus tremula, L.): a comparison of neutral markers, candidate genes and phenotypic traits. Evolution 61: 2849–2860. Hardy OJ, Vekemans X, Cartwright RA. 2009. SPAGeDi 1.3: a program for spatial pattern analysis of genetic diversity [WWW document]. Last modified on 22 March 2009 [accessed 24 April 2012]. Helyar SJ, Hemmer-Hansen J, Bekkevold D, Taylor MI, Ogden R, Limborg MT, Cariani A, Maes GE, Diopere E, Carvalho GR et al. 2011. Application of SNPs for population genetics of nonmodel organisms: new opportunities and challenges. Molecular Ecology Resources 11: 123–136. Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust

Research 191 Hendry AP, Taylor EB. 2004. How much of the variation in adaptive divergence can be explained by gene flow? An evaluation using lakestream stickleback pairs. Evolution 58: 2319–2331. Heuertz M, Teufel J, Gonz
a lez-Mart
ınez SC, Soto A, Fady B, Al
ıa R, Vendramin GG. 2010. Geography determines genetic relationships between species of mountain pine (Pinus mugo complex) in western Europe. Journal of Biogeography 37: 541–556. Holliday JA, Ralph SG, White R, Bohlmann J, Aitken SN. 2008. Global monitoring of autumn gene expression within and among phenotypically divergent populations of Sitka spruce (Picea sitchensis). New Phytologist 178: 103–122. Holliday JA, Ritland K, Aitken SN. 2010b. Widespread, ecologically relevant genetic markers developed from association mapping of climate-related traits in Sitka spruce (Picea sitchensis). New Phytologist 188: 501–514. Holliday JA, Yuen M, Ritland K, Aitken SN. 2010a. Postglacial history of a widespread conifer produces inverse clines in selective neutrality tests. Molecular Ecology 19: 3857–3864. Husson F, Josse J, Le S, Mazet J. 2012. FactoMineR: Multivariate Exploratory Data Analysis and Data Mining with R. version 1.20. 2012-10-02 06:51:50. Joost S, Bonin A, Bruford MW, Depr
es L, Conord C, Erhardt G, Taberlet P. 2007. A spatial analysis method (SAM) to detect candidate loci for selection: towards a landscape genomics approach to adaptation. Molecular Ecology 16: 3955–3969. Keller SR, Levsen N, Olson MS, Tiffin P. 2012. Local adaptation in the flowering-time gene network of balsam poplar, Populus balsamifera L. Molecular Biology Evolution 29: 3143–3152. Keller SR, Soolanayakanahally RY, Guy RD, Silim SN, Olson MS, Tiffin P. 2011. Climate-driven local adaptation of ecophysiology and phenology in balsam poplar, Populus balsamifera L. (Salicaceae). American Journal of Botany 98: 99–108. K€orner C. 2003. Alpine plant life – functional plant ecology of high mountain ecosystems, 2nd edn. Heidelberg, Germany: Springer. ISBN 3-540-00347-9. Koz
a kov
a R, Samonil P, Nov
ak J, Kunes P, Koc
a r P, Koc
arov
a R. 2011. Contrasting local and regional Holocene histories of Abies alba in the Czech Republic in relation to human impact: evidence in forestry, pollen and anthracological data. Holocene 21: 431–444. Kremer A, Le Corre V, Petit RJ, Ducousso A. 2010. Historical and contemporary dynamics of adaptive differentiation in European oaks. In: DeWoody A, Bickham J, Michler C, Nichols K, Rhodes G, Woeste K, eds. Molecular approaches in natural resource conservation. Cambridge, UK: Cambridge University Press, 101–117. Kremer A, Ronce O, Robledo-Arnuncio JJ, Guillaume F, Bohrer G, Nathan R, Bridle JR, Gomulkiewicz R, Klein EK et al. 2012. Long-distance gene flow and adaptation of forest trees to rapid climate change. Ecology Letters 15: 378–392. Lee CR, Mitchell-Olds T. 2011. Quantifying effects of environmental and geographical factors on patterns of genetic differentiation. Molecular Ecology 20: 4631–4642. Leimu R, Fischer M. 2008. A meta-analysis of local adaptation in plants. PLoS One 3: e4010. Lewontin R, Krakauer J. 1973. Distribution of gene frequency as a test of theory of selective neutrality of polymorphisms. Genetics 74: 175–195. Liepelt S, Cheddadi R, Debeaulieu J, Fady B, G€om€ori D et al. 2009. Postglacial range expansion and its genetic imprints in Abies alba (Mill.) – a synthesis from palaeobotanic and genetic data. Review of Palaeobotany and Palynology 153: 139–149. McCune B, Dylan K. 2002. Equations for potential annual direct incident radiation and heat load. Journal of Vegetation Science 13: 603–606. Meier I, Brkljacic J. 2010. The Arabidopsis nuclear pore and nuclear envelope. The Arabidopsis Book 8: e0139. Mimura M, Aitken SN. 2007. Adaptive gradients and isolation-by-distance with postglacial migration in Picea sitchensis. Heredity 99: 224–232. Mimura M, Aitken SN. 2010. Local adaptation at the range peripheries of Sitka spruce. Journal Evolutionary Biology 23: 249–258. Mosca E, Eckert AJ, Di Pierro EA, Rocchini D, La Porta N, Belletti P, Neale DB. 2012b. The geographical and environmental determinants of genetic diversity for four alpine conifers of the European Alps. Molecular Ecology 21: 5530–5545. New Phytologist (2014) 201: 180–192 www.newphytologist.com

New Phytologist

192 Research Mosca E, Eckert AJ, Liechty JD, Wegrzyn JL, La Porta N, Vendramin GG, Neale DB. 2012a. Contrasting patterns of nucleotide diversity for four conifers of Alpine European forests. Evolutionary Applications 5: 762–775. M€ uller-Starck G, Baradat P, Bergmann F. 1992. Genetic variation within European tree species. New Forests 6: 23–47. Namroud MC, Beaulieu J, Juge N, Laroche J, Bousquet J. 2008. Scanning the genome for gene single nucleotide polymorphisms involved in adaptive population differentiation in white spruce. Molecular Ecology 17: 3599–3613. Nosil P, Egan SP, Funk DJ. 2008. Heterogeneous genomic differentiation between walking-stick ecotypes: “isolation by adaptation” and multiple roles for divergent selection. Evolution 62: 316–336. Panagos P, Jones A, Bosco C, Senthil Kumar PS. 2011. European digital archive on soil maps (EuDASM): preserving important soil data for public free access. International Journal of Digital Earth 4: 434–443. Pavy N, Namroud M-C, Gagnon F, Isabel N, Bousquet J. 2012. The heterogeneous levels of linkage disequilibrium in white spruce genes and comparative analysis with other conifers. Heredity 108: 273–284. Petit RJ, Hu FS, Dick CK. 2008. Forests of the past: a window to future changes. Science 320: 1450–1452. Pleuss AR. 2011. Pursuing glacier retreat: genetic structure of a rapidly expanding Larix decidua population. Molecular Ecology 20: 473–485. Prentice IC, Cramer W, Harrison SP, Leemans R, Monserud RA, Solomon AM. 1992. A global biome model based on plant physiology and dominance, soil properties and climate. Journal of Biogeography 19: 117–134. Prentice IC, Harrison SP, Bartlein PJ. 2011. Global vegetation and terrestrial carbon cycle changes after the last ice age. New Phytologist 189: 988–998. Prunier J, Laroche J, Beaulier J, Bousquet J. 2011. Scanning the genome for gene SNPs related to climate adaptation and estimating selection at the molecular level in boreal black spruce. Molecular Ecology 20: 1702–1716. Quantum GIS Development Team. 2011. Quantum GIS Geographic Information System. Open Source Geospatial Foundation Project. R Development Core Team. 2011. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. [WWW document] URL http://www.R-project.org [accessed 8 July 2011]. Rebetez M, Reinhard M, Buttler A. 2004. Forests, tree physiology and climate. Encyclopedia of forest sciences. London, UK: Academic Press. Rousset F. 2004. Genetic structure and selection in subdivided populations. Princeton, NJ, USA: Princeton University Press. Sartori G, Mancabelli A. 2009. Carta dei suoli del Trentino: scala 1:250.000. Museo Tridentino di Scienze Naturali. Savolainen O, Pyh€a j€arvi T, Kn€ urr T. 2007. Gene flow and local adaptation in trees. Annual Review of Ecology and Systematics 38: 595–619. Schl€otterer C, Harr B. 2002. Single nucleotide polymorphisms derived from ancestral populations show no evidence for biased diversity estimates in Drosophila melanogaster. Molecular Ecology 11: 947–950. Schlotterer C. 2002. A microsatellite-based multilocus screen for the identification of local selective sweeps. Genetics 160: 753–763. Semerikov VL, Lascoux M. 1999. Genetic relationship among Eurasian and American Larix species based on allozymes. Heredity 83: 62–70. Slatkin M, Voelm L. 1991. F(ST) in a hierarchical island model. Genetics 127: 627–629. Sork VL, Davis FW, Westfall R, Flint A, Ikegami M, Wang H, Grivet D. 2010. Gene movement and genetic association with regional climate gradients in California valley oak (Quercus lobata N
ee) in the face of climate change. Molecular Ecology 19: 3806–3823. Theurillat J-P, Guisan A. 2001. Potential impact of climate change on vegetation in the European Alps: a review. Climate Change 50: 77–109. Tollefsrud MM, Sonstebo JH, Brochmann C, Johnsen O, Skroppa T, Vendramin GG. 2009. Combined analysis of nuclear and mitochondrial markers provide new insight into the genetic structure of North European (Picea abies). Heredity 102: 549–562. Tsumura Y, Uchiyama K, Moriguchi Y, Ueno S, Ihara-Ujino T. 2012. Genome scanning for detecting adaptive genes along environmental gradients in the Japanese conifer, Cryptomeria japonica. Heredity 109: 349–360. Vitasse Y, Delzon S, Bresson CC, Michalet R, Kremer A. 2009. Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden. Canadian Journal of Forest Research 39: 1259–1269. New Phytologist (2014) 201: 180–192 www.newphytologist.com

Vogt KA, Grier CG, Meier CE, Keyes MR. 1983. Organic matter and nutrient dynamics in forest floors of young and mature Abies amabilis stands in Western Washington, as affected by fine root input. Ecological Monographs 53: 139–157. Weir BS, Cockerham CC. 1984. Estimating F-statistics for the analysis of population structure. Evolution 38: 1358–1370.

Supporting Information Additional supporting information may be found in the online version of this article. Fig. S1 Plot of the nine environmental groups in Abies alba and Larix decidua. Fig. S2 LD plots based on 231 and 233 SNP markers in Abies alba and Larix decidua, respectively. Fig. S3 Detection of outlier loci in Abies alba and in Larix decidua using the Bayesian approach. Table S1 List of the sampling sites for Abies alba and Larix decidua assigned to the geographical groups (Geo-Group, G1:G4) and to environmental groups (Env-Group) Table S2 All variable contributions (AFDM analyses) to the individual distribution in Abies alba and Larix decidua; variables that mainly contributed to the first two axes in A. alba and L. decidua Table S3 Expected (He) and observed (Ho) heterozygosities calculated in each sampling site (population) for Abies alba and Larix decidua Table S4 Outlier loci found in Abies alba and Larix decidua with the island model and the hierarchical island model with two regions or three/four geographical groups and with nine environmental groups Table S5 Functional annotation of all outlier loci detected with the island model (M1) and the hierarchical island model with two regions (M2) or three/four geographical groups (M3), with nine environmental groups (M4) and the Bayesian approach (M5), in Abies alba and Larix decidua Table S6 Multiple regressions of distance matrices (MRMs), estimating the association of genetic distance (F) considering all populations and only populations within geographical groups with elevation (E), annual mean temperature (T), annual cumulative precipitation (P) and soil type (S) Notes S1 All studied genes in Abies alba (ABAL) and Larix decidua (LADE) with their putative function. Please note: Wiley Blackwell are not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office. Ó 2013 The Authors New Phytologist Ó 2013 New Phytologist Trust