Login using

You can login by using one of your existing accounts.

We will be provided with an authorization token (please note: passwords are not shared with us) and will sync your accounts for you. This means that you will not need to remember your user name and password in the future and you will be able to login with the account you choose to sync, with the click of a button.

Tree species responses to climate change will be greatly influenced by their evolutionary potential and their phenotypic plasticity. Investigating tree-rings responses to climate and population genetics at the regional scale is crucial in assessing the tree behavior to climate change. This study combined in situ dendroclimatology and population genetics over a latitudinal gradient and compared the variations between the two at the intra- and inter-population levels. This approach was applied on the northern marginal populations of Thuja occidentalis (eastern white-cedar) in the Canadian boreal forest. We aimed first to assess the radial growth variability (response functional trait) within populations across the gradient and to compare it with the genetic diversity (microsatellites). Second, we investigated the variability in the growth response to climate at the regional scale through the radial growth-climate relationships, and tested its correlation with environmental variables and population genetic structure. Model selection based on the Akaike Information Criteria revealed that the growth synchronicity between pairs of trees of a population covariates with both the genetic diversity of this population and the amount of precipitation (inverse correlations), although these variables only explained a small fraction of the observed variance. At the regional scale, variance partitioning and partial redundancy analysis indicate that the growth response to climate was greatly modulated by stand environmental variables, suggesting predominant plastic variations in growth-response to climate. Combining in situ dendroclimatology and population genetics is a promising way to investigate species' response capacity to climate change in natural stands. We stress the need to control for local climate and site conditions effects on dendroclimatic response to climate to avoid misleading conclusions regarding the associations with genetic variables.

Introduction

Climate change will alter ecological gradients, notably via pole-ward shifts of isotherms and locally-dependent changes in precipitation patterns (IPCC, 2013). These modifications are likely to change the selection pressure patterns acting on plant distribution, such as cold temperatures and dryness (Gienapp et al., 2008; Chapin et al., 2011). Recent tree growth reductions in the boreal zone put into question species' capacity to cope with the rapid climate change (Beck et al., 2011; Berner et al., 2011). Induced impacts are expected to be most visible at the limits of a species distribution, where populations are further from the species optimum in ecological gradients (Guo et al., 2005; Sexton et al., 2009). It is now admitted that species response to climate change will combine migration with evolutionary and plastic responses (Davis and Shaw, 2001; Savolainen et al., 2007; Alberto et al., 2013). For long-lived sessile species such as trees, the local fate of a population will rely on its ability to either adapt or acclimate (Aitken et al., 2008). The evolutionary response of a population requires the existence of genetic variability for traits involved in responses to climate (Hartl and Clark, 1997; Savolainen et al., 2007; Polechová et al., 2009). Anticipating the impact of climate change on tree species therefore requires a good characterization of the spatial variations of those traits and of genetic features across large biogeographic gradients.

Tree radial growth is a functional trait that reflects trees' responses to climate. Cambial activity is an integrative indicator accounting for physiological mechanisms such as response to drought stress, dormancy phenology, and resistance to frost injuries (Schweingruber, 1996; Fritts, 2001). By investigating the relationships between climatic variability and annual radial growth, dendroclimatology precisely determines the past, and current local climatic constraints on trees (Fritts, 2001; Filion and Payette, 2010). Dendroclimatic studies conducted over large latitudinal transects are therefore a suitable tool to investigate the plastic response of trees to climate (Huang et al., 2010; Lapointe−Garant et al., 2010; Lloyd et al., 2011; Housset et al., 2015). However, extrapolating the current dendroclimatic response to future climate is impeded—among other challenges—by our lack of knowledge on the genetic control of cambial activity (e.g., Huang et al., 2013).

Consequently, there is a growing interest in research bridging tree-ring and genetic (Franks et al., 2014). A few dendroclimatic studies investigating the among-population variations in the growth-climate relationships in common-garden experiments evidenced genetic variation among provenances for the growth-climate relationships and a high plastic component (Savva et al., 2002, 2008; McLane et al., 2011; Taeger et al., 2013; Montwé et al., 2016). While transplant experiments in common gardens are the most direct way of controlling for environmental differences among populations, these experiments are not available for most non-commercial species. In such cases, investigations of tree responses to climate in natural stands are required. Moreover, investigation in natural stands brings complementary information to common-gardens, which are usually limited to even-aged and monospecific stands in a limited set of soil and climatic conditions.

Recent studies combined in situ genotyping and tree-rings data to test for genetic effects on synchronicity (e.g., King et al., 2013; Latutrie et al., 2015) or on growth-climate relationships (e.g., Bosela et al., 2016) across broad spatial transects. For instance, an effect of the postglacial genetic lineages on the growth-climate relationships for European silver fir (Abies alba Mill.) is reported (Bosela et al., 2016). However, this study did not control for the effect of possible confounding environment features (e.g., soil, elevation, etc). Comparing in situ genetic data and phenotypic data is very challenging because phenotypes results from a complex interplay between genetic control and plastic responses to environmental conditions over an ecological gradient (Gienapp et al., 2008; Franks et al., 2014). Numerous dendroclimatic studies report that growth-climate relationships are modulated by environmental variables such as regional precipitation or temperature regimes, soil variables and forest stands characteristics (e.g., Babst et al., 2013; Gewehr et al., 2014; Housset et al., 2015; Brienen et al., 2016). The effect of those ecological gradients could be confounded with a genetic effect on the growth responses to climate and need to be controlled.

In this study, we test for the correlation between population genetic structure and tree-rings responses to climate, while controlling for ecological gradients through an approach combining model selection and variance partitioning. To test this approach, we took advantage of two existing datasets on the boreal tree Thuja occidentalis L. (eastern white-cedar, hereafter “cedar”). The first dataset results from a study of stem annual radial growth carried out at its northern boreal margin evidencing an effect of tree size, soil, regional precipitation, and regional temperature (Housset et al., 2015). The second is a genetic study conducted in the same study area that investigated the cedar neutral genetic structure based on neutral markers (Xu et al., 2012). The objective is to test for a correlation between neutral population genetic structure and tree-rings, while accounting for environmental variables, at both the intra-population and the inter-population levels.

Materials and Methods

Study System

Cedar is an evergreen long-lived and shade-tolerant species native to North America, with a distribution (Figure 1) ranging from James Bay, Canada, in the north to the southern Great Lakes area, USA, in the south (Little, 1971). Our study area is located from the center to the northern leading edge of cedar distribution, in the eastern Canadian boreal zone (Figure 1B). A latitudinal transect was established from 47.3 to 50.0°N and divided into three zones based on cedar abundance: The continuous zone (CZ), where cedars are common, the discontinuous zone (DZ), which marks the northern edge of the continuous distribution where it becomes less common in the forest matrix and, the marginal zone (MZ) where only a few isolated stands were found (Figure 1A). The relative abundance of cedars in each zone was estimated at 55, 9, and 3%, respectively (Figure 1B; Xu et al., 2012; Beaudoin et al., 2014). A total of 24 cedar stands were sampled along this gradient, with 8 sites in the CZ, 7 in the DZ, and 9 in the MZ, with 4 in the east (MEZ) and 5 in the west (MWZ). Those 24 sampling sites will be referred to as populations in the following sections. Cedar species show plastic tolerance to moisture edaphic conditions by occupying poorly drained lowland sites to xeric or rocky sites (Hofmeyer et al., 2009). All cedar stands in the MZ were found in forest swamps or along lake and river shores. Consequently, to avoid sampling bias in the CZ and the DZ, trees were sampled in the same edaphic situations as in the MZ. Over the 1953–2010 period, climate was colder in the MZ, with a mean (±sd) annual temperature (MAT) of −0.03±0.37°C, than in the DZ and CZ, which had a MAT of 1.40 ± 0.25°C and 2.40 ± 0.33°C, respectively (Table 1). The annual sum of precipitation was on average (±sd) 908 ± 53 mm.year−1 in the MZ, 886 ± 27 mm.year−1 in the DZ and 937 ± 56 mm.year−1 in the CZ (Table 1).

Genetic Data

Genetic structure variables were inferred from a recent study (Xu et al., 2012), which analyzed the cedar neutral genetic structure of our 24 populations. For each population, 15–30 trees (Table 1) were randomly selected and used for microsatellite genotyping. The genetic structure was assessed by genotyping at four polymorphic microsatellite loci derived from the phylogenetically closed species Thuja plicata Donn ex D. Don (TP9, TP10, TP11, and TP12). Xu et al. (2012) selected the optimal combination of markers (available at that time) delivering accurate estimates of cedar genetic diversity. The four microsatellite loci showed high discriminating power, with a combined probability of exclusion (PE) greater than 99% and a joint probability of genetic identity (PI) smaller than 10−5 (Peakall and Smouse, 2012).

Several indicators of intrapopulation genetic variability were used in this study. The expected heterozygoty (He) and the average within-population pairwise relatedness between each pair of trees in a population (Genrelatedness) statistics were computed with the software Genalex 6.5 (Peakall and Smouse, 2012). He is a proxy of gene diversity, with values ranging from 0 to 1 for low to high diversity. The allelic richness (AR), representing the alleles number corrected for sampled population size effect, was calculated using the software FSTAT2.9.3. The inbreeding coefficient (Fis) was calculated for each population using the software FSTAT2.9.3.

Genetic variability among populations at the landscape scale was analyzed through a clustering approach using the STRUCTURE v.2.3.2 software (Pritchard et al., 2000). Individuals were assigned to a number of assumptive clusters K, ranging from 1 to 15, with an admixture model and the option of correlated allele frequency (Falush et al., 2003). All parameters were set following the users' manual. To choose an appropriate run length, we performed a pilot run showing that burn-in and MCMC (Markov chain Monte Carlo) lengths of 300,000 each were sufficient to obtain consistent data. Increasing the burn-in or MCMC lengths did not improve the results significantly. Ten replicate runs for each value of K were performed. The most likely value of K was selected by plotting ΔK (see Supplementary Figure 2), following ad-hoc statistics (Evanno et al., 2005). This analysis identified three clusters (K = 3), which are hereafter depicted using three colors: Orange, yellow, and blue (Table 1). The average proportion of assignment of each individual to each cluster in a given population was used as a proxy to describe among population genetic differences at the regional scale.

Tree-Ring Data and Growth-Climate Relationships

The present study focused on the variability in growth sensitivity to climate among populations, but not on the growth itself. As response variables, we used the correlation matrix between climate and growth relationships described in Housset et al. (2015). This matrix was established on the correlations between monthly climatic variability and annual radial growth of cedar from 27 sampling sites, including the aforementioned 24 populations. Annual radial growth increments were obtained and analytical approaches were conducted using classical dendroclimatological procedures. The following paragraph summarizes Housset et al.'s approach used to derive the correlation matrix, and the main variables selected.

During summer 2010, up to 30 cedars were sampled in each site for dendroclimatic analysis. The genetic and dendrochronological samplings were conducted separately but on the exact same locations. Only dominant and codominant trees were sampled to lower the biases due to competition on tree-ring series. Two cores per tree were collected at 1.3 m above ground, sanded, measured, and cross-dated. To maximize the high-frequency growth responses associated with year-to-year climatic variability, each raw ring-width time-series was detrended using 60-year splines with a 50% frequency response (Cook and Kairiūkštis, 1990). After removing the first order autocorrelation, a mean Tree Growth Index (TGI) chronology was computed for each population using the R package “dplR” (Bunn, 2008). These residual chronologies were statistically tested against climatic time-series using bootstrapped correlation and response function analyses over the 1953–2010 period covered by meteorological observations (Biondi and Waikul, 2004). Correlation coefficients are computed between a site residual chronology and the matrix of monthly climatic data in order to represent the average relationship for the studied period. On the other hand, PCR (also known as response function analysis) is a multiple regression technique that uses the principal components of monthly climatic data to assess growth-climate relationships. Climatic time-series were monthly mean temperatures and total monthly precipitation, interpolated from Environment Canada meteorological records at each of the sampling sites using the software BioSIM for 1953–2010 (Régnière and Bolstad, 1994). This analysis revealed that cedar interannual growth variability was positively correlated with May temperature of the current year, but negatively correlated with June temperatures of the current year and with August temperatures of the year preceding ring formation. Climatic conditions of the previous year can affect tree growth, for instance through the process of carbohydrates reserve accumulation (Fritts, 2001). Precipitation also influenced annual radial growth, with a negative correlation with previous October and current-year May precipitation, and a positive correlation with previous June and current-year August precipitation. Figure 2 sums up the correlation scores of cedar radial growth with the seven monthly climatic variables, which will be referred to as the “growth-climate correlation matrix” in the following sections.

FIGURE 2

Figure 2. Growth-climate correlation matrix for the seven main monthly climatic drivers of radial growth interannual variations (columns) and the 24 sites (rows). Months in capital letters code for climate variables of the current year and, small letters for the year preceding ring formation. Colors represent the correlation coefficients scores computed between Thuja occidentalis residual tree-ring index chronologies and monthly mean temperature and monthly precipitation. Significant correlations are represented by dots (•) for the bootstrapped response function and by circles (○) for the correlation function. Sites are presented in order of increasing number of mean annual growing degree-days over the 1953–2010 period, the coldest sites being at the top.

Variability within a Population

Growth synchronicity was evaluated between the trees of each population. A mean chronology was computed for each tree from the spline-detrended chronologies of all its cores. The average pairwise correlation between single-tree chronologies within a population (r¯bt) was calculated (Cook and Kairiūkštis, 1990). This statistic, representing the average correlation between trees, was used as a proxy of growth synchronicity. It is inversely linked to growth variability. To remove biases due to population size differences, we resampled eight trees in each population and estimated this parameter by averaging it over 1000 repetitions. For each resampling, the r¯bt was computed using the R package “dplR” (Bunn, 2008). Growth synchronicity between trees could be influenced by various factors such as differences in DBH or by genetic variability within a population. Variability in DBH was estimated using the standard deviation of tree DBH for each sampled population (sdDBH). Three complementary proxies of intrapopulation genetic diversity were used: Allelic richness (AR, see Supplementary Figure 1), expected heterozygoty (He), and average within-population pairwise relatedness (Genrelatedness). The inbreeding coefficient Fis (see Supplementary Figure 1) was also included in the model selection. The effects of regional climate on growth synchronicity along the latitudinal gradient was tested by adding the annual growing degree-days >5°C (GDD5) and the annual sum of precipitation from April to September (i.e., during the active vegetation stage) averaged over the 1953–2010 period.

We used model selection based on the corrected Akaike Information Criterion (AICc) statistic to determine which variable best explained growth synchronicity. In the model selection approach, each hypothesis is associated with a model that can be compared with other model hypotheses through AICc (Burnham and Anderson, 2010). The model with the lowest AICc is the one with the optimal trade-off between the best fit and the greatest parsimony. A model is considered to be more parsimonious when the number of estimated parameters k is lower. This method selects the hypothesis with the greatest likelihood to fit the observed data while avoiding the risk of overfitting them. A general rule of thumb is to consider only models with a difference in AICc lower than two as compared with the best model (Burnham and Anderson, 2010). Here, each model had growing variability r¯bt as dependent variable and one of the following independent variables: sdDBH, AR, He, Genrelatedness, Fis, annual GDD1953−2010 or Seasonal P1953−2010.

Growth-Climate Relationships Variability among Populations

We performed another model selection procedure to test for covariation between the growth-climate relationships and among population genetic differences (STRUCTURE clustering) at the landscape scale. Sensitivity to climate is defined here as the coefficients of the aforementioned growth-climate correlation matrix. Alternatively, variability in sensitivity to climate may be due to the environmental conditions of each population. Therefore, the analysis included a model with genetic variables and a model with environmental variables having an effect on sensitivity to climate. A full model including both genetic and environmental variables was also included in the model selection. In the genetic model, we used the composition of each STRUCTURE groups (clusters) as explanatory variables to describe among population variations at the regional scale. Four environmental variables having an effect on variability among populations of the aforementioned growth-climate correlation values were identified in the Housset et al.' study (2015) and were included in the environmental model. Those four variables were: Mean DBH of the cored trees, carbon to nitrogen ratio (C:N) in the upper 50 cm, GDD5 and seasonal precipitation (mean sum of precipitation from April to September, i.e., during the current growing season) averaged over the 1953–2010 period. To sum up, the three following models were compared:

where ri is the correlation coefficient of growth with one of the aforementioned monthly climatic variables (columns in Figure 2). Furthermore, for each monthly correlation coefficients, a variance partitioning was computed between the three models to examine the amount of variability (adjusted R2) explained by each group of variables. The overall comparison of environmental vs. genetic variables was analyzed through the calculation of multivariate variance partitioning (Legendre and Legendre, 2012) on the seven monthly climatic variables using the R package “vegan” (Oksanen et al., 2013). Significance of the variance components were assessed through redundancy analysis (RDA) and partial redundancy analysis (pRDA) with 99 999 permutations (Legendre and Legendre, 2012). RDA are multivariate regressions where a set of response variables (here the growth-climate correlation matrix) are regressed against a set of explanatory variables. pRDA are RDA where the effect of the explanatory variables is assessed while controlling for another set of variables. pRDA can be particularly useful to test for the effect of genetic variables while controlling for the effect of local climate or soil conditions for example.

Results

Genetic Variability and Growth Synchronicity

The growth synchronicity between the trees of a given population (r¯bt) was smaller in the MZ than in the DZ (p = 0.025; Figure 3A). No statistical differences in r¯bt were observed between the CZ and the DZ, nor between the CZ and the MZ. Genetic pairwise relatedness was smaller in the MZ than in the CZ (p = 0.049; Figure 3B). Based on model selection (Table 2), the best model to explain growth synchronicity between the trees of a population is the one with Seasonal P1953−2010 as an explanatory variable. Trees from a site receiving more precipitation had a lower synchronicity. Nevertheless, this relation was barely significant (p = 0.053) and poorly explains interindividual growth variability (adjusted R2 = 0.122). Four other models had a difference in AICc (Δ AICc) < 2, with the following explanatory variables: He, annual GDD1953−2010 (Δ AICc = 0.56), AR (1.39) and Genrelatedness (1.68). None of those models had a significant relation with r¯bt. The major part of growth variability between trees was not captured by our variables.

Table 2. Results of the linear model selection explaining the radial growth synchronicity between trees in a given site (r¯bt).

Variability among Populations

Model selection revealed that the model including environmental variables best explained variability among populations in the correlation with temperatures of previous August and current June, as well as with precipitation of previous October, current May and current August (Table 3). According to variance partitioning, environmental variables captured most of the explained variance for those monthly climatic variables (Figure 4). Nevertheless, the genetic model (Equation 2) was the most likely to explain variations among populations in the correlation with current May temperature and previous June precipitation at the landscape scale (Table 3). Although this model was more parsimonious (k = 4) than the environmental one (k = 6), the Δ AICc (5.2 and 5.4, respectively) largely exceeded the difference in k. According to our data, populations with a greater percentage of trees from the blue genetic cluster were significantly (p = 0.019) less positively correlated with May temperature. A greater percentage of “blue” trees in the population was also significantly (p = 0.012) linked to a more positive correlation with previous June precipitation. The genetic structure accounted for the most important part of the explained variance for those two dependent variables (Figure 4), with an adjusted R2 of 0.17 for the correlation with May temperature and of 0.21 with previous June precipitation. The complete model with both the environmental and genetic variables was classified as the best model for none of the monthly climatic correlation coefficients included in the analysis (Table 3), although the adjusted R2 was greater (Figure 4). Overall, multivariate variance partitioning computed on the seven dependent variables (growth correlation with climate variables) revealed that environmental variables only explained 24% of the variance, while genetic structure only explained 5% (Figure 5). The part of variance that was jointly explained by both environmental and genetic data accounted for 43% of the total variance (with 13% overlap), while 57% of it remained unexplained.

TABLE 3

Table 3. Results of the model selection explaining variability in growth-climate correlation coefficients among sites for each of the main climatic variables (columns from Figure 2).

FIGURE 4

Figure 4. Univariate variance partitioning between environmental variables (Equation 1) and genetic variables (Equation 2) to explain variability in growth-climate correlation coefficient among sites for each monthly climatic variable of Figure 2. Black bars represent the percentage of variance explained by genetic variables only, dark gray bars represent the variance explained jointly by environmental and genetic variables, and pale gray bars represent the variance explained by environmental variables only. The sum of the three represents the total variance explained by all the variables (Equation 3).

Discussion

Our combined analysis of the association between population genetics and the dendroclimatic signal of cedar revealed a correlation between radial growth and genetics at both the intrapopulation (growth synchronicity) and the among-population levels (growth-climate relationships). However, our results evidenced a strong confounding effect by environmental variables in both cases, largely because environmental drivers and the neutral genetic structure are inter-correlated.

Intrapopulation Variability

Growth synchronicity between trees of a site (r¯bt) was lower in the MZ, meaning that cedar growth patterns were more variable from one tree to another toward the northern margin of its distribution. Based on model selection, this result was neither due to a greater spatial variability between the trees nor to size heterogeneity. Growth synchronicity was mostly influenced by seasonal precipitation (best model), a drier climate being associated with an increase in synchronicity. This finding is consistent with other studies showing that growth synchronicity was generally stronger under the most stressful climatic conditions (Schweingruber, 1996; King et al., 2013). It may seem surprising that precipitation exerts a common constraint on tree growth in our study system. However, although cedar is found in cold wetlands and shores, Housset et al. (2015) reported a growth limitation due to drought therein owing to the combination of a superficial root system and a lowering of the water table during the summer. Growth synchronicity between trees can be affected by stochastic processes such as abrupt change in competition, injuries or disturbances, especially in uneven aged natural stands, which could not be controlled in our study set and are thus accounted for as “noise.” Moreover, microtopographic variations among trees within a same site may cause differential responses to the same climatic conditions. It is therefore understandable that the explained variance is quite low (R2 = 0.122, p = 0.053) for the relation between r¯bt and Seasonal P1953−2010. The model with the expected heterozygoty (He) had an AICc score very close to that of the model with Seasonal P1953−2010 (Δ AICc = 0.07). Both p (0.055) and R2 (0.119) values were also close to those of the first model. Higher He, meaning higher allelic diversity in the population, is negatively correlated with the intrapopulation growth synchronicity (r¯bt). Growth synchronicity was generally negatively correlated with other indicators of genetic diversity (AR, which is in line with the findings of Latutrie et al., 2015) for Populus tremuloides Michaux. Disentangling the effects of environmental variables and genetics on growth synchronicity would require further investigations. Nevertheless, it remains that the amount of seasonal precipitation appears to primarily drive growth synchronicity, which supports the idea that plastic responses to environmental variables accounted for a great part of cedar intrapopulation growth variations across our transect.

Predominant Effect of the Environment on Growth-Climate Relationships

Model selection (Table 3) and variance partitioning (Figures 4, 5) conducted at the among-population level confirms the predominant effect of stand environmental conditions on growth response to climate. Environmental model had a better fit with the data for the correlation with summer temperatures (August of the previous year), with previous October precipitation and with current May and August precipitation. Unfortunately, there is no common garden for T. occidentalis cedar species to further experiment our empirical results. A provenance trial has shown that the height growth of T. plicata differed significantly among provenances (Jokela and Cyr, 1977; Russell et al., 2003). However, it did not test the effects of provenances on growth response to climate. Indeed, very few dendroclimatic studies have been conducted on trees in common gardens. For instance, a dendroclimatological study of different provenances of Pinus banksiana Lamb. in a common garden experiment in the eastern Canadian boreal forest concluded that there was no difference in the growth-climate relationships between provenances (Savva et al., 2008). In another common garden experiment conducted with Pinus sylvestris L. in Siberia, provenances accounted for only 15% of the explained variation in the dendroclimatic signal, while 85% was due to the environment (Savva et al., 2002). This is consistent with our results: Variance partitioning revealed that the variables of the environmental model accounted for 79% of the explained variance (Figure 5).

Genetic vs. Growth Climate Relationships

The multivariate analysis revealed that genetic variables were significantly associated with among-populations variations in the growth climate relationships (p = 0.004, Figure 5). Yet, the partial RDA revealed that, when controlling for the environmental variables, this relation was no longer significant (p = 0.074). When considering the growth-climate relationships independently, it appeared that the genetic model best fitted the variability in growth correlation with May and with previous June precipitation. According to variance partitioning, the genetic model's unique contribution exceeded that of the environmental model for those two growth-climate relationships (Figure 4). The genetic model was the second best model (Δ AICc < 2) to explain variations in the correlation with June temperature and with May and August precipitation. These findings suggest that environmental variables and the neutral genetic structure may be correlated, as was also suggested by the high amount of shared variance between the two (Figure 5).

Effect of Population Genetic Composition on Growth-Climate Relationships

We observed that the among-population variation in growth-responses to May temperature and to previous June precipitation (Δ AICc = 0) covariated with the genetic composition (K = 3 clusters) of cedar populations. This is convergent with the study of Bosela et al. (2016), who found contrasting growth responses to climate between two fir genetic lineages. Interestingly, those two growth-climate relationships are associated with mechanisms that are known to be under selection such as (1) dormancy phenology and (2) resistance to drought. First, it is generally accepted that dormancy phenology has a stronger genetic control than growth itself in several tree species (e.g., Howe et al., 2003; Alberto et al., 2013). The timing of growth initiation bears an important evolutionary signification, since there is a trade-off between the decreasing risk of frost damage and the optimization of the growing period length (Nienstaedt, 1974; Chuine, 2010). A common garden experiment evidenced a genetic variation in the date of radial growth initiation for Douglas fir Pseudotsuga menziesii (Mirb.) Franco (Li and Adams, 1994; Gould et al., 2012) and confirmed the existence of a genetic control of cambial reactivation to spring temperature. Second, for the response to drought stress, trees that use water more efficiently would be less affected by hydric stress the year before ring formation, because it allows them to store more carbohydrates or to build more cambial cells for the next growing year. Many experiments have already shown a genetic control for traits relative to drought resistance, such as hydraulic conductivity, stomatal number, stomatal density, and belowground/aboveground biomass ratio (Raj et al., 2011; Richter et al., 2012; Alberto et al., 2013; Moreira et al., 2014). Further, Fan et al. (2008) reported local adaptation to precipitation for T. plicata, for both height increments and transpiration efficiency traits. Because our genetic dataset only included neutral markers, (as in Bosela et al., 2016), further investigations including adaptive genetic markers are needed to test for causal relations between genetic and growth responses to climates. However, our findings provide interesting hypotheses to be tested in future research on cedar adaptation to climate.

Implications for Studies Combining Dendrochronology and Population Genetics

Combined in situ investigation of dendroclimatology and population genetics proved to be a promising approach to study the response to climate of tree species such as cedar, for which no common garden trial is available. It made it possible to estimate both the genetic and growth variabilities within natural populations, as well as differences in growth sensitivity to climate among populations at the landscape scale. Our approach based on variance partitioning and model selection evidenced that most of the current variability in growth response to climate was associated with site environmental conditions, suggesting predominant plastic variations across our latitudinal transect. Partial redundancy analysis revealed that the effects of genetic variables are confounded with environmental effects, although stronger correlation with genetic variables were observed for responses to summer drought and spring temperature. This study was an attempt to demonstrate that combining in situ population genetic variable and tree-rings data without controlling for environmental effects may bring to misleading results and ultimately affect studies' conclusions. We argue that combining tree-rings data and population genetics across vast ecological gradients should be developed to assess the influence of demographic and evolutionary histories of non-model tree species, and to anticipate their responses to climate change. Future efforts should include adaptive genetic markers to better capture signals of genetic control and should comprise relevant environmental descriptors that could blur the effects of genetic background on growth-climate relationships.

Author Contributions

JH, CC, MG, FT, and YB conceived the idea of this manuscript. JH and HX collected the data and conducted the lab measurements. JH and MG conducted the statistical analyses. All authors contributed to the interpretation. JH wrote the manuscript with the aid of all authors.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Aurélie Terrier and Nathalie Isabel for their useful comments on preliminary versions of this manuscript, as well as Cécile Leroy, Mélanie Desrochers, Danielle Charron, Véronique Paul, Pamela Cheers, Isabelle Lamarre, and Daniel Lesieur for their valuable help. This study was conducted with the financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through its CREATE program in Forest Modeling Complexity and of the French program Paris Nouveau Monde of the heSam consortium (Ph.D. grants to JH). This research was also funded by a NSERC strategic grant (STPGP336871) to FT. The research was carried out within the framework of the MONTABOR International Associated Laboratory (LIA France-Canada CNRS-UM2-EPHE-UQAT-UQAM-UQAC) and the European Union IRSES NewForest program.