Post-glacial re-colonization and natural selection have shaped growth responses of silver fir across Europe

• Long-term growth and its climate responses have a clear genetic basis. • Natural selection and past demography explain intra-specific differences in


Introduction
Forest managers in Europe are today facing an enormous challenge with ongoing climate change to maintain sustainable forest ecosystem services and healthy forests for the next decades. In particular, the steady increase in temperature as well as the increase in frequency and magnitude of extreme drought events and heat waves may threaten forest species as reflected in the recent massive dieback occurring worldwide (Allen et al., , 2015. The two recent extreme summer droughts of 2003 and 2018 that occurred in Central and western Europe have largely impacted tree growth and vitality (Vitasse et al., 2019a;Brun et al., 2020;Buras et al., 2020;Schuldt et al., 2020), and may reflect what will likely become normality in the near future. Young trees that are naturally establishing today will face increasing temperature and more frequent extreme droughts. Phenotypic plasticity, migration and/or genetic adaptation will allow trees to partly cope with these novel conditions in-situ, but the adaptive capacity of trees may be limited due to the rapidity of ongoing climatic changes (Aitken et al., 2008;Dauphin et al., 2021). Because forest trees are long-lived and slow to reach reproductive maturity, studies assessing shifts in species distribution of forest trees and accounting for demographic processes and competition suggest that they will mostly lag behind the shift of their climatic niche (Fréjaville et al., 2020;Scherrer et al., 2020).
In this context, forest managers rightly question whether to promote natural regeneration or to apply assisted migration/gene flow using seeds or seedlings of alternative species/provenances that are adapted to warmer and drier conditions (Aitken and Whitlock, 2013;Aitken and Bemmels, 2016). For a given species, non-local populations might have a higher potential to tolerate drought than the local populations due to their history of migration since the last glaciation (re-colonization history) and/or due to the adaptation to their contemporary, local climate (natural selection). However, it remains to be tested to what extent populations from these genetically differentiated lineages are able to tolerate future extreme droughts and warmer climate.
In Europe, silver fir (Abies alba Mill.) is a very valuable tree species providing important economic, ecological and cultural services. It is an evergreen conifer that grows in a wide variety of soil types with different nutrient contents and pH values, forming pure or mixed stands with other temperate tree species. Silver fir is a shade-tolerant (Dobrowolska et al., 2017) and sensitive species to vapor pressure deficits and atmospheric drought (Aussenac, 2002). In Central Europe, after a period of growth decline due to air pollution prior to the 1980s (see Elling et al., 2009), growth of silver fir has substantially increased, particularly after 2000 (Büntgen et al., 2014). In contrast, growth decline has been reported since the 1980s towards its southern distribution limit (e.g., the Pyrenees region) due to the increasing occurrence of dry summers (Macias et al., 2006;Peguero-Pina et al., 2007;Camarero et al., 2011;Linares and Camarero, 2012;Büntgen et al., 2014;Gazol et al., 2015;Davi and Cailleret, 2017;Latreille et al., 2017). More recently, triggered by the extreme 2018 summer drought, growth decline symptoms have been observed in Central Europe and also affecting many other tree species (Schuldt et al., 2020).
The current distribution area of silver fir is limited mainly to the mountainous regions of eastern, western, southern and Central Europe. As for many European tree species (Petit et al., 2003), its distribution reflects both (i) the history of its migration after the last glaciation (Lang, 1994;Liepelt et al., 2009), and (ii) the human influence since the Neolithic age that has largely contributed to fragment and reduce its distribution (Tinner et al., 1999;Vitasse et al., 2019a;Martínez-Sancho and Gutiérrez, 2019). Palaeobotanical and molecular genetic data suggested that there were four main refugia during the last glaciation period: the Calabrian region in southern Italy, the Apennine mountains in northwestern Italy, the Balkan peninsula in northwestern Greece, and the Pyrenees in northeastern Spain and southwestern France (Konnert and Bergmann, 1995;Muller et al., 2007;Cheddadi et al., 2014). Populations from the Calabrian region and the Pyrenees do not seem to have contributed to the recolonization of silver fir during the Holocene in Central Europe, but rather the populations originating from the refugia of northwestern Italy and the southern Balkans (Muller et al., 2007;Cheddadi et al., 2014). In contrast, the populations from the Pyrenees remained isolated from the main distribution range and may have been subject to severe bottlenecks resulting in reduced genetic diversity (Vendramin et al., 1999), whereas recent studies suggest a potential trans-Adriatic gene flow via pollen between the populations from the Calabrian region and the Balkans (Piotti et al., 2017). Thus, the post-glacial recolonization history of silver fir populations has strongly shaped the genetic diversity of the species and created several separate lineages with likely distinct responses to environmental fluctuations (Bosela et al., 2016). Thus, to assess the full potential of a species to thrive in a given region in future warmer and drier conditions, it is necessary to examine the response of genetically differentiated populations from different lineages to climatic conditions using fitness-related traits such as growth or resilience to extreme drought events.
Growing multiple seed sources under the same environment (so called 'provenance trials' or 'common garden experiments') reveals genetically based population differences in tree fitness (e.g., Vitasse et al., 2009;Housset et al., 2018). Provenance trials are of high value but still under-exploited to investigate intra-specific variability of growth-related traits for forest trees as it requires i) important resources for installation and maintenance of trials over the years; ii) long-term observations that ideally includes trees' adult stages; and iii) enough local climatic variability, including extreme events, to compare responses among provenances. Consequently, most of the studies conducted in common gardens have been restricted to a specific life stage (e.g., juvenile stage), and/or a short time period (e.g., one to two growing seasons). In this context, dendrochronology is an excellent tool to overcome some of these issues because i) tree-ring data allow for the retrospective analysis of growth responses to climate at tree, population, and forest ecosystem levels and ii) growth integrates most physiological processes related to resource acquisition and use (Fritts, 1976). Although some studies have started using dendrochronological techniques in provenance trials located in North America (Isaac-Renton et al., 2018;Sang et al., 2019;Depardieu et al., 2020), this approach still remains largely unexplored for most of the European tree species (but see Trujillo-Moya et al., 2018). Moreover, there are only a few studies that combine dendrochronological and genetic analyses to assess the roles of demographic processes and local adaptation in long-and short-term growth performance (e.g., Depardieu et al., 2020). Identifying populations that might be better adapted to the future climatic conditions of a given site based on long-term growth and extreme-event response patterns is, however, crucial for optimizing assisted gene flow strategies.
Here we evaluate the ability of 18 provenances of silver fir from different regions across Europe to tolerate increasing temperature and extreme droughts using trees of two common gardens established at the end of the 1980s in Switzerland. Based on tree-ring width series and genetic data, we investigated the impact of climatic conditions on various growth traits including the resistance, recovery, and resilience of the provenances to the exceptional 2003 summer drought. Specifically, we addressed the following questions: (i) Do growth responses to climate including resistance, recovery, and resilience to the exceptional 2003 summer drought differ among provenances? (ii) Are differences associated with the recolonization history of the provenances after the last glacial maximum? (iii) Is there evidence that these differences are a result of natural selection?

Experimental design and field measurements
Plant material (increment cores for the tree-ring series and needles for the genetic samples) were collected in summer 2018 from two distinct common garden experiments located in Switzerland: Bourrignon (47°23′N, 7°15′E, 800 m a.s.l.) and Lenzburg (47°24′N 8°9′E, 400 m a.s.l.). Both common gardens were established during the late 1980s and included silver fir trees of a total of 18 provenances from across the species' natural distribution range (Fig. 1, Table 1). The climate at the location of the common gardens is moderate continental (Dfb, Köppen classification) but Lenzburg (total annual precipitation of 978 mm and mean annual temperature of 10.9°C) is substantially drier and warmer than Bourrignon (total annual precipitation of 1439 mm and mean annual temperature 8.5°C). Precipitation is generally evenly distributed throughout the year at both locations ( Fig. S1), but high values peak in autumn. The experimental design is a randomized block design with three blocks in Bourrignon and two in Lenzburg.
Up to five individuals per provenance and block were sampled, resulting in a total of 253 trees across both common gardens (Table 1).
Tree height and diameter at breast height (DBH) were measured in the field. In the case of the site at the higher elevation (Bourrignon), one increment core per tree was extracted at breast height (1.3 m) using an increment borer, perpendicular to the slope direction to avoid reaction wood. In the case of Lenzburg, stem discs were cut at breast height for some individuals. Needle tissue for DNA extraction was stored at −20°C until further processing. Cores and stems were air dried and then sanded with progressively finer sanding paper until wood cells were clearly visible under a binocular microscope. Tree-ring widths were measured with a precision of 0.01 mm using a binocular microscope connected to a LINTAB measuring table (Rinntech, Heidelberg, Germany). Two radii were measured and averaged for each stem disk.

Climatic data
Time series of monthly maximum and minimum temperature as well as total precipitation from the two common gardens were extracted from the CHELSA high-resolution climate dataset (Karger et al., 2017) for the period 1993-2016. Similarly, we also extracted time series for all provenances (climate of seed origin) from the same dataset for  Mauri et al. (2017), and the corresponding climate was extracted from CHELSA (Karger et al., 2017). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 1 Geographic origin of the analyzed provenances of silver fir. B and L indicate the common garden where the provenance was planted (Bourrignon and Lenzburg, respectively). the climate normal 1961-1990. A total of 19 bioclimatic parameters were then calculated for each provenance using the R package dismo (Hijmans et al., 2011): annual mean temperature, mean diurnal range, isothermality, temperature seasonality, max temperature of warmest month, min temperature of coldest month, temperature annual range, mean temperature of wettest quarter, mean temperature of driest quarter, mean temperature of warmest quarter, mean temperature of coldest quarter, annual precipitation, precipitation of wettest month, precipitation of driest month, precipitation seasonality (coefficient of variation), precipitation of wettest quarter, precipitation of driest quarter, precipitation of warmest quarter and precipitation of coldest quarter.

Tree-ring analyses
The age-related trends from each individual tree-ring width series were removed by fitting a cubic smoothing spline with a 50% frequency cut-off at 20 years. The detrended series were then averaged at common garden level by using a Tukey's biweight robust mean (Fig. S2). To assess the main climatic drivers of tree growth in each common garden, bootstrapped Pearson correlations were calculated between the two chronologies and monthly mean temperature and precipitation of the growing season (March to September) for the period 1993-2016 using the R package treeclim (Zang and Biondi, 2015). Only significant correlations at common garden level were tested for differences among provenances and genetic clusters. Correlations between detrended series and climate conditions were also calculated for all the individuals. We also calculated individual resilience components (Lloret et al., 2011) to the 2003 extreme drought (Fig. S3). The resistance (Rt) quantifies the ratio between growth during the drought period and growth during the preceding period, thus representing the capacity of the trees to buffer drought stress. The recovery (Rc) quantifies the growth reaction following the drought period and is defined by the ratio between growth during the post-drought period and growth during the drought period. The resilience (Rs) quantifies the ratio between growth during the post-drought period and growth during the pre-drought period, which represents the capacity of trees to regain the growth of the pre-drought period. We delimited preand post-drought periods to two years to ensure that the lag effect of a given extreme climatic event does not overlap with another extreme event (Anderegg et al., 2015). The levels of first-order autocorrelation of the series (Ac) were also calculated since this parameter can be associated with the "biological memory" of trees (Speer, 2010). This parameter can be used to assess the level of sensitivity to climate fluctuations since high levels of autocorrelation are usually linked to complacent growth.

DNA extraction and SNP genotyping
We lyophilized 16-26 mg of needle tissue per tree for 18 h and subsequently grinded the samples three times for 1 min using metal beads in a Retsch mill (Retsch, Haan, Germany). DNA was then extracted with a KingFisher 96 system (Thermo Fisher Scientific, Waltham, MA, USA) using the sbeadex maxi plant kit originally optimized for oak tree leaves (LGC Genomics, Berlin, Germany), adding 1% Thioglycerol and 0.8% RNAse to the lysis buffer. DNA concentrations were measured with Quantus (Promega Corporation, Madison, WI, USA) and Picogreen (Thermo Fisher Scientific) on an Infinite M200 microplate reader (Tecan, Männedorf, Switzerland). DNA quality was assessed with a NanoDrop spectrometer (Thermo Fisher Scientific).
We used the Kompetitive Allele Specific PCR (KASP) technology to genotype all individuals in 150 putatively neutral single nucleotide polymorphisms (SNPs) from Csilléry et al. (2020). Genotyping and subsequent bi-allelic SNP scoring was performed by LGC Genomics (Hoddesdon, UK). Four SNPs and three individuals were discarded from the dataset due to a too high proportion of missing data (>10%), resulting in 146 SNPs and 250 genotypes for further genetic analyses.

Population genetic analyses
For each provenance, we pooled genotype data from both common gardens, resulting in 11-24 trees per provenance. We then calculated expected heterozygosity (He) using the R-package adegenet (Jombart, 2008) and the inbreeding coefficient (F IS ) using the R-package hierfstat (Goudet, 2005) for each provenance. Genetic differentiation (Nei's F ST ) among provenances was also assessed with hierfstat. For all F IS and F ST values, we computed 95% confidence intervals with 1000 bootstraps over loci to test whether values are significantly higher than 0 (i.e. lower limit of confidence interval is above 0). Thereafter, we tested for patterns of isolation by distance (IBD) using linearized genetic (F ST /1-F ST ) and transformed geographic (ln) distances (Rousset, 1997) with a Mantel test using the R-package ecodist (Goslee and Urban, 2007) with 1001 permutations. Finally, we used Structure 2.3.4 (Pritchard et al., 2000) to describe neutral population structure among provenances using ten independent runs with different seeds for each number of genetic clusters (K) from 1 to 10, with 1,000,000 repetitions after a burn-in period of 100,000 runs, admixture model, correlated allele frequencies, and no prior population information. Summary statistics for each K were then calculated with Structure Harvester (Earl and vonHoldt, 2012) and average assignment probabilities were computed with Clumpp 1.1.2 (Jakobsson and Rosenberg, 2007) for K = 1-7 using the Greedy algorithm. Because Structure is sensitive to differences in samples size, we cross-checked its results with a principal component analysis (PCA) on genetic data with adegenet, replacing missing values with the average value of the whole dataset.

Statistical analyses
We selected a total of seven traits from three distinct trait types: tree growth performance including tree height, DBH; growth memory, namely Ac; growth climate responses including responses to the 2003 extreme event (Rs, Rc and Rt) and long-term climate responses namely the correlation between tree growth and May temperature. The latter trait was selected since it was the only climate-growth correlation significant in both common gardens. Linear models were used to test for differences in the selected traits among provenances. Provenance as well as block nested in site were used as fixed effects. To disentangle the differences among the identified genetic clusters, linear mixedeffects models (Zuur et al., 2009) were applied using the R package nlme (Pinheiro et al., 2013). Cluster and block were treated as fixed effects, and provenance was included as a random effect on the intercept. Diagnostic plots were checked for heteroscedasticity and non-normal distribution of the residuals, and log-transformations were applied whenever needed. Least square means of provenances and clusters were extracted from the respective models and post-hoc pairwise comparisons among provenances and clusters were performed with a Tukey adjustment using the R package emmeans (Lenth, 2016).
To explore potential signs of adaptation to original climate conditions, Pearson correlations between the least square means of provenances of the traits and bioclimatic variables from the seed origin (1961( -1990 . S4) were performed adjusted for Holm false discovery rate for reporting significances. To check whether trait differences among provenances are a result of natural selection, we also used a Q ST -F ST comparison. While F ST depicts the neutral genetic differentiation among provenances, which can be considered as a null expectation of population differentiation solely due to drift and migration, Q ST is a quantitative analogue to F ST that measures the amount of genetic variation among provenances in relation to the total genetic variance in a trait (reviewed in e.g., Leinonen et al., 2013). If selection is directional, then trait differentiation exceeds neutral expectations, i.e.; Q ST > F ST . In this case, larger genetic differentiation of the phenotypic traits among populations can be interpreted as a consequence of natural selection (Whitlock, 2008). We conducted Q ST -F ST comparisons for those traits that exhibited a significant provenance or cluster effect in the analyses described above, using only the 246 samples for which both trait and genetic data were available. Because the pedigree of the investigated trees is unknown, we estimated genetic relatedness among individual trees after Wang (2002) with the R-package related (Pew et al., 2015) using the 146 SNPs described above. We estimated the additive genetic variance (V A ) of each trait with a linear mixedeffect model (i.e. animal model) implemented in ASReml-R v4 (Butler et al., 2017) with site and block as fixed effects, and the provenance and genetic relatedness matrix as random effects. The animal model is commonly used to partition phenotypic variance into genetic and environmental components (Wilson et al., 2010). We used provenance, and not cluster (as in the analysis above) as sample unit, because the genetic clustering done by Structure is essentially also performed in the animal model. We used the genetic relatedness matrix to estimate the additive genetic variance (V A ) of the phenotypic traits. The between-provenance genetic variance (V B ) was estimated by the variance of the provenance random effect in the animal model. Q ST was then estimated as V B / (V B + 2 × V A ). From V A , we could also estimate the narrow-sense heritability (h 2 ) of the traits. Standard errors were obtained from the model using the varcomp function in R. From these standard errors, confidence intervals were calculated for each trait and compared to the confidence intervals of pairwise F ST . If Q ST > F ST and confidence intervals did not overlap, we considered the trait under divergent selection. This approach is a rough evaluation of the significance of the Q ST -F ST contrast. We could not use the more adequate bootstrap approach of Whitlock and Guillaume (2009) because of the lack of pedigree information. Finally, we performed Pearson correlations to test whether Q ST -values and absolute correlation coefficients of trait values with the climate of seed origin are correlated.
All statistical analyses were performed in the R environment (version 3.5.2, R Development Core Team, 2015) and the graphing was implemented using the R package ggplot2 (Wickham, 2009).

Genetic diversity and structure
Expected heterozygosity (H e ) was on average 0.28 (range = 0.20-0.31, Table S1), with highest values found in central provenances (Switzerland) and the lowest values in peripheral provenances (Ca, La, Py, Ri, St). A significant level of inbreeding was only present in three provenances (Ca, St, To, Table S1). Average pairwise genetic differentiation (F ST ) among provenances was 0.086 (range = −0.001-0.224, Table S2), with provenances Py and Ri showing the highest F ST values. All pairwise F ST values were significantly higher than zero. There was a strong and significant isolation by distance (Mantel r = 0.64, p = 0.001): provenances that come from locations close to each other tended to be genetically similar (Fig. S5).
In the Structure analysis, the mean estimated probability of the data reached a plateau at K = 7 (Fig. S6), indicating that our provenance set consists of seven genetic clusters (Pritchard et al., 2000). However, in K higher than three, additional clusters did not result in individuals being completely assigned to these new clusters (Fig. S7). Rather, these additional clusters illustrate further substructure and the high genetic admixture in Central European provenances, with the exception of the provenance Le which was assigned to a new cluster with K ≥ 5. We therefore concentrated on the three main clusters (K = 3) in further analyses, assigning provenances completely to the cluster with the highest average membership probability (Figs. 2, S7). This was not possible for the highly admixed provenance Si which was therefore Fig. 2. Spatial distribution of neutral genetic population structure. Simplified geographical distribution of the three main genetic clusters of the study system (top). Dark grey areas correspond to the natural silver fir distribution obtained from EUFORGEN (www.euforgen.org). Due to admixture, provenance Si was not assigned to a cluster and is coloured in black. Detailed assignment probabilities to the three main clusters on which the maps are based on (bottom). Each individual tree (trials pooled) is represented by a thin vertical line, which is partitioned into three coloured segments that represent the assignment probabilities of the individuals to the three main clusters determined by Structure (for information on other K values, see Fig. S7). Provenances are ordered longitudinally. excluded from all cluster-based analyses. These three clusters showed a clear longitudinal pattern, with provenance Py as sole member of the westernmost cluster, the Central European provenances representing the central cluster and all the eastern European provenances forming the eastern cluster. The Structure results were confirmed by the PCA, where PC1 separated the eastern and western provenances (with provenance Si in the middle), and PC2 provenance Py from the rest (Fig. S8).

Climate-growth relationships
Climate-growth correlations were performed at common garden level to depict the main climatic factors governing tree growth at both sites, and to identify common drivers between both common gardens (Fig. 3). Three significant correlations were found in the common garden located in the western part of Switzerland (Bourrignon). At this site, the tree-ring width chronology was correlated positively with May temperature and July precipitation, and negatively with July temperature. In the common garden located at lower elevation in the northeast of Switzerland (Lenzburg), the tree-ring width chronology was only significantly and negatively correlated with May temperature. High variability among provenances was observed in these climategrowth correlations.

Differences in growth performance and climatic responses among provenances and genetic clusters
Provenance and common garden had a significant effect (p < 0.05) on traits related to tree growth performance and growth memory such as DBH, tree height and Ac (Table S3). Significant pairwise differences in DBH and height among provenances were observed: Py and Ca together with La stood out as the provenances with the lowest and highest DBH and height values, respectively ( Fig. 4A and B, left). Similarly, cluster had a significant effect when assessing differences in traits among clusters, except for Ac (Table S4). Overall, trees from the eastern cluster tended to have higher DBH and tree height values than trees from the Pyrenees and the central cluster. Consistently, the eastern cluster showed the highest values of DBH and height, significantly higher than the ones from the Pyrenees and the central cluster ( Fig. 4A and B , right). Trees from the Pyrenees clearly showed the lowest long-term performance in terms of DBH and tree height. Levels of autocorrelation were equal among pairwise comparisons of provenances (Fig. 4C left), but trees from Py showed significantly lower levels of autocorrelation than trees from the central cluster, and lower (but not significant) than trees from the eastern cluster (Fig. 4C right).
Provenance had a significant effect on recovery and resistance to the extreme event of 2003. The site effect was only significant for recovery whereas resilience was not significantly affected by either of them (Table S3). There was no significant pairwise difference among provenances when assessing resilience and recovery of tree growth in response to the extreme event of 2003 ( Fig. 4D and E, left), but Ve (from the central cluster) and La (from the eastern cluster) showed significantly lower and higher levels of resistance in comparison to the other provenances, respectively (Fig. 4F left). Cluster had a significant effect on resilience and recovery but not on resistance (Table S4). Trees from Central Europe (central cluster) showed significantly higher levels of resilience than trees from the Pyrenees and higher levels of recovery in comparison to trees from eastern cluster. Trees from eastern Europe displayed slightly, but not significantly, higher levels of resistance than trees from the other two clusters, mirroring, nevertheless, a longitudinal pattern.
Provenance and site as well as cluster had a significant effect on correlations between tree growth and temperature of May (Tables S3 and  S4). Correlations between tree growth and temperature of May at the common gardens showed significant pairwise differences among provenances (Fig. 4G left): trees from Le displayed lower correlation compared to trees originating from Py, Ca, and St. At the cluster level, correlations between tree growth and May temperature were found to be higher for trees from the Pyrenees and eastern cluster than trees from central cluster. Within the higher elevation site (Bourrignon), trees from Py showed significantly lower correlations between treering width and July temperature than trees from Ba and Ri (Fig. S9), but no significant pairwise difference in correlations with July precipitation was observed among the provenances. However, at the cluster level, growth of trees from the Pyrenees was significantly more Fig. 3. Climate-growth correlations of the 18 studied provenances. Bootstrapped monthly Pearson correlations between mean temperature and precipitation with i) common garden chronologies (black dots) and ii) mean site provenance chronologies (coloured dots). Red dots overlaying black dots represent significant correlations only for common garden chronologies (p < 0.05). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) correlated to both climatic parameters compared to growth of trees from the other two clusters.

Signs of climate adaptation
To assess potential adaptation of the provenances to the climate of seed origin , Pearson correlations between the selected traits and 19 bioclimatic variables were calculated. Only seven of the explored bioclimatic variables displayed significant correlation with the traits. DBH and tree height showed a negative correlation with mean temperature of the wettest quarter, maximum temperature of the warmest month, and isothermality (ratio between mean diurnal range of temperature and temperature annual range). This means that provenances from areas with relatively cool growing seasons showed higher DBH and heights than trees from warm provenances (Fig. 5A). Levels of autocorrelation were unrelated to the investigated bioclimatic variables. Few significant correlations between resilience components and climate of seed origin were found (Fig. 5B). Growth resilience and recovery in response to the extreme event of 2003 were negatively correlated to precipitation seasonality, meaning that provenances from regions with less variable precipitation have higher resilience and recovery values. Additionally, recovery was positively correlated to precipitation of the driest quarter, which indicates that provenances with a rather well distributed precipitation over the year recovered better from extreme events than provenances with less variable precipitation. Resistance was negatively associated with isothermality. Correlations between tree growth and May temperature were positively and negatively correlated with precipitation of the driest quarter and precipitation seasonality, respectively (Fig. 5C). Within the higher-elevation site (Bourrignon), correlations between tree growth and July temperature or precipitation were not significantly associated with any of the bioclimatic parameters (Fig. S10).
Despite the rather low sample sizes (i.e. number of trees per provenance), which resulted in high standard errors and confidence intervals of Q ST (Table 2, Fig. S11), Q ST -F ST comparisons indicated that most trait differentiation among provenances may be a result of directional selection. Tree height had the strongest adaptive signal with the largest Q ST value (0.72) and non-overlapping CIs with F ST . Q ST values of all other traits were higher than overall F ST (0.086) but CIs of Q ST and F ST were overlapping (Fig. S11). In particular, the Q ST -F ST differences for autocorrelation and resistance to the extreme drought event in 2003 were negligible. Heritability (h 2 ) was similar in all traits (range: 2.6-10.5%, Table 2). Interestingly, Q ST values were significantly correlated to the correlation of traits with long-term climate at the seed origin ( Fig. S12) for mean annual temperature and total annual precipitation, and very close to significance in isothermality, maximum temperature of the warmest month, and precipitation to the driest quarter (Fig. S12).

Discussion
Overall, our results reveal clear genetically-based differences in growth traits among provenances and particularly among genetic clusters from the entire natural distribution range of silver fir. Patterns of growth performance followed a West-East pattern similar to the neutral population genetic structure. Provenances from eastern Europe showed higher long-term growth performance but similar levels of resistance and resilience to the extreme drought of summer 2003 compared to the provenances from Pyrenees and Central Europe. However, the provenances from Central Europe displayed the highest recovery rate after the 2003 extreme drought. The provenance from Pyrenees showed lower long-term growth performance as well as a reduced resilience to cope with extreme droughts. Altogether, our results suggest that the differences found among provenances in growth-related traits likely result from genetic differentiation processes that took place several millennia ago within the European glacial refugia and during the post- Provenances and clusters sharing the same lowercase letter are not different at p < 0.05. DBH, diameter at breast height; Correlation, Pearson correlation between tree-ring detrended series and May temperature. Besides DBH and tree height, the remaining traits are unitless. Provenances and clusters are ordered longitudinally. P, C and E name Pyrenees, central and eastern clusters, respectively. Due to admixture (Figs. 2, S7), provenance Si was not assigned to a cluster and is coloured in black. glacial recolonization, as well as from genetic adaptation to the contemporary local climate.

Neutral population genetic structure and diversity
Our results clearly identified genetically differentiated clusters (as identified by Structure) following a longitudinal West-East pattern that can be associated with the different genetic refugia reported in the literature (Tinner and Lotter, 2006;Liepelt et al., 2009;Cheddadi et al., 2014). The provenance Si represents an exception, since it could not be clearly assigned to one of the clusters, also at higher K-values (Fig. S7). This provenance, situated between the provenances of the central and eastern clusters, seems to include admixed individuals from these two clusters, which is in contrast to the literature mentioned above. Our results are also in agreement with other studies showing that silver fir populations from the Pyrenees are genetically distinct and isolated from other populations in Europe (Konnert and Bergmann, 1995;Terhürne-Berson et al., 2004); and with those studies highlighting genetic similarities between populations from southern Italy and the Balkans possibly caused by trans-Adriatic gene fluxes (Piotti et al., 2017). Similarly, our results show that provenances from Central Europe are genetically different from the eastern lineage as they may result from colonization from the Apennine refugia (Muller et al., 2007;Cheddadi et al., 2014). We did not find differences in the genetic structure of provenances from the Carpathians (To, St, Ri and La). These results contrast with previous studies indicating that this region is a suture zone for Balkan and Apennine lineages (Gömöry et al., 2012;Bosela et al., 2016). A better representation of provenances from those regions (Central Italy and Balkans) would be useful to better understand the details of the admixture zones. However, the goal of our study is to highlight differences in growth responses of clusters and provenances across the study range, rather than to resolve detailed colonization patterns since the last glacial maximum. Our neutral genetic analysis also suggests high genetic diversity and low inbreeding within provenances (Table S1). This is, together with low genetic differentiation (Table S2), a typical pattern of wind-pollinated trees with large census population sizes (Petit and Hampe, 2006). Indeed, it has been demonstrated that silver fir is able to establish a highly efficient and far-reaching pollen-mediated gene flow between populations and even refugia (e.g., Liepelt et al., 2002). Slightly lower genetic diversity was evident in five rather peripheral provenances (Ca, La, Py, Ri, and St, Table S1); in three cases (La, Py, Ri) this did not coincide with an elevated and significant inbreeding coefficient, pointing towards historic bottlenecks and genetic drift in these provenances.

Traits under natural selection
The observed differences among provenances suggest a genetic basis of the differentiation in growth performance (tree height and DBH). In contrast to what is expected for locally adapted provenances, we observed that trees from provenances located geographically and climatically near the provenance trials (i.e. Central European seed sources) did not display the highest long-term growth performance. Instead, the provenances from eastern Europe showed the highest values of tree height and DBH. These results suggest that the membership to a specific genetic lineage (as a result of glacial refugia and post-glacial migration) should be accounted for when assessing growth performance of different provenances in addition to accounting for adaptation to contemporary climate conditions . However, our results also suggest that such a pattern of variation is partly governed by temperature-related drivers, since both traits (height and DBH) were negatively correlated to temperature conditions from the region of origin. These results were confirmed by our Q ST -F ST analyses indicating that recent natural selection has also contributed to the differences found among provenances, at least for height. These findings disagree with other studies where silver fir provenances from cold regions had lower total height (Csilléry et al., 2020) or, in case of other conifer species, northern provenances were growing at slower rates than provenances originating from warmer climates (Montwé et al., 2016;Sang et al., 2019)    . Pearson correlations between climate of the seed origin and (A) growth performance including diameter at breast height (DBH), tree height and levels of autocorrelation (Ac); (B) responses to the 2003 extreme event including resilience (Rs), recovery (Rc) and resistance (Rt) of the extreme drought event in 2003, and (C) long-term climate response namely correlation between tree-ring indices and temperature of May. Asterisks indicate significant correlations at p < 0.05 level.

Table 2
Results of Q ST analyses. Shown are estimates and standard errors (SE) for Q ST , heritability (h 2 ), and additive genetic variance (V A ) for the seven tested traits diameter at breast height (DBH), autocorrelation (Ac), resilience (Rs), recovery (Rc), resistance (Rt), correlation between tree growth and May temperature (TMAY). Trees from the Pyrenees displayed lower levels of autocorrelation in comparison to provenances from the central cluster, meaning that previous year growth has a low influence on current year growth. High levels of autocorrelation are usually linked to complacent growth, meaning that trees are less sensitive to climate fluctuations (Fritts, 1976) as found here for the Central and eastern European provenances. The fact that there were no signs of natural selection when assessing correlations of autocorrelation levels with climate from seed origins and performing Q ST -F ST comparisons suggests that year-to-year growth memory is not a trait under strong natural selection. These results indicate that the lower autocorrelation levels displayed by the trees from the Pyrenees are likely due to the poor genetic pool due to isolation processes rather than natural selection processes. To our knowledge, levels of autocorrelation have never been tested for signs of natural selection, but our findings suggest that, in the case of silver fir, autocorrelation might not be a trait under selection itself but rather a consequence of other physiological traits that could be under selection (e.g., earlywarning signal to detect tree mortality, McDowell et al., 2010;Camarero et al., 2015;Cailleret et al., 2019).
Signs of natural selection of the resilience components to extreme climatic events have been studied for different tree species (Taeger et al., 2013;Montwé et al., 2016;Isaac-Renton et al., 2018;Sang et al., 2019). However, most of these studies did not consider genetic data to confirm their findings (but see Heer et al., 2018;Trujillo-Moya et al., 2018;Depardieu et al., 2020). Our correlation analysis revealed that resilience might have been under natural selection. However, patterns of divergent selection were not detectable based on Q ST -F ST comparisons, although further analyses with larger sample size are needed to confirm such findings. Indeed, the non-significant comparison of Q ST -F ST of resilience slightly disagrees with Heer et al. (2018) where the authors found significant associations of traits with SNPs in plant defense genes, but note that such genotype-phenotype associations do not directly test for signs of natural selection. Other studies assessing adaptive traits in silver fir also highlighted a reduced genetic differentiation in drought resilience among provenances (George et al., 2015), in line with our results. A distinct pattern was observed regarding the resistance capacity of the provenances to withstand the 2003 extreme event. Such pattern also mirrors a West-East gradient, although the differences among genetic clusters were not significant. The negative relationship of resistance values with isothermality could be a potentially confounding result, since isothermality also displays the same West-East gradient as the distribution of the studied provenances. The fact that the Q ST value for resistance is higher (but not significantly) than F ST indicates that further studies might be required to disentangle the role of selection during the post-glacial migration and in contemporary habitat.
Diverse physiological mechanisms underlie the observed differences among provenances in climate-growth correlations. The differences observed in the correlations between growth and May temperature among the provenances suggest a genetic basis of cambial phenology. Although studies assessing the environmental cues governing cambial resumption are relatively scarce compared to the amount of studies investigating leaf phenology, there is strong evidence that temperature is the main driver of xylogenesis among other factors such as photoperiod or water availability (Moser et al., 2009;Rathgeber et al., 2016;Rossi et al., 2016). Indeed, the timing of phenological events is often associated with the delicate trade-off between growth-wise benefits and risk of frost events (Vitasse et al., 2014). Although studies have demonstrated that trees have evolved to leaf out at the time when the statistical recurrence of frost is almost nil (Lenz et al., 2016), this risk is currently increasing in Europe with climate warming (Zohner et al., 2020), and spring frost is assumed to cause massive growth reduction due to the shortening of the growing season and the cost of additional carbon resources to produce new leaves (Augspurger, 2009;Vitasse et al., 2019b). Our results suggest that the timing of growth resumption might have a genetic basis. This is based on the fact that silver firs from the Central European cluster displayed no significant correlation with May temperature suggesting a perfect tuning of the cambial activity with current climate. In contrast, trees from warmer provenances may start their growing season earlier due to a lower thermal requirement as suggested by the significant correlations between growth and May temperature.
Although our Q ST -F ST comparisons were limited by the low sample size, they suggest that some of our measured traits may be under recent natural selection with overall Q ST mostly being higher than overall F ST . The fact that Q ST values and correlation coefficients of traits with longterm climate were often correlated strongly supports this finding. However, other evolutionary processes such as the degree of plasticity in the responses of the provenances also play a role in the capacity of the species to withstand future projected climatic extremes (Valladares et al., 2014). In order to infer such effects, a better representation of the same provenances in both common gardens as well as in additional common gardens located in more extreme and diverse climatic conditions will be required for further investigations. Further studies including reciprocal transplants in a wide geographic scale will also greatly help to strengthen our findings.

Implications for assisted gene flow
Because more frequent, longer lasting, and more severe droughts are forecasted in Europe (IPCC, 2013), long-term growth of several tree species is threatened for the coming decades. Silver fir has been suggested as one of the preferred conifer species in Central Europe for future management strategies partly replacing Norway spruce (Vitali et al., 2017), since genetic maladaptation is unlikely to be a problem for silver fir (Frank et al., 2017) and its deep root system allows to cope with severe droughts. Except for Py and Ve that displayed resilience values <1 after the 2003 summer drought, all provenances were able to regain and maintain their previous growth after the extreme summer drought 2003, suggesting that the species as a whole has the ability to withstand dry spells similar to the one from 2003, even provenances from more humid and cooler climatic conditions. However, this relatively high resistance and resilience to extreme droughts might not be enough to cope with higher frequency of such events as expected in the future. As a matter of fact, recent dieback of silver fir and also of other species such as European beech or Norway spruce was observed in recent years at lowlands in Switzerland, Germany and Austria as a consequence of the extreme summer droughts that occurred in 2015 and 2018 in Central Europe (Schuldt et al., 2020).
It is still unclear which provenances will be able to best cope with projected climate change. Our study provides some insights to support assisted gene flow as an adaptation strategy for climate change. Against the popular conviction that provenances from the southern part of species distribution are potentially more able to withstand forthcoming climate conditions, we observed a reduced growth performance of the individuals from the Pyrenees, which are likely linked to their poor genetic pool due to isolation processes rather than natural selection processes. As a consequence, numerous episodes of tree mortality and forest decline have been observed and projected in the Pyrenees region Gazol et al., 2015Gazol et al., , 2020Sánchez-Salguero et al., 2017). In contrast, assisted gene flow using provenances from the Balkan lineage might be considered as they showed the best longterm growth performance. In particular, the provenance from southern Italy (Ca) was among the best growth performers. This provenance has been earlier identified as preferred material for assisted gene flow (see Larsen and Mekic, 1991;Hansen and Larsen, 2004;Kerr et al., 2015) due to its high growth vigor. However, the higher sensitivity of these provenances to spring temperature suggests a relatively early start of the growing season, which might be detrimental in case of late damaging spring frosts, whose risk remains high in Central Europe (Zohner et al., 2020). Moreover, the eastern provenances showed similar levels of resistance and resilience but reduced recovery in response to the investigated extreme drought year 2003 compared to the local Central European provenances. This indicates that despite higher long-term growth, this genetic cluster might not be an optimal seed source in the long-term assuming that extreme drought events become more severe and more frequent in the coming decades. Further investigations would be necessary to test the ability of these provenances to cope with multiple severe droughts. Overall, our results stress the importance of considering not only the growth performance of the provenances but also the intra-population genetic diversity and the demographic history which may have strongly influenced their ability to cope with climate change.
CRediT authorship contribution statement YV, CR and TW designed the experiment. EM-S, CR, YV and FG carried out the statistical analyses. EM-S, YV, CR, CB and PF helped in the interpretation of the results. EM-S, YV and CR led the writing of the manuscript and all co-authors contributed to the text with substantial inputs.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.