Statistical Analyses
To investigate the within- and among-population and genotype variation in phenotypic traits, each garden was modeled separately using linear random effects models fit by maximum likelihood in the lme4package in R (R Core Team 2014; Bates et al. 2015). The tree traits were modeled as response variables, while population and genotype nested within population were random effects. Garden plot was included as a random variable to help account for within-garden environmental variance. Statistical significance was calculated using likelihood ratio tests for the random effects using the package lmerTest(Kuznetsova et al. 2015).
To model population and genotype variation in phenotypic plasticity, we first needed to obtain replicated estimates of plasticity for each genotype and each trait. For each genotype with at least one tree in each of all three gardens, we randomly assigned all available trees into genotype triplets, with one tree from each garden. We then calculated plasticity for each triplet as the absolute value of the maximum difference in those three trait values. This produced a number of estimates of plasticity equal to the lowest number of trees available for that genotype in any garden. We repeated this random triplet assignment 100 times to obtain a set of possible plasticity datasets. For each dataset, we estimated the variance components necessary to calculate QST using a linear random effects model as described above, where trait plasticity was the response variable, and population and genotype nested within population were the random effects.
For each trait and trait plasticity, we compared the quantitative trait variation (QST) with genetic variance at neutral loci (FST). To calculate QST we used the following formula:
where σ2P is the additive genetic variance among populations and σ2G is the additive within-population variance (Spitze 1993; McKay & Latta 2002), i.e., the variance among genotypes within populations. Each trait or plasticity was analyzed using the models described above, and population and genotype variances were extracted to calculate QST. Parametric bootstrap and Bayesian estimation are considered the best methods to obtain a precision estimate around QST (O’Hara & Merilä 2005). We performed parametric bootstrapping to obtain a 95% confidence interval for QST, resampling the 16 populations with replacement 1000 times, and estimating QST for each bootstrapped data set. Resampling over the highest level in a hierarchical experimental design (here the population) is considered best practice (O’Hara & Merilä 2005). Variance in QST becomes quite large as the number of populations decreases (< 20), especially if populations are highly differentiated (O’Hara & Merilä 2005; Goudet & Büchi 2006). Goudet & Büchi (2006) recommend sampling many populations relative to the number of families. Our design of 16 populations with 12 genotypes per population comes close to their recommended sampling design of upwards of 20 populations with 10 families (O’Hara & Merilä 2005; Goudet & Büchi 2006). In using clonally replicated genotypes, our estimate of σ2G includes both additive and non-additive genetic effects, an approach that has been shown to lower QST estimates and is thus a conservative test of QST > FST (Cubry et al. 2017). Conversely, lower QST estimates derived from non-additive genetic effects contribute to a more liberal test of convergent selection (QST < FST) (Whitlock 2008; Cubry et al. 2017). To determine whether QST was significantly different from FST, we compared the 95% confidence intervals for both, which provides much stronger inference than simply comparing QST to the mean FST value (Whitlock 2008; Leinonen et al. 2013). Broad-sense heritability (H2) point estimates and confidence intervals were also calculated for each trait in each garden using the equation, H2 = σ2G/(σ2G +σ2E), where σ2E includes both the plot variance and the error variance. Calculations of heritability for plasticity did not include plot-level variance since plasticity was measured across gardens.
In order to test whether phenotypes showed strong climatic relationships, we regressed population trait means and trait plasticities against the first principal component (PC1) from the environmental PCA. We tested these regressions and calculated an adjusted R2 using a linear model in R (R Core Team 2014). Here, we used a single estimate of plasticity for each genotype. Using all available replicates for each genotype, we first calculated the mean trait value for each genotype in each garden, and then calculated plasticity as the maximum difference between gardens. Systematic differences among populations seen in these trait-climate correlations are another, stronger test for evidence of divergent selection acting over genetic drift (Whitlock 2008).
Results: