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: