Tree rings provide a new class of phenotypes for genetic associations that foster insights into adaptation of conifers to climate change

Summary Local adaptation in tree species has been documented through a long history of common garden experiments where functional traits (height, bud phenology) are used as proxies for fitness. However, the ability to identify genes or genomic regions related to adaptation to climate requires the evaluation of traits that precisely reflect how and when climate exerts selective constraints. We combine dendroecology with association genetics to establish a link between genotypes, phenotypes and interannual climatic fluctuations. We illustrate this approach by examining individual tree responses embedded in the annual rings of 233 Pinus strobus trees growing in a common garden experiment representing 38 populations from the majority of its range. We found that interannual variability in growth was affected by low temperatures during spring and autumn, and by summer heat and drought. Among‐population variation in climatic sensitivity was significantly correlated with the mean annual temperature of the provenance, suggesting local adaptation. Genotype–phenotype associations using these new tree‐ring phenotypes validated nine candidate genes identified in a previous genetic–environment association study. Combining dendroecology with association genetics allowed us to assess tree vulnerability to past climate at fine temporal scales and provides avenues for future genomic studies on functional adaptation in forest trees.


Introduction
In boreal and temperate ecosystems, a substantial increase in temperature and in the frequency of extreme events (e.g. severe droughts or heat waves) is projected for the upcoming century (IPCC, 2013). Poleward climatic niche shifts are expected (e.g. McKenney et al., 2014), therefore raising questions about the ability of trees to cope with those rapid changes (Price et al., 2013). Migration and long-term genetic adaptation are expected to be too slow to keep pace with climatic niche shifts (Jump & Peñuelas, 2005;Aitken et al., 2008;Jump et al., 2009;Jezkova & Wiens, 2016). Nevertheless, populations may persist in their current location depending on species sensitivity and adaptive capacity to future environmental changes (Aubin et al., 2016). To better predict the fate of tree species, knowledge is needed about the climatic constraints affecting their growth/development and the key genetic and physiological mechanisms involved in the response to those constraints (Aubin et al., 2016;Urban et al., 2016). These knowledge gaps challenge the forecast of species' responses to shifting climatic conditions and the implementation of forest management strategies such as assisted gene flow from pre-adapted populations (Aitken & Whitlock, 2013). Without a genuine tree-centered approach based on time-series data to evaluate sensitivity to fluctuating environmental conditions, elucidating the genetic architecture of adaptation to climate throughout a tree's lifespan will remain difficult (Alberto et al., 2013).
Local adaptation in tree species has been documented through a long history of common garden experiments composed of georeferenced genetic material (provenances), where dendrometric traits such as diameter or height growth were used as proxies for fitness (Langlet, 1971;Alberto et al., 2013). The rapidity of climate change advocates the use of those existing common garden experiments (which have already been exposed to warming during the last few decades) to investigate responses to climate change (e.g. Rehfeldt et al., 1999;Sork et al., 2013). In combination with genomic approaches, common gardens represent a unique opportunity to decipher the genomic architecture of local adaptation, and to identify putative genes or genomic regions involved in tree species' responses to climate (Lepais & Bacles, 2014;de Villemereuil et al., 2016).
Two genomic strategies are commonly used to search for genes involved in local adaptation to climate (Sork et al., 2013). The first is the genotype-environment association (GEA) method that looks for correlations between genetic markers and environmental variables from where populations/individuals originate (Joost et al., 2007;Coop et al., 2010). However, GEA does not provide information about the phenotypic traits controlled by the candidate genes. Another strategy is to test for associations between genotypes growing in a common environment and traits of interest (genotype-phenotype association; GPA). GPA has the advantage of establishing more direct links between genes and traits potentially under selection (e.g. Eckert et al., 2015). However, in tree species, the list of functional traits that closely reflect adaptation to climate is still rather limited and labor-intensive to evaluate (for a review, see Aitken & Bemmels, 2016). In particular, assessment of dendrometric traits usually consists of punctual measurements (e.g. height or diameter at a certain age) that are probably the cumulative result of many different climatic events and constraints that occurred throughout the tree's lifespan. The use of innovative tools and methods to measure traits that reflect how and when climate exerts selective constraints (Rahaman et al., 2015;Urban et al., 2016) is thus needed to improve our understanding of species sensitivity and adaptive capacity to climate.
Through the investigation of interannual tree-ring variations, dendroecology allows the quantification of climatic constraints that are exerted on trees. For decades, dendroecologists have developed methods to link climate with variations in wood anatomy characteristics (Fritts, 1976;Cook & Kairi uk stis, 1990;Schweingruber, 1996;Gessler et al., 2014;Housset et al., 2015;Girardin et al., 2016;Hartmann & Trumbore, 2016). Indeed, cambial activity is an integrative indicator that accounts for physiological mechanisms such as responses to drought stress, phenology of dormancy and resistance to freezing injuries. Those wood characteristics can be measured using nondestructive techniques to build a retrospective time-series of growth-related traits throughout a tree's lifespan. The influence of abiotic stress can be estimated by quantifying annual growth responses to punctual extreme climatic events (e.g. Montw e et al., 2016). Several recent reports have demonstrated the intraspecific genetic variation in growth response to punctual extreme climatic events (Montw e et al., 2016) or in cumulative growth responses to the average climate (over many years) in common garden experiments (McLane et al., 2011;Taeger et al., 2013;Leland et al., 2016). Most importantly, analysis of the relationships between growth time-series and monthly climatic variables makes it possible to assess sensitivity to climatic constraints on a year-to-year basis. Therefore, it enables the identification of growth-limiting factors and the time of the year (month or season) when trees are the most affected (Fritts, 1976;Girardin et al., 2012). Such tree-ring traits reflecting sensitivity to climatic constraints have never been used as phenotypic traits in the context of quantitative genetics.
Here, we develop a new class of phenotypes that reflect responses to interannual climatic variability across trees' lifespan. We then combine these new climate-related phenotypes with association genetics to establish a link between the genotypes, phenotypes and climatic constraints affecting trees (Fig. 1). The present study demonstrates how tree-ring traits derived from dendroecology can be used to identify and quantify climatic constraints influencing tree growth; study local adaptation in growth responses to such constraints; and identify and validate the genes putatively involved in such adaptation processes. We illustrate our approach using a Pinus strobus (eastern white pine, hereafter P. strobus) common garden experiment located in Quebec that includes provenances from a wide range of environments across the species' distribution area in eastern North America. This species is representative of many temperate and boreal tree species for which the magnitude and speed of climate change are expected to induce widespread range shifts by 2060 (Joyce & Rehfeldt, 2013;Pedlar & McKenney, 2017). Local adaptation in P. strobus was made evident in provenance trial studies, which detected significant among-population variation in classic quantitative traits (e.g. height, bud phenology: see Li et al., 1997;Joyce & Rehfeldt, 2013). More recently, GEA studies found several single nucleotide polymorphisms (SNPs) associated with climate (Nadeau et al., 2016;Rajora et al., 2016).
As a first step, a tree-ring analysis was conducted to provide three different classes of phenotypic traits (Fig. 1). Traits included average annual ring width and density (hereafter type A traits), relative growth variation for a given year following an extreme climatic event (type B traits) and climatic sensitivity (CS) traits, that is, the strength of relationships between annual growth and climatic variation such as temperature and drought (type C traits). In a second step, those traits were correlated with mean annual temperature of the provenance (MAT P ) to search for local adaptation in growth responses to climate (e.g. genetic-environment clines). Third, GPA analyses were performed on those treering traits to identify putative loci involved in growth responses to climatic constraints and to further cross-validate candidate genes that had been identified in a previous GEA study (Nadeau et al., 2016).

Biological material and experimental design
Throughout this paper, we will use the term 'provenance' when referring to the geographic and climatic origin of a population, and the term 'population' when referring to the trees grown from seed sampled at the provenance (one population per provenance; McLane et al., 2011). A list of definitions for abbreviations used throughout this paper is given in Table 1. We studied P. strobus L. trees from 38 provenances planted in a common garden experiment. Provenances were selected to cover the majority of the species' range from Canada and the US ( Fig. 2; Supporting Information Table S1) to maximize the representativity of the intraspecific genetic variability. Nadeau et al. (2015) identified two main genetic groups (southern and northern) probably resulting from separate glacial refugia. Thirty-six provenances were part of the northern genetic group and two provenances were from the southern group (provenances 12 and 13). Seeds were collected between 1976 and 1986 in natural stands, stored at À18°C, germinated in a glasshouse and cultivated for 3 yr in a tree nursery before their establishment in the common garden experiment at the Valcartier arboretum (Quebec, Canada; 46.94886°N, 71.4962°W) in 1989. The experiment was designed as two randomized complete blocks, with up to n = 16 trees per seedlot and blocks arranged in four-tree row plots. Spacing between trees was 2.0 9 1.5 m (thinned to 2.0 9 3.0 m after the 2002 growing season) in block 1, and 2.0 9 2.0 m in block 2. On average, three trees (SD = 0.96) per population were sampled in each block. A total of 233 trees were sampled, representing two to ten trees per population (mean = 6.13, SD = 1.51), depending on the population survival rate (Table S1).

Meteorological data
Monthly means of daily temperature (T ), mean monthly precipitation (P ) and monthly means of daily drought code (DC) timeseries at the Valcartier experimental site were obtained for the 1993-2014 period using the BIOSIM v.10 software (R egni ere & Bolstad, 1994). Daily temperature and precipitation data were interpolated from the eight closest weather stations, adjusted for elevation and location differentials with regional gradients, and averaged using inverse distance weighting (1/d 2 , where d is distance). Monthly values were computed from the average of the daily quantities. DC is an index for the amount of moisture stored in deep organic layers and soils and takes into account snow arrival in autumn and thawing in spring based on modeled snow accumulation (Terrier et al., 2013). The drought season is defined as the snow-free period, which ranges from c. May to October. DC is unitless and ranges from 0 to c. 800, with low values indicating low drought risk and high values indicating high drought risk. Collinearity amongst these different monthly climate variables (temperature, precipitation and drought) is   (Table S2). Mean annual temperature and precipitation of each provenance (MAT p and MAP p , respectively) were also computed and averaged over the 1961-1990 period (i.e. the period corresponding to seed collection) using BIOSIM v.10.

Tree-ring data
At the end of October 2014, a 5-mm-diameter core was taken (from bark to pith) from each tree at 1.3 m above ground using a Pressler increment borer. Cores were kept frozen until being airdried, with no risk of developing mold. They were later sawn to a 1.7-mm thickness longitudinally using a twin blade sleeve to obtain a smooth, thin surface. The same sample was used for dendrochronological and densitometric analyses. Annual rings of each core were visually cross-dated with skeleton-plots and pointer-year identification under 4.59 magnification. Samples were then conditioned at a temperature of 20°C and relative humidity of 40% to maintain an equilibrium moisture content of 8%. For dendrochronological analyses, ring-widths were measured from high-resolution scanned images of the cores using the software Coo-Recorder (Cybis, Saltsj€ obaden, Sweden) at a precision of 10 À2 mm. Cross-dating was validated using the software CDENDRO (Cybis) and COFECHA (Holmes, 1983). The ringwidth measurement series were converted into basal area increment (BAI) chronologies using the R library 'DPLR' 1.6.3 (Bunn, 2008). Average annual BAI for the 1993-2014 period was calculated for each population from the raw data to obtain absolute levels of growth.
Two BAI metrics were used in downstream analyses: the raw absolute BAI and standardized residual BAI scores. Eliminating the intrinsic age-and size-related growth trends in BAI makes it possible to address annual growth variability independently of biological age or tree size. Here, raw BAI chronologies, which were not normally distributed, were transformed using a logarithmic scaling and detrended using modified generalized www.newphytologist.com exponential functions ('Hugershoff' type; Warren, 1980) with the software ARSTAN v.41d XP (Cook & Holmes, 1986). Additionally, an autoregressive model was applied to remove temporal autocorrelation in the time-series, and residuals were extracted. Annual averages of these standardized residual BAI scores were computed for each population using a robust biweight mean, which removes outliers attributable to disturbance or endogenous site factors (R library 'DPLR' 1.6.3; Bunn, 2008). The mean interseries correlation (rbar) computed from all of the standardized residual tree-level BAI time-series (Wigley et al., 1984) was rbar = 0.80 (SD AE 0.17). The expressed population signal obtained for the residual tree-level series was 0.99, which was higher than the commonly used threshold value (0.85) for growth-climate correlation analysis (Wigley et al., 1984). The density of each thin core slice was measured using a Quintek X-ray measuring system at Universit e Laval (Quebec, Canada). Cores were scanned in 0.04-mm steps, producing highresolution density profiles. Average annual density was determined for each tree ring and the correct identification of annual ring borders was verified using the previously cross-dated ringwidth measurements. Mean ring density was calculated for each population by averaging the annual ring density of all trees in a population for the 1993-2014 period. Then, for each core, average ring-density time-series were detrended using modified generalized exponential functions and temporal autocorrelation was removed as described above. Finally, annual averages of these standardized residual density scores were computed for each population as described above.
Three classes of tree-ring traits derived using a dendroecological approach Traits were classified according to the following nomenclature (A, B and C; see description in Fig. 1c): (A) mean xylem BAI (meanBAI 1993(meanBAI -2014 ) and average density (avDens 1993-2014 ) over the 1993-2014 period, corresponding to classic measures used in forestry; (B) growth variation for a specific year following an extreme climatic event; and (C) climatic sensitivity (CS) traits. Type A traits are long-term averages that probably result from many climatic events or constraints on growth throughout a tree' lifespan, whereas type B and type C traits reflect growth responses to climate on a year-to-year basis. To illustrate type B traits, we analyzed the standardized residual BAI scores for the year 2003 (BAI 2003 ) because we observed large among-population variation during that year (Fig. 3), which could be related to lower temperatures observed in autumn 2002 and spring 2003 (Figs S1, S2). The standardized residual BAI 2003 scores indicate the relative growth variation of a population in comparison to its average growth level. Therein, a value of BAI 2003 < 1 indicates a lower growth in 2003 as compared with the population's average growth level for the 1993-2014 period or, conversely, a BAI 2003 > 1 indicates a higher growth than average.
CS (type C) traits were determined as follows. For each population, correlation analyses were computed between the annual standardized residual BAI scores, averaged at the population level, and the matrix of monthly climatic data (from 1993 to 2014) interpolated at the location of the common garden (Biondi & Waikul, 2004). Long-term trends in climatic series were linearly detrended to obtain unbiased data of interannual climatic variations. Correlation analyses were performed for the period spanning from June of the year preceding ring formation to September of the year concurrent with ring formation to account for lagged effects from year to year (e.g. previous autumn freezing injury). The coefficients of correlation between growth metrics and climatic variables hereafter constitute the CS traits. Significance in correlation analyses was assessed using bootstrapped randomization (n = 1000 simulations): when the 95% confidence interval included zero, the null hypothesis of 'no relationship' could not be rejected. CS traits were also calculated at the individual-tree level to be used as phenotypes in GPAs (see below). As is routinely done in dendroecology, a response function analysis was also used as a complement to quantify the influence of climate on tree growth (Cook & Kairi uk stis, 1990). This analysis is a multiple regression technique that uses the principal components of monthly climatic data to estimate growth-climate associations. This complementary analysis, made in support of the correlation analysis, controls for collinearity of monthly climate predictors; results are presented in the Supporting Information section and are not discussed herein (Figs S3, S4). All analyses were conducted using the R library 'TREECLIM' 1.0.11 (Zang & Biondi, 2015). CS traits of the average density metric were also determined by repeating all of the above analyses on the annual standardized residual density scores, averaged at the population level.

Local adaptation in growth responses to climate
We screened tree-ring traits (type A, B and C) exhibiting amongpopulation variation that could be correlated with the climate of the provenance, an indication of local adaptation in growth responses to climate (Sork et al., 2013). We used MAT p as an indicator because it was highly correlated with the number of growing degree-days above 5°C (r = 0.97), which was found to be the best predictor of genetic variation in growth potential among P. strobus populations (Joyce & Rehfeldt, 2013). MAT p is correlated with MAP p in our study system (Spearman r = 0.42, with P = 0.030 and df = 24, taking into account spatial autocorrelation; Dutilleul, 1993), and hence it also partially reflects a moisture availability gradient. First, Mantel tests were conducted to determine if similarities between population-averaged timeseries of standardized residual BAI scores, as measured using pairwise correlations between time-series, and geophysical characteristics of the provenances, were related ('mantel.rtest' function in the ADE4 R package; Dray & Dufour, 2007). Next, for type A and B traits, the relationship with provenance MAT p was tested using both linear and second-degree polynomial regressions. We tested for the effect of block sample size by including a covariate in these models, which was defined as the sample size ratio of block 1 vs block 2; the effect of the covariate was not significant. Model selection was based on the Akaike information criterion (Burnham & Anderson, 2010

Research
New Phytologist traits (type C), a linear regression model provided the best fit. To reduce the risk of generating spurious results, we only tested the CS traits for which growth-climate correlations were significant in at least five populations. This threshold was arbitrarily set to allow sufficient interpopulation variance in the CS traits for the detection of genetic associations.

Genetic data acquisition
For genotyping, DNA from each of the 233 sampled trees was extracted from foliage using the protocol described by Nadeau et al. (2015). Leaf sampling was conducted independently from the sampling of Nadeau et al. (2015Nadeau et al. ( , 2016: only 21 out of 133 populations sampled in Nadeau et al. (2015Nadeau et al. ( , 2016 are common to both datasets (16%), and samples within populations are largely different (c. 4% in common). The 233 trees were genotyped on the same set of 153 SNP markers as in Nadeau et al. (2015) previously developed in a Sequenom iPlex Gold MassARRAY (Agena Bioscience, San Diego, CA, USA). This array included SNPs from 52 candidate genes for growth and phenology previously identified in Picea glauca and 68 noncandidate genes randomly distributed across the genome (for more details see Nadeau et al., 2015Nadeau et al., , 2016. Using this array and an independent set of populations distributed across the range of P. strobus, Nadeau et al. (2016) identified 26 candidate SNPs for local adaptation using GEAs and F ST outlier tests. In this study, eight trees and six SNPs failed during the genotyping step and were discarded from subsequent analyses. SNPs with a minor allele frequency < 0.03 were also discarded. This resulted in 225 trees and 128 SNPs successfully genotyped (SNP genotyping call rate: min. = 90.7%, mean = 98.9%), including 21 of the candidate SNPs identified by Nadeau et al. (2016). All of the 21 candidate SNPs were detected by at least one GEA method, and two SNPs were F ST outliers (Nadeau et al., 2016;Table S3). No trees were discarded based on the genotyping call rate (min. = 78.9%; mean = 98.9%). GenBank accession numbers of the genes containing the 128 SNPs are given in Table S3.

Genotype-phenotype associations
We searched for genes involved in growth-sensitivity to climate by performing GPA analyses on tree-ring traits. GPAs were tested at the individual-tree level for the 128 SNPs using TASSEL v.5.0 software (Yu et al., 2006;Zhang et al., 2010). We only tested the tree-ring traits that showed a highly significant correlation with MAT p at the population level (P < 0.01, see Results). TASSEL tests for associations between genotypes and phenotypes for each SNP-trait combination. The full linear mixed model used can be expressed as follows: where y is the vector of observations for a tree-ring trait (i.e. the phenotypes), b is the vector of fixed block effects, a is the vector of fixed SNP (genotypes) effects, v is the vector of fixed population structure effects (described below), l is a vector of additive or 'background' polygenic effects, which takes into account relatedness among individuals, and is the vector of random residuals. X, S, Q and Z are incidence matrices relating y to b, a, v and l, respectively. Random additive effects and errors are assumed to be normally distributed with l~N(0, 2Kr 2 a ) and ~N(0, Ir 2 ), where r 2 a is the additive genetic variance, r 2 is the residual variance, K is the kinship matrix and I is the identity matrix. The kinship matrix K was computed from the 128 SNPs using the method of Ritland (1996) implemented in the program SPAGeDi v1.5a (Hardy & Vekemans, 2002).
Population structure can create spurious associations between genotypes and phenotypes (false positives). Nevertheless, when population structure covaries with adaptive genetic variation as in P. strobus (Nadeau et al., 2016), over-adjustment for population structure may decrease power and result in false negatives (Sork et al., 2013). Therefore, three models were tested with TASSEL: no correction (i.e. omitting the population structure and kinship effects); correction for population structure (Q model, i.e. omitting the kinship effects); and correction for population structure and for kinship (Q + K model, i.e. full model). To obtain the population membership matrix (Q), we performed a discriminant analysis of principal components (DAPC; Jombart et al., 2010) using all 225 genotyped trees from this study and 821 additional trees (133 provenances) sampled across the entire range and genotyped on the same set of 128 SNPs taken from Nadeau et al. (2015). DAPC maximizes among-population genetic variation while minimizing within-group genetic variation. Therefore, this method is well suited for GPA analyses because population structure is removed by including the group memberships identified by DAPC as covariates in the TASSEL model (Jombart et al., 2010). A complete description of the DAPC method can be found in Fig. S5. We used the membership probabilities for k = 3 clusters (k À 1 = 2 covariates) as the Q matrix to correct for population structure in the TASSEL analyses (Fig. S5b). The three clusters detected were similar to those in Nadeau et al. (2015): a southern and a northern group with weak west to east differentiation (Fig. S5c). To summarize results from the three models (no correction, Q and Q + K), each SNP-phenotype association was classified as follows: 'very likely' when significant (P < 0.05) for all three models; 'likely' when significant for two out of the three models; 'uncertain' when significant for only one model; and 'no evidence' when not significant in the three models (Girardin et al., 2016). While summarizing results across models may help controlling for errors, the likelihood of false discoveries is significantly raised under multiple testing. Because we tested a small number of SNPs (128), we preferred to reduce the risk of false negatives (i.e. not detecting true associations) by using P-values that were not corrected for multiple testing. Instead, the cross-validation of the SNPs detected by TASSEL with the 21 candidate SNPs found in Nadeau et al. (2016), which used independent samples and different methods, should reduce false positives and make a stronger case for those candidate genes. GPAs excluding the two southern group populations (n = 221) gave similar results (Table S4).
To look for functional categories, genes containing SNPs involved in significant GPAs were annotated using BLASTX against the RefSeq database (https://www.ncbi.nlm.nih.gov/ refseq) using a threshold E-value of 1 e À10 . Given that two P. strobus sequences had no significant hit on the RefSeq database, BLAST analyses were performed against the P. glauca gene catalog (tBLASTX, E-value of 10 À10 ; GCAT 3.3; Rigault et al., 2011) and the P. strobus transcriptome (BLASTN, E-value of 10 À20 ; Ingo Ensminger (University of Toronto), personal communication). Complementary information regarding their putative function was then deduced from P. glauca and P. strobus best-ortholog TAIR annotations (TAIR, Arabidopsis thaliana database, https://www.arabidopsis.org/index.jsp). SNP annotations (synonymous, nonsynonymous) were deduced as previously described by Nadeau et al. (2016).

Results
High variability was observed for annual averages of raw-absolute BAI among P. strobus populations (Fig. 3a). Globally, time-series of population-averaged standardized residual BAI scores were well intercorrelated; the mean of the pairwise correlations was 0.66. Yet, the values of these pairwise correlations tended to be higher when computed between provenances that had smaller geographic distances (Mantel test observed correlation r = À0.30, P = 0.019) and smaller MAT p differences (r = À0.39, P = 0.003).

CS traits
Analysis of interannual variability of population-averaged standardized residual BAI and density scores in relation to climatic fluctuations enables the identification of climatic constraints acting upon ring formation. The coefficients of correlation between growth metrics and the climatic variables hereafter constitute the CS traits (type C; Figs 4, 5; see Figs S3, S4 for response function analyses). The positive influence of a lengthening of the growing season on BAI was made evident by the positive correlations between population averages of standardized residual BAI scores (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) and October temperatures during the autumn of the year preceding ring formation (CS BAI -Temp.Oct (tÀ1) ), and with May temperature of the year concurrent with ring formation (CS BAI -Temp.May (t) ) (Fig. 4a). Furthermore, standardized residual BAI scores were negatively correlated with the numbers of freezing days occurring in October of the year preceding growth (CS BAI -Freeze.Oct (tÀ1) ) and in June of the year of ring formation (CS BAI -Freeze.Jun (t) ) (Fig. S6a). In addition, the effect of heat

Research
New Phytologist stress on BAI was revealed by the negative correlations between standardized residual BAI scores and July temperature of the current growth year (CS BAI -Temp.Jul (t) ; Fig. 4a). Negative correlations with July and August DCs (CS BAI -Drought.Jul (t) and CS BAI -Drought.Aug (t) ) were also observed, highlighting the influence of summer drought stress on BAI (Fig. 5a). In parallel, standardized residual BAI scores were positively correlated with current-year July precipitation (CS BAI -Prec.Jul (t) ; Fig. S7a). Individual-tree level CS traits were consistent with the CS traits computed at the population level (visual inspection, Fig. S8). As for ring density, a negative effect of summer drought during the year preceding growth was detected (Fig. S9). Tree-ring density was positively influenced by high temperatures during the growing season (May, June and July; Fig. S10).

Local adaptation in tree-ring traits
The raw BAI averaged by population over the entire period (absolute meanBAI 1993-2014 , type A trait), expressed as a function of MAT p , exhibited a second-order curve ( Fig. 3d; adj R 2 = 0.46, P < 0.001). This revealed a high degree of variation among populations and suggested local adaptation to the mean annual temperature of the provenance. According to that model, the bestperforming population was nonlocal. The MAT p of the population that had the highest growth performance was 6.6°C, which is higher than the MAT in the Valcartier common garden experiment (3.95°C). Populations from both warmer and colder provenances exhibited lower growth. No significant relationship was observed between absolute avDens 1993-2014 and MAT p (Fig. 3c). However, note that the two populations from the southern group have a much higher tree-ring density (Fig. 3f). Among-population variation was also observed for climaterelated tree-ring traits (types B and C). The reduction in growth in 2003(BAI 2003 , type B trait) was stronger for the populations from warmer provenances and absent in some cold-climate provenances (Fig. 3e). Growth in 2003 was particularly reduced for the two southernmost provenances, corresponding to the southern genetic group. Among-population genetic variation was found for CS traits using growth-climate correlations (type C traits). Among traits related to heat and drought stresses, CS BAI -Temp.Jul (t) (Fig. 4) and CS BAI -Drought.Aug (t) (Fig. 5) were positively correlated with MAT p , and CS BAI -Prec.Jul (t) (Fig. S7) was negatively correlated with MAT p . Among traits reflecting cold tolerance, CS BAI -Freeze.Oct (tÀ1) and CS BAI -Freeze.Jun (t) (Fig. S6) were negatively correlated with MAT p , whereas CS BAI -Temp.Oct (tÀ1) (Fig. 4) was positively correlated with MAT p .

CS and genotype-phenotype associations
We searched for genes related to local adaptation in growth responses to climate by testing those tree-ring traits that showed a significant (P < 0.01) correlation with MAT p : meanBAI 1993-2014 (type A), BAI 2003 (type B) and five CS traits (type C): CS BAI -Drought.Aug (t) , CS BAI -Prec.Jul (t) , CS BAI -Temp.Jul (t) , CS BAI -Freeze.Oct (tÀ1) and CS BAI -Freeze.Jun (t) . Diameter at breast height (DBH 21 ) and height measured at age 21 (H 21 ), referred to as dendrometric phenotypic traits in the present study, were also tested. GPA analysis revealed significant associations (P < 0.05) with dendrometric traits, as well as with the three types

Research
New Phytologist of tree-ring traits for both candidate and noncandidate SNPs (Table S3). For dendrometric traits (DBH 21 and H 21 ) and type A tree-ring traits (meanBAI 1993(meanBAI -2014 ), the proportion of significant SNPs was similar between candidate and noncandidate SNPs (Table S5). For all the climate-related type B (BAI 2003 ) and type C traits (CS traits), we found a greater proportion of significant candidate SNPs as compared to noncandidate SNPs, although the proportion of SNPs classified as very likely was greater for noncandidate SNPs for CS BAI -Prec.Jul (t) and CS  .
Among the 21 candidate SNPs tested (Nadeau et al., 2016), 12 were confirmed in this study (10 genes; Fig. 6; Table 2): three SNPs (three genes) were significantly associated with dendrometric traits DBH 21 and H 21 , while 11 SNPs (nine genes) were associated with tree-ring traits (types A, B and C). As many as 10 SNPs (eight genes) were associated with the climate-related type B and type C traits (seven SNP-trait associations classified as very likely, five likely and eight uncertain). Traits reflecting responses to specific freezing events (BAI 2003  The annotation of the 12 candidate SNPs confirmed in this study revealed a good match with the known functions of homologous genes in other plant models (Table 2). Eight SNPs were located in coding regions and three were located in a noncoding region (intron). Among those located in coding regions, three SNPs were identified as synonymous and five were identified as nonsynonymous. All the putative gene products had functions ranging from transcription factor activity to regulation of signaling molecules.

Discussion
This study illustrates how the combined analysis of tree-ring traits derived from dendroecology and genotypes from a common garden experiment can contribute to understanding the complex relationship between genotypes, phenotypes (growth responses embedded in the annual rings) and interannual climatic fluctuations. Indeed, the new class of phenotypes obtained provided useful information to identify climatic constraints for a given tree species. We illustrated this approach by an examination of the CS traits shaping the growth of an important conifer in eastern North America's temperate forests. We quantified the sensitivity of populations to climatic-related constraints across their lifespan and identified the critical times of the year (or previous year) when those constraints affected growth. We then highlighted the clinal variation along MAT p for those traits. This approach made it possible to identify candidate genes for local adaptation in growth responses to climatic variability, and to Fig. 6 Tests of genotype-phenotype associations (GPAs) between traits and candidate loci from a previous genotype-environment association (GEA) analysis (Nadeau et al., 2016): (a) dendrometric traits (DBH 21 and H 21 = diameter at breast height and height at 21 yr), and (b) tree-ring traits of types A, B and C. For each SNP-trait pair, the association was tested with different levels of corrections: without correction, with correction for the genetic structure (Q), and with correction for both genetic structure and kinship (Q + K). The significance of the associations was classified as very likely (when significant for three levels of correction), likely (when significant for two levels) or uncertain (when significant for only one level of correction). BAI, basal area increment; CS, climate sensitivity trait; Temp, temperature; Prec, precipitation; Drought, drought code; Freeze, number of freezing days during the month. Signaling and plant defense (Afzal et al., 2008) Diverse abiotic stresses in herbaceous plants (Afzal et al., 2008) Dendrometric, Signaling and plant defense (Afzal et al., 2008) Diverse abiotic stresses in herbaceous plants (Afzal et al., 2008) New Phytologist (2018)

Research
New Phytologist cross-validate candidate loci obtained from an independent GEA study.

Tree sensitivity to climate
The dendroclimatic analysis pointed to traits reflecting CS for each population (type C traits). At our study site, growth was mainly limited by very low/freezing temperature during the dehardening phase in spring and cold acclimation phase in autumn, when trees are probably the most vulnerable (probably causing freeze injuries; Howe et al., 2003); and by summer heat and drought stress.
The positive influence of a longer growing season on growth is shown by the positive effect of previous October and current May temperatures (Fig. 4), as well as the negative effect of previous October and current June number of days below 0°C (Fig. S6). Warm spring conditions result in an earlier cessation of dormancy and earlier cambium reactivation, which leads to a longer growth period (Begum et al., 2013). Similarly, a warm autumn during the year before ring formation translates into more time available for accumulation of nonstructural carbohydrate reserves (Hartmann & Trumbore, 2016) and for the division of leaf primordia meristematic cells in buds (Fritts, 1976;Macey & Arnott, 1986). This influence of temperature in combination with photoperiod on the growth of the tree populations under study could be confirmed by studying bud and cambium phenology (e.g. Perrin et al., 2017).
Growth was limited by summer heat or drought, with a negative correlation with July temperature (Fig. 4) and July and August drought severity (Fig. 5), and a positive correlation with July precipitation (Fig. S7). Warm summers can have a direct negative effect on growth by increasing the ratio of carbohydrates allocated to maintenance respiration, while hydric stress reduces both water intake and photosynthesis through stomata closure (Hartmann & Trumbore, 2016). The correlations with climate that we found for this common garden were consistent with those reported for natural stands in other regions (e.g. Marchand & Filion, 2011;Chhin et al., 2013;Girardin et al., 2016).
The study of tree rings makes it possible to precisely determine the influence of yearly climatic anomalies on growth during a tree's life cycle; this information cannot be obtained from conventional approaches used in tree genetics, such as the study of tree height or diameter at a given age. It also has the advantage of accounting for time-lagged effects of climate, which is particularly relevant for conifer species showing a determined growth pattern. For instance, autumn temperatures may influence the timing of budburst of the subsequent year, which in turn can impact tree fitness (Heide, 2003).

Tree adaptation potential to climate change
Comparing the populations' responses to climate (type B and type C traits) in the common garden revealed existing genetic variation for tolerance to freeze and drought. The positive influence of a warm autumn during the year preceding ring formation (CS BAI -Temp.Oct (tÀ1) ) was stronger for populations from southernmost provenances and significantly increased with MAT p Anthocyanin biosynthesis, cell proliferation and differentiation (Yamada et al., 2011;Xiang et al., 2015) Diverse abiotic stresses (Ji et al., 2016) B, C N-040 2_4107_01 S Photosynthetic electron transfer (Ishihara et al., 2007) Diverse abiotic stresses in herbaceous plants (Gururani et al., 2015) Known functions and responses to abiotic stresses of homologous genes in other plant species are reported in the two right-hand columns of the  (Fig. 4a,b). Furthermore, the negative effect of early autumn freeze (CS BAI -Freeze.Oct (tÀ1) ) and late spring freeze (CS BAI -Freeze.Jun (t) ) was greater for warm provenances (Fig. S6). The timing of bud set and bud flush are two important traits involved in the adaptation to local climate that are known to be under strong genetic control in forest trees (e.g. Billington & Pelham, 1991;Savolainen et al., 2004;Bennie et al., 2010;McKown et al., 2014). Populations from colder provenances have probably evolved to set buds and enter dormancy earlier, thereby minimizing the risk of early autumn freeze injuries. Clinal variations in cold hardiness were also reported for P. strobus (Lu et al., 2003), thus supporting our findings that warmer provenances were more likely to be impacted by freeze injuries and to exhibit reduced growth. Genetic variation related to cold hardiness was further demonstrated by among-provenance variations in the standardized residual growth of 2003 (BAI 2003 ) following lower than average temperatures (e.g. À2°C) during autumn 2002 and spring 2003 (type B traits). The relative growth reduction following this event was significantly more pronounced for populations from warmer provenances (Fig. 3e), with the two southern group provenances exhibiting a very severe growth reduction. Local adaptation in response to drought was suggested by the correlation between MAT p and among-population variation of CS BAI -Temp.Jul (t) (Fig. 4b), CS BAI -Drought.Aug (t) (Fig. 5b) and CS BAI -Prec.Jul (t) (Fig. S7). Provenances from warmer locations were less sensitive to summer drought in the common garden, suggesting that those provenances were more adapted to droughts (in agreement with Savva et al., 2002;McLane et al., 2011;Taeger et al., 2013;Montw e et al., 2016). Our study demonstrates the great potential of using CS derived from tree rings sampled in common gardens to investigate trees' adaptive capacity to climate change, especially when testing provenances for suitability in future warming climates.

Genetic basis of climate sensitivity
GPA analysis conducted on a combination of dendrometric and tree-ring traits allowed the detection of significant associations for both candidate and noncandidate SNPs from Nadeau et al. (2016). While the noncandidate SNPs may also contain relevant genes for responses to climate, our main goal was to see whether we could confirm those candidate SNPs that had previously been detected in genetic-environment associations (Nadeau et al., 2016) and gain further insights about their putative roles in adaptation to climate. For climate-related traits (types B and C), we found a greater proportion of significant SNPs in candidate than in noncandidate SNPs. In addition, the independent validation of candidate SNPs from a previous study using different samples and different methods should reduce false positives and so we further discuss the candidate genes here. GPA analyses confirmed 12 of the 21 candidate SNPs (six very likely, two likely, four uncertain), and 10 of them (eight genes) were associated with climaterelated traits (types B and C). The five strongest candidate SNPs in Nadeau et al. (2016; namely N-029, G-014, M-015, M-016 and M-017), which had previously been detected by at least two GEA or F ST outlier methods, were also significantly associated with climate-related tree-ring traits in this study. Although false positives are possible, many loci were associated with several traits and with different levels of correction for population structure and kinship, thus increasing our confidence in the results.
The associations between candidate SNPs and type B and C tree-ring traits provided information about the climatic constraints affecting P. strobus growth, and thus helped bridge the gap between the genotypes, the phenotypes and the environmental selection pressures. Taken together, our results suggest that the putative functions obtained for the genes detected in GPA associations were plausible, because information found in the literature supports the established associations between the candidate SNPs and the measured traits reflecting either drought tolerance, freeze tolerance or both. Several SNPs were associated with both freeze and drought tolerance traits (three very likely and one likely; Fig. 6). Among those, SNP N-029 was located within a bHlH (basic helix-loop-helix) transcription factor, which is known to confer tolerance to multiple abiotic stresses in higher plants (Table 2). Similarly, the gene containing SNP N-040 codes for a thylakoid luminal protein involved in photosystem II, one of the most susceptible components of the photosynthetic machinery that bears the brunt of abiotic stresses. Other SNPs were specifically associated with either freeze or drought tolerance traits. Interestingly, the GAUT (SNPs M-017, very likely; O-002, uncertain) and the endoglucanase genes (SNP G-004, uncertain), which are involved in carbohydrate metabolism, were specifically associated with growth responses to freezing events (Table 2). By preserving and stabilizing cell structures, carbohydrates associated with the cell wall can enhance freezing tolerance by controlling water movement and ice formation under cold conditions (Tenhaken, 2015). In addition, one SNP (G-014; uncertain), exclusively associated with drought tolerance traits (Table 2), encoded a PLASTID MOVEMENT IMPAIRED1 (PMI1), a protein that mediates the abscisic acid (ABA)-dependent pathway under drought conditions in Arabidopsis thaliana (Rojas-Pierce et al., 2014). Contrary to coldresponsive genes, many water-stress-inducible genes are upregulated by ABA under drought conditions (Yamaguchi- Shinozaki & Shinozaki, 2006). Further gene expression and functional studies aimed at identifying coexpressed gene networks involved in specific biological processes, such as freeze tolerance, would be needed to determine the exact role of those genes in P. strobus and in other conifers (for a review, see Serin et al., 2016). Multi-SNP GPA models (e.g. Segura et al., 2012) should also be tested given that highly metric traits such as growth are probably controlled by a large number of interacting genes, which was beyond the scope of this study given the small number of SNPs tested.

Conclusion
Accurately determining the periods of the year when trees are most vulnerable to abiotic and biotic factors is crucial to understanding the genetic architecture of local adaption. In this paper, we showed that tree-ring traits measured in common garden experiments offer great potential to unravel local adaptation in forest trees through a fine-scale, retrospective and dynamic view

Research
New Phytologist of tree responses to past environmental variations. For the first time, we proposed using a metric (i.e. climate sensitivity traits; CS) that incorporates both the phenotypic trait under study (i.e. growth or wood density) and the response to climate in association genetic studies. We identified important climatic constraints and the critical times of the year when those constraints affect growth, information that cannot be derived from classic measurements such as height and DBH. We believe that this novel approach will provide a valuable contribution to the study of differences in climate sensitivity between provenances in common garden experiments, and more generally to examine forest adaptation to environmental changes.          Table S1 Summary of provenance location, mean annual temperature of provenance, mean annual precipitation and number of sampled trees per provenance