Temperature and Rainfall Are Separate Agents of Selection Shaping Population Di ﬀ erentiation in a Forest Tree

: Research highlights : We present evidence indicating that covariation of functional traits among populations of a forest tree is not due to genetic constraints, but rather selective covariance arising from local adaptation to di ﬀ erent facets of the climate, namely rainfall and temperature. Background and Aims : Traits frequently covary among natural populations. Such covariation can be caused by pleiotropy and / or linkage disequilibrium, but also may arise when the traits are genetically independent as a direct consequence of natural selection, drift, mutation and / or gene ﬂow. Of particular interest are cases of selective covariance, where natural selection directly generates among-population covariance in a set of genetically independent traits. We here studied the causes of population-level covariation in two key traits in the Australian tree Eucalyptus pauciﬂora. Materials and Methods : We studied covariation in seedling lignotuber size and vegetative juvenility using 37 populations sampled from throughout the geographic and ecological ranges of E. pauciﬂora on the island of Tasmania. We integrated evidence from multiple sources: (i) comparison of patterns of trait covariation within and among populations; (ii) climate-trait modelling using machine-learning algorithms; and (iii) selection analysis linking trait variation to ﬁeld growth in an arid environment. Results : We showed strong covariation among populations compared with the weak genetic correlation within populations for the focal traits. Population di ﬀ erentiation in these genetically independent traits was correlated with di ﬀ erent home-site climate variables (lignotuber size with temperature; vegetative juvenility with rainfall), which spatially covaried. The role of selection in shaping the population di ﬀ erentiation in lignotuber size was supported by its relationship with ﬁtness measured in the ﬁeld. Conclusions : Our study highlights the multi-trait nature of adaptation likely to occur as tree species respond to spatial and temporal changes in climate. eigenvector map (MEMs) variables (MEM 1–8) were additionally ﬁtted to account for spatial processes and unmeasured environmental variation and were derived using an inverse-weighted Euclidean distance matrix among provenance home-sites. MEM1 reﬂected a signiﬁcant positive latitudinal cline in similarity among provenances ( r = 0.78, p < 0.001); however, MEM2 to MEM8 were independent of the provenance home-site geographic co-ordinates and climate after Bonferroni adjustment for multiple testing (data not shown), and any importance is likely to reﬂect unmeasured environmental or historic facets of the variation. The importance of the independent variables was assessed through the percentage increase in mean-squared error (IncMSE) of the extended forests model, using the conditional permutation method of Strobl et al. (2008). Independent variables causing a decrease in IncMSE were sequentially dropped during model development. The lower and upper 95% conﬁdence intervals (red bar) associated with each independent variable were derived from iterative bootstrapping (1000 bootstraps) of the extended forests model.


Introduction
A knowledge of trait covariation is important to an understanding of phenotypic integration, modularisation, trade-offs and adaptive syndromes which may occur in natural populations [1][2][3][4][5]. Such trait covariation may be phenotypic, environmental or genetic, and occurs at multiple levels [2][3][4]6,7]. Of particular relevance is the genetic-based trait covariation which exists among populations of a species. This is particularly evident in forest trees, where information from decades of provenance trials has revealed substantial population differentiation in northern [8][9][10] and southern [11]

Common-Garden Experiments and Measurements
The present study used data from two common-garden trials located in a greenhouse and in the field, and comprising families (progeny based on open-pollinated seed collection from a single tree) from a sample of 37 native populations (hereafter also referred to as "provenances") of E. pauciflora, covering the full geographic and altitudinal ranges of the species on the island of Tasmania [30]. Apart from this range constraint and accessibility, these populations can be considered a random sample of the Tasmanian E. pauciflora. From the total of 275 open-pollinated families tested in the greenhouse [26], only 257 were used in the current study, as some families were not planted in the field trial. There was an average of 7 families per provenance. In both trials, the experimental layout was effectively a randomized complete block design, with 3 and 8 replicates for the greenhouse and field trials, respectively; the families were randomised within each replicate irrespective of their provenance, and were planted as single-tree plots. All provenances were represented across all replicates for both experiments, as were most of the families in the greenhouse trial (average number of trees per family = 3) and field trial (average number of trees per family was 5 for the field trial at the assessment age). Further details on the geographic location and climatic conditions of the provenance home-sites, sample collection of families within provenances, and establishment of the greenhouse and field trials are provided in Gauli et al. [26]. In brief, seedlings for both trials were initially germinated and grown as family lots in 40-cell seedling trays (BCC Hiko V93) in the same commercial nursery, but in different years. In each case, families were randomly arranged in the nursery with respect to provenance. For the greenhouse trial, 28 weeks after sowing, when seedlings had 2 to 3 nodes expanded, three seedlings from each family were pricked out into individual forestry tubes (120 mm × 50 mm). They were then arranged in their experimental design and moved to an unheated greenhouse, where they were maintained and measured 2 and 6 months later. Seedlings used for the field trial studied were from the same seed collection, but grown three years later. In this case, families were maintained in seedling trays until the experimental design was set up, just prior to field trial establishment, when seedlings were 16-to 17-months-old. The field trial was established on a geographically central (latitude 41.828 • S; longitude 147.138 • E), low altitude (185 m above sea level) site at the drier limits of the species distribution in Tasmania. The soil was ripped and mounded with rip lines 3 m apart. The trial was planted as a grid in August 2014, with the nursery-grown seedlings spaced 2.5 m apart within rip lines.
Two seedling traits measured in the greenhouse trial were considered in the current study, as they potentially have adaptive value in hot and dry climates. These functional traits represent facets of resource allocation (lignotuber size) and leaf ontogeny (development of the petiolate leaf). As with many plant species, eucalypts often exhibit dramatic ontogenetic variation in leaf morphology (i.e., they are heteroblastic). The vegetative transition from ontogenically juvenile to adult leaf forms is genetically independent of the onset of reproductive maturity and its timing may vary markedly between populations of a single species ( [25] and references therein). In the case of most eucalypts, including E. pauciflora, ontogenetically juvenile leaves are opposite and sessile, whereas mature/adult leaves are alternate and petiolate [26]. Leaf ontogenic variation in E. pauciflora was thus measured as the percentage of nodes expanded at 6 months of age which had petiolate leaves (PET; %), with the onset of early vegetative maturity (i.e., less vegetative juvenility) indicated by greater PET. Lignotuber size was measured as the lignotuber diameter at the cotyledonary node (LIGD; mm) [26]. Above-ground growth was used as a fitness surrogate, as measured by total tree height (HT; m) in the field trial 44 months after planting. The use of above-ground growth as a fitness surrogate has been suggested in long-lived species, such as forest trees, particularly at young stages when plants have not reached reproductive maturity [31,32]. Within eucalypt species, radial stem size (e.g., diameter or basal area) and tree height are usually highly positively correlated at the genetic and phenotypic levels [33] and for trees of the same cohort grown under competition, survival increases with radial stem size [34][35][36], consistent with size-dependent mortality [37]. Nevertheless, although HT growth 44 months after planting may be considered as a fitness surrogate in an early stage of a tree's life cycle, it will not necessarily represent a major component of lifetime fitness.

Data Analysis
Within-and among-provenance patterns of (co)variation were studied through linear mixed model analyses of greenhouse and field trial data for the two focal traits (LIGD and PET) and fitness surrogate (HT), respectively. We then explored selective agents likely acting on the focal traits, by using machine-learning algorithms to model the associations between provenance variation in LIGD and PET measured in the greenhouse trial and the climate of the provenance home-sites. Finally, to evaluate putative fitness consequences of variation in the focal traits, field HT growth was related to LIGD and PET via a regression framework for estimating selection. The following analyses were undertaken by using ASReml [38], SAS [39] and R [40] statistical software.

(Co)Variance Components and Genetic Parameters
Estimates of (co)variance components were obtained within the framework of the general linear mixed model: where y is the vector of observations, b and u are vectors of fixed and random effects (respectively), e is the vector of random residual terms, and X and Z are incidence matrices relating the observations to the model effects. The terms in b comprised fixed effects for the overall mean and replicates; in the analyses involving LIGD, the vector b also included a covariate (i.e., stem diameter measured at the cotyledonary node after 6 months' growth, hereafter denoted as STMD) to adjust the variation in the focal trait for differences among trees in radial growth. The vector u had sub-vectors for provenances and additive genetic effects within provenances; although initially included in u, interaction effects between provenances and replicates were subsequently dropped from the parameterization of a final model, as the related variance estimate was never found to be statistically significant (see hypothesis tests below). Univariate models were applied for the estimation of variance components for all characters, and a bivariate model was used for the estimation of covariance components between LIGD and PET (and again adjusting the variation in LIGD for the effects of STMD). In the Supplementary Materials (see Methods S1), details are provided for the definition of the (co)variance matrices associated with the random terms in u and e, under the univariate and bivariate models. The (co)variance components were estimated by restricted maximum likelihood (REML), using the average information REML algorithm ( [41]). Narrow-sense heritability estimates (ĥ 2 ) were then calculated for the fitness surrogate and focal functional traits as e are phenotypic, additive genetic, provenance and residual variance estimates, respectively. Correlation estimates between LIGD and PET were obtained according to standard formulae at the additive genetic, provenance and total genetic (i.e., additive plus provenance) levels. Standard errors for the heritability and correlation estimates were based on approximations using Taylor series expansion [42].
One-tailed [43] and two-tailed likelihood-ratio (LR) tests were pursued to evaluate whether the estimated variances and correlations differed significantly from zero, respectively, for additive genetic and provenance effects. A Wald F-test was conducted to determine whether differences among replicate fixed effects were statistically significant, with approximate denominator degrees of freedom being calculated as defined by Kenward and Roger [44]. The verification of model assumptions is described in the Supplementary Materials (see Methods S2).

Estimation of Marginal Means
Estimates of marginal means (hereafter denoted as "EMMeans") were obtained for a given character at the provenance and family levels, and to be used as input data for the climate-trait modelling and the regression framework for estimating selection (see below), respectively. This was pursued by fitting a linear fixed-effects model that included terms for replicate, provenance and family within provenance, as well as STMD used as a covariate in the model for LIGD. The variance of the estimated family EMMeans comprises both additive genetic and provenance effects (akin to the individual phenotypic variance defined in Equation (2)), in addition to the family-mean error variance.

Climate-Trait Associations
To explore selective agents likely acting on the focal traits, we modelled the provenance variation in LIGD and PET (i.e., the dependent variables) against the climate of the provenance home-sites. Climate data for each provenance home-site and the trial site were derived from contemporary climate surfaces , centred on 1990), obtained from the BIOCLIM module of the software package ANUCLIM (version 6.1 [45]) using site latitude, longitude and altitude (Table S1; Supplementary Materials). The provenance home-site climate data were the mean of the BIOCLIM data obtained for each mother tree sampled, and 19 bioclimatic variables (eleven temperature and eight precipitation variables) were chosen for analysis. We applied a principal component (PC) analysis on the correlation matrix of these variables to summarise the contemporary climate variation among provenance home-sites and to represent the climate differences under which provenances have been locally adapted [25]. Principal components with an eigenvalue greater than one were retained. To calculate the incremental change (i.e., climate anomaly) in the growing period (2014-2018) relative to the contemporary climate period at the trial site, daily minimum and maximum temperature and daily precipitation values were obtained from the Australian Bureau of Meteorology (http://www.bom.gov.au/jsp/awap/, accessed on 12 August 2018 [46]) using the AUSClim package (unpublished R package). Daily temperature and precipitation values were used to calculate the same bioclimatic variables (Table S1) at yearly steps, and then aggregated to obtain the growing period mean. The growing period climate mean was then subtracted from the contemporary climate mean calculated in AUSClim for each bioclimatic variable, to give the incremental change in climate at the trial site (Table S1).
For each focal trait, we related the provenance EMMeans with the provenance home-site climate represented by: (i) the main principal components obtained as described above; and (ii) the individual bioclimatic variables. This modelling was pursued using the nonparametric, machine-learning regression tree algorithm random forest [47] implemented through the extendedForest package of R. The extendedForest package extends the randomForest package [48] by incorporating the conditional permutation approach of Strobl et al. [49] to account for biased variable importance selection arising due to intercorrelation among predictor variables (see below). The two focal traits were modelled separately using the univariate extendedForest (http://gradientforest.r-forge.r-project.org/Conditionalimportance.pdf; hereafter denoted as "extended forests model"), and jointly by using the multivariate random forest approach implemented in the gradientForest (hereafter denoted as "gradient forests model"; [50]) package of R. Following the approach outlined by Fitzgerald and Keller [20], we accounted for spatial processes and unmeasured environmental variation (i.e., latent variables) using the Moran's eigenvector map variables (MEMs) [51], with 50% of the non-negative MEMs included in the regression trees. Positive MEMs describe the broad-scale 'global' spatial structure in the provenance level data [51]. In our case, eight non-negative MEMs were retained, none of which were significantly correlated with the climate-derived principal components nor with the individual contemporary climate variables (data not shown). The MEMs were estimated through eigen decomposition of the inverse-distance weighted Euclidean matrix calculated from a Gabriel's graph, derived by converted latitude and longitude coordinates of each provenance home-site into UTM mapping units using the adspatial package of R [52].
The importance of each independent variable was assessed through the percentage increase in the mean-squared error (IncMSE) of the extended forests model using a conditional permutation approach [49]. This approach proceeds for each tree within the forest by firstly estimating the out-of-bag (OOB, which concerns the data not contained within the bootstrapped subsample) prediction accuracy following Breiman [47]. Under this procedure, the subset of independent variables that are correlated with a given independent variable X i with a Pearson correlation equal to or greater than 0.5 are identified, and an expanded grid that bisects the sample space is generated using the cut-points that split each correlated variable with X i in a given tree. Within this expanded grid, the values of X i are permuted and the OOB prediction accuracy is estimated. The difference between the prediction accuracy before and after permutation is thus the importance of X i for a given tree. The average permutated importance over all trees is the overall importance of X i in predicting the dependent variable. Independent variables contributing to a decrease in IncMSE (indicative of low importance) were sequentially dropped during model development. Partial dependency plots from the final model were obtained to visualise the effect of an independent variable on the dependent variable whilst averaging the effect of all other independent variables. Iterative bootstrapping (1000 bootstraps) of the extended forests model was applied to calculate 95% confidence intervals for: (i) the predictive ability of the model (pseudo-r 2 ); (ii) the conditional permutated variable importance; and (iii) the marginal effect (depicted in the partial dependency plots) of an independent variable on a dependent variable, following the approach detailed in Costa e Silva et al. [25].

Selection Analysis
A regression framework for estimating selection was applied to determine whether the selection patterns expected for the two focal traits were consistent with the directionality of provenance transfer regarding the difference between the provenance home-site contemporary climate and (i) the trial site contemporary climate; and (ii) the growing period climate (i.e., climate anomaly). This selection analysis was performed by modelling the fitness consequences of variation in the focal traits (LIGD and PET), considering HT as a fitness surrogate. Assessment of the form, strength and direction of selection followed analytical procedures similar to those described by Lande and Arnold [21] except that, instead of data consisting of individual phenotypic values, family EMMeans were used as observation units in the selection analysis [53,54]. Furthermore, these family EMMeans encapsulate variation both within and between provenances and thus, the selection analysis is not restricted to a particular provenance but it refers to an extended inference space that better represents the range of phenotypic variation found across the natural distribution of the species (see the Discussion Section).
Prior to modelling the fitness surrogate as a function of the focal traits, data were pre-processed as follows (e.g., as suggested by Lande and Arnold [21]): for HT, a given family EMMean was divided by the overall trial mean to obtain a relative fitness value (with mean equal to one); for each focal trait, the family EMMeans were converted to have a mean of zero and a variance equal to one. After this centring and standardization, the squares and cross-product of the input independent variables were obtained for the estimation of non-linear selection gradients (see below); in addition, to estimate appropriate quadratic selection gradients, squared trait observations were halved before modelling (see also [55]).
We first obtained standardized directional selection differentials for each focal trait to evaluate the total (i.e., direct and indirect) effects of linear selection, estimated by using a single-trait regression model relating relative fitness (dependent variable) and a variance-standardized trait (independent variable). In this context, estimates of linear regression coefficients are equal to the covariances of relative fitness with variance-standardized trait values and thus, refer to standardized directional selection differentials [56].
Selection gradients are more appropriate than selection differentials for measuring the direct consequences of variation in a trait's value on relative fitness. We estimated standardized selection gradients by using a multiple-trait regression model relating relative fitness with both of the variance-standardized focal traits. Following the suggestion of Lande and Arnold [21], a first-order polynomial model was used to estimate standardized linear selection gradients and then, a separate second-order polynomial model (including also a cross-product term) was applied to estimate standardized non-linear selection gradients (for further details, see also Methods S3; Supplementary Materials). The complete multiple-trait regression model with linear, quadratic and cross-product terms provides the best quadratic approximation to the fitness surface [22]. The vector of standardized linear selection gradients (β) determines the slope of the fitness surface, and it measures the strength of positive (β > 0) or negative (β < 0) directional selection on a trait; for a variance-standardized trait,β quantifies the change in relative fitness due to increasing the trait value by one standard deviation unit from its mean. The matrix of standardized non-linear selection gradients (γ) has information regarding the curvature and orientation of the best quadratic approximation to the fitness surface, with the quadratic estimates in the diagonal elements pertaining to convex (γ ii < 0) or concave (γ ii > 0) forms of selection on an individual trait, and the cross-product estimate in an off-diagonal element describing correlational selection on a pairwise combination of traits (i.e.,γ ij > 0 andγ ij < 0 indicate that traits i and j are selected to become positively and negatively correlated, respectively) [21,22]. Statistically significant convex and concave forms of quadratic selection can be interpreted as actual stabilizing and disruptive selection, respectively, only when relative fitness has a maximum or a minimum at a stationary point corresponding to a trait value located at an intermediate position within the range of the observations [22,57]. Diagnostic measures indicated that the fitted single-and multiple-trait regression models were adequate for parameter estimation and related statistical hypothesis testing (see Methods S2; Supplementary Materials).

(Co)Variance Components and Genetic Parameters
Estimates of overall means and standard deviations are given Table S2 (Supplementary Materials) for the functional traits and fitness surrogate. Table 1 presents the additive genetic and provenance variances, along with narrow-sense heritabilities, estimated for the focal functional traits and fitness surrogate. Differences among the fixed replicate effects were significant at the 5% level of significance for HT in the field trial, but not for LIGD and PET in the greenhouse trial. In general, highly significant (p < 0.001) variance estimates were detected for all characters at the additive genetic and provenance levels; both focal traits and the fitness surrogate appeared to be under moderate additive genetic control, with narrow-sense heritabilities varying from 0.203 ± 0.085 for LIGD to 0.425 ± 0.127 for PET (Table 1). = narrow-sense heritability. Estimates of variance components and heritabilities (± approximate standard errors) are presented for: LIGD = lignotuber diameter at the cotyledonary node; PET = percentage of nodes expanded at six months of age which had petiolate leaves; HT = total tree height at age 44 months. In parenthesis, the significance probability obtained from a one-tailed likelihood-ratio test is given for the estimated variance components.
The estimated correlations between the focal traits at the additive genetic, provenance and total genetic (i.e., additive plus provenance) levels were −0.21 ± 0.25, −0.69 ± 0.14 and −0.37 ± 0.14, respectively, with the consistent negative sign indicating that large lignotuber size is associated with greater vegetative juvenility. However, the additive genetic correlation estimate was not significantly different from zero (p > 0.05), a result that does not provide support for a putative pleiotropic relationship between LIGD and PET, and indicates that the two traits are not likely to be co-dependent. Conversely, the estimated provenance correlation was highly significant (p < 0.001), hence was the component that had the most important contribution to the total genetic correlation estimate. The significant provenance correlation was partly reflected in the tendency for provenances in the central regions of the distribution to have larger LIGD and below average PET ( Figure 1). additive genetic control, with narrow-sense heritabilities varying from 0.203 ± 0.085 for LIGD to 0.425 ± 0.127 for PET (Table 1). The estimated correlations between the focal traits at the additive genetic, provenance and total genetic (i.e., additive plus provenance) levels were −0.21 ± 0.25, −0.69 ± 0.14 and −0.37 ± 0.14, respectively, with the consistent negative sign indicating that large lignotuber size is associated with greater vegetative juvenility. However, the additive genetic correlation estimate was not significantly different from zero (p > 0.05), a result that does not provide support for a putative pleiotropic relationship between LIGD and PET, and indicates that the two traits are not likely to be co-dependent. Conversely, the estimated provenance correlation was highly significant (p < 0.001), hence was the component that had the most important contribution to the total genetic correlation estimate. The significant provenance correlation was partly reflected in the tendency for provenances in the central regions of the distribution to have larger LIGD and below average PET ( Figure 1).  . The provenance estimates of marginal means (EMMeans) for LIGD are plotted such that the larger the brown circle the greater the LIGD mean is above the grand mean, and the larger the green triangle the more the mean is below the grand mean. For PET, the plotting is reversed such that the smaller values of the EMMeans for PET (thus plants are more juvenile) are depicted with larger brown circles. Seedlings from provenances in the central region of the distribution tend to have larger lignotubers and to retain juvenility longer (i.e., low proportion of nodes with petiolate leaves).

Climate-Trait Associations
The variable importance and partial dependency plots from the final extended forests models using four climate principal components (PC1-PC4, accounting for 95% of the climate variation; Table S1) and eight spatial MEMs to explain the variation among provenance EMMeans for LIGD and PET are shown in Figures 2 and 3, respectively. These models explained 76% of the variation in LIGD, but only 37% of the variation in PET. The first climate principal component (PC1) was the most important predictor of LIGD and PET ( Figure 2). PC1 accounted for 54% of the total variation in the home-site contemporary climate among the sampled provenances; it was negatively correlated with altitude (r = −0.77; p < 0.001), and reflected an "aridity" gradient with positive values along PC1 associated with increasing temperature and decreasing rainfall (Table S1). The major changes in the partial dependency plots were supported by 90% of the data and showed that LIGD increased ( Figure 3a) and PET decreased (Figure 3b) with increasing values along PC1 (increasing aridity); however, the bootstrapped 95% confidence interval was larger for the predictions of PET at the positive extreme of PC1 (Figure 3b).
Forests 2020, 11, x FOR PEER REVIEW 9 of 22 seedlings from 37 Eucalyptus pauciflora provenances distributed within the island of Tasmania (reanalysed data from Gauli et al. (2015). The provenance estimates of marginal means (EMMeans) for LIGD are plotted such that the larger the brown circle the greater the LIGD mean is above the grand mean, and the larger the green triangle the more the mean is below the grand mean. For PET, the plotting is reversed such that the smaller values of the EMMeans for PET (thus plants are more juvenile) are depicted with larger brown circles. Seedlings from provenances in the central region of the distribution tend to have larger lignotubers and to retain juvenility longer (i.e., low proportion of nodes with petiolate leaves).

Climate-Trait Associations
The variable importance and partial dependency plots from the final extended forests models using four climate principal components (PC1-PC4, accounting for 95% of the climate variation; Table  S1) and eight spatial MEMs to explain the variation among provenance EMMeans for LIGD and PET are shown in Figures 2 and 3, respectively. These models explained 76% of the variation in LIGD, but only 37% of the variation in PET. The first climate principal component (PC1) was the most important predictor of LIGD and PET (Figure 2). PC1 accounted for 54% of the total variation in the home-site contemporary climate among the sampled provenances; it was negatively correlated with altitude (r = −0.77; p < 0.001), and reflected an "aridity" gradient with positive values along PC1 associated with increasing temperature and decreasing rainfall (Table S1). The major changes in the partial dependency plots were supported by 90% of the data and showed that LIGD increased (Figure 3a) and PET decreased (Figure 3b) with increasing values along PC1 (increasing aridity); however, the bootstrapped 95% confidence interval was larger for the predictions of PET at the positive extreme of PC1 (Figure 3b).  Table S1, which accounted for 95% of the climate variation among  Table S1, which accounted for 95% of the climate variation among provenances. Moran's eigenvector map (MEMs) variables (MEM 1-8) were additionally fitted to account for spatial processes and unmeasured environmental variation and were derived using an inverse-weighted Euclidean distance matrix among provenance home-sites. MEM1 reflected a significant positive latitudinal cline in similarity among provenances (r = 0.78, p < 0.001); however, MEM2 to MEM8 were independent of the provenance home-site geographic co-ordinates and climate after Bonferroni adjustment for multiple testing (data not shown), and any importance is likely to reflect unmeasured environmental or historic facets of the variation. The importance of the independent variables was assessed through the percentage increase in mean-squared error (IncMSE) of the extended forests model, using the conditional permutation method of Strobl et al. (2008). Independent variables causing a decrease in IncMSE were sequentially dropped during model development. The lower and upper 95% confidence intervals (red bar) associated with each independent variable were derived from iterative bootstrapping (1000 bootstraps) of the extended forests model. reflect unmeasured environmental or historic facets of the variation. The importance of the independent variables was assessed through the percentage increase in mean-squared error (IncMSE) of the extended forests model, using the conditional permutation method of Strobl et al. (2008). Independent variables causing a decrease in IncMSE were sequentially dropped during model development. The lower and upper 95% confidence intervals (red bar) associated with each independent variable were derived from iterative bootstrapping (1000 bootstraps) of the extended forests model.  Table S1) on the provenance estimates of marginal means for: (a) lignotuber diameter at the cotyledonary node (LIGD); and (b) percentage of nodes expanded at six months of age which had petiolate leaves (PET; lower values of PET pertain to seedling leaves that are more juvenile). The red tick marks show the deciles of the data. The grey shading indicates the 95% confidence intervals derived by iterative bootstrapping (1000 bootstraps). Increasing scores on PC1 corresponds to increasing home-site aridity (i.e., high summer temperatures and low summer rainfall; see Table S1). The images show the development of the basal lignotuber in the top panel (yellow arrow), and the ontogenetic change from opposite, sessile (bottom) to alternate, petiolate (top) leaves in the lower panel. The yellow arrow at the bottom panel corresponds to the node of the first petiolate leaf (node 4, left; node 6, right).  Table S1) on the provenance estimates of marginal means for: (a) lignotuber diameter at the cotyledonary node (LIGD); and (b) percentage of nodes expanded at six months of age which had petiolate leaves (PET; lower values of PET pertain to seedling leaves that are more juvenile). The red tick marks show the deciles of the data. The grey shading indicates the 95% confidence intervals derived by iterative bootstrapping (1000 bootstraps). Increasing scores on PC1 corresponds to increasing home-site aridity (i.e., high summer temperatures and low summer rainfall; see Table S1). The images show the development of the basal lignotuber in the top panel (yellow arrow), and the ontogenetic change from opposite, sessile (bottom) to alternate, petiolate (top) leaves in the lower panel. The yellow arrow at the bottom panel corresponds to the node of the first petiolate leaf (node 4, left; node 6, right).
Consistently with the observed negative sign of the estimated provenance correlation, LIGD and PET showed an increasing and decreasing response to increasing home-site aridity (positive values along PC1; Figure 3), respectively. However, when fitting the extended forests model as a function of specific climate variables and the eight MEMs, there was support for the variation among the provenance EMMeans to be associated with different facets of the climatic variation described by PC1 (Figure 4). In this sense, temperature variables were among the most important for predicting provenance variation in LIGD, particularly maximum temperature of the warmest week (TMXWW), compared with the MEMs that were among the least important (Figure 4a). In contrast, the extended forests model fitted for PET indicated that precipitation variables (in particular the precipitation of the driest quarter; RDRYQ) and MEMs tended to be more important than temperature variables in explaining the provenance variation in EMMeans (Figure 4b). Joint modelling with gradientForest ( Figure S1; Supplementary Materials) confirmed the importance of TMXWW for LIGD and RDRYQ for PET, but suggested greater independence of PET to changes in home-site TMXWW than implied by the extended forest model ( Figure S1a vs. Figure 5c). compared with the MEMs that were among the least important (Figure 4a). In contrast, the extended forests model fitted for PET indicated that precipitation variables (in particular the precipitation of the driest quarter; RDRYQ) and MEMs tended to be more important than temperature variables in explaining the provenance variation in EMMeans (Figure 4b). Joint modelling with gradientForest ( Figure S1; Supplementary Materials) confirmed the importance of TMXWW for LIGD and RDRYQ for PET, but suggested greater independence of PET to changes in home-site TMXWW than implied by the extended forest model ( Figure S1a versus Figure 5c).  Table S1. Independent variables causing a decrease in IncMSE were sequentially dropped during model development. The most important independent variables in the models for LIGD and PET were the provenance home-site maximum temperature of the warmest week (TMXWW) and the precipitation of the driest quarter (RDRYQ), respectively. The MEMs are as detailed in the Figure 2 caption. Figure 5 shows the partial dependency plots from the extended forests model using the home-site climate variable that was found to be the most important in explaining the provenance variation for a given trait (i.e., TMXWW for LIGD, and RDRYQ for PET). The main change in provenance EMMeans for PET occurred at the lower rainfall end of the RDRYQ gradient, with provenances originating from sites receiving less than 150 mm precipitation during the driest quarter tending to have a greater vegetative juvenility (decreasing PET; Figure 5d). However, this change in RDRYQ had little effect on  Table S1. Independent variables causing a decrease in IncMSE were sequentially dropped during model development. The most important independent variables in the models for LIGD and PET were the provenance home-site maximum temperature of the warmest week (TMXWW) and the precipitation of the driest quarter (RDRYQ), respectively. The MEMs are as detailed in the Figure 2 caption. Figure 5 shows the partial dependency plots from the extended forests model using the home-site climate variable that was found to be the most important in explaining the provenance variation for a given trait (i.e., TMXWW for LIGD, and RDRYQ for PET). The main change in provenance EMMeans for PET occurred at the lower rainfall end of the RDRYQ gradient, with provenances originating from sites receiving less than 150 mm precipitation during the driest quarter tending to have a greater vegetative juvenility (decreasing PET; Figure 5d). However, this change in RDRYQ had little effect on LIGD (Figure 5b). Rather, the change in provenance EMMeans for LIGD was most aligned with a positive monotonic response to increasing provenance home-site TMXWW (Figure 5a). While the extended forests model showed a small monotonic decrease in PET with increasing TMXWW above 19 • C (Figure 5c), this climate variable was only of moderate importance in the fitted model (Figure 4b).
The trial site was located at the hot, dry extreme of the aridity gradient for the natural distribution of E. pauciflora in Tasmania ( Figure 6). Over the growing period assessed (2014-2018), TMXWW was estimated to be 1.1 • C warmer and RDRYQ to be 9 mm less than the contemporary average  for the trial site (see Table S1). The TMXWW of the growing period was also, on average, 1.5 • C warmer than the contemporary home-site climate for any of the provenances translocated onto the trial site. In contrast, the growing period was less extreme for RDRYQ, with seven of the 37 provenances originating from home-sites with average contemporary RDRYQ even drier than the growing period (i.e., > 111 mm; Figure 6). Accordingly, if RDRYQ and TMXWW are indeed the main climatic drivers of provenance differentiation in PET and LIGD, respectively, a greater proportion of the gene pool at the trial site would be expected to experience directional selection for larger LIGD rather than for lower PET.  (Figure 5b). Rather, the change in provenance EMMeans for LIGD was most aligned with a positive monotonic response to increasing provenance home-site TMXWW (Figure 5a). While the extended forests model showed a small monotonic decrease in PET with increasing TMXWW above 19 °C (Figure 5c), this climate variable was only of moderate importance in the fitted model (Figure 4b). The trial site was located at the hot, dry extreme of the aridity gradient for the natural distribution of E. pauciflora in Tasmania ( Figure 6). Over the growing period assessed (2014-2018), TMXWW was estimated to be 1.1 °C warmer and RDRYQ to be 9 mm less than the contemporary average  for the trial site (see Table S1). The TMXWW of the growing period was also, on average, 1.5 °C warmer than the contemporary home-site climate for any of the provenances translocated onto the trial site. In contrast, the growing period was less extreme for RDRYQ, with seven of the 37 provenances originating from home-sites with average contemporary RDRYQ even drier than the growing period (i.e., > 111 mm; Figure 6). Accordingly, if RDRYQ and TMXWW are indeed the main climatic drivers of provenance differentiation in PET and LIGD, respectively, a greater proportion of the gene pool at the trial site would be expected to experience directional selection for larger LIGD rather than for lower PET. The importance of the climate variables in the different models is given in Figure 4.   Table S1), such that mesic provenances are dark blue (more negative values on PC1) and arid provenances are dark red (more positive values on PC1). The TMXWW and RDRYQ values for the contemporary (LTA: 1976-2005, centred on 1990) and growing (GP: 2014-2018) climate periods at the trial site are indicated. All provenances are from home-sites with a contemporary climate cooler than the growing period TMXWW, but seven provenances originate from home-sites with a contemporary climate drier than the growing period RDRYQ.

Selection Differentials and Gradients
The standardized selection differentials (representing the total-direct and indirect-effects of selection) and gradients (describing the direct effects of selection) estimated for the focal traits are shown in Table 2. We found highly significant (p < 0.001) standardized linear selection differentials for Figure 6. The provenance home-site maximum temperature of the warmest week (TMXWW) plotted against the precipitation of the driest quarter (RDRYQ). Increasing circle size pertains to increasing provenance home-site altitude (meters above sea level). The colour of each circle reflects the provenance score along the climate principal component 1 (PC1; see Table S1), such that mesic provenances are dark blue (more negative values on PC1) and arid provenances are dark red (more positive values on PC1). The TMXWW and RDRYQ values for the contemporary (LTA: 1976-2005, centred on 1990) and growing (GP: 2014-2018) climate periods at the trial site are indicated. All provenances are from home-sites with a contemporary climate cooler than the growing period TMXWW, but seven provenances originate from home-sites with a contemporary climate drier than the growing period RDRYQ.

Selection Differentials and Gradients
The standardized selection differentials (representing the total-direct and indirect-effects of selection) and gradients (describing the direct effects of selection) estimated for the focal traits are shown in Table 2. We found highly significant (p < 0.001) standardized linear selection differentials for both traits and suggesting total directional selection favouring increased LIGD but reduced PET: the selection estimates indicated that an increase of 14% and a decrease of 6% in relative fitness would result from increasing the values of LIGD and PET by one standard deviation unit from their mean, respectively. For LIGD, the estimated standardized linear selection gradient (β) was highly significant and positive, hence indicating direct directional selection for larger LIGD; in addition, for this trait,β was similar in magnitude to the estimated linear selection differential. Conversely, theβ estimate for PET was not statistically significant (p > 0.05) and, albeit also negative in sign, its magnitude was much smaller in absolute value than the corresponding linear selection differential.
In other words, when compared with the correspondingβ estimate, the more negative estimate of the linear selection differential for PET may reflect the effects of: (i) a strong direct directional selection for increased LIGD; and (ii) a negative correlation of LIGD with PET (i.e., a statistically significant family-mean correlation of −0.32, which compares well with the total genetic correlation of −0.37). This suggests that the estimated linear selection differential for PET is mainly driven by the effect of indirect directional selection operating through LIGD (and possibly also via unmeasured traits) rather than direct directional selection on PET per se (see Methods S3 in Supplementary Materials for further details). Table 2. Standardized linear selection differentials, and standardized linear and non-linear selection gradients, estimated for the functional traits measured in E. pauciflora. Functional traits: LIGD = lignotuber diameter at the cotyledonary node; PET = percentage of nodes expanded at six months of age which had petiolate leaves. The regression framework for estimating selection used relativized fitness values (i.e., scaled to have a mean of one), and trait observations that have been centred (mean of zero) and standardized (unit variance). Standardized selection differentials refer to linear (directional) selection estimates.

Linear Selection Differentials Linear and Non-Linear Selection Gradients
Standardized selection gradients pertain to: β = vector containing estimates (β) of linear (directional) selection gradients for the two focal traits; γ = matrix of non-linear selection gradients, with the diagonal elements (γ ii ) referring to estimates of quadratic selection for a given trait, and the off-diagonal element (γ ij ) corresponding to an estimate of correlational selection between the two traits. Estimates of standardized selection differentials and selection gradients are provided with their standard errors, as well as significance probabilities in parenthesis.
The estimated standardized non-linear selection gradients did not support statistically significant direct effects of quadratic selection (γ ii ) for either LIGD or PET; in addition, as indicated by a non-significant correlational selection estimate (γ ij ), the interaction between LIGD and PET did not appear to determine relative fitness (although the negative sign of the estimate would have been indicative of selection against a positive covariance between the focal traits). Thus, a non-significant γ ij indicates that the effect of one trait on fitness does not depend on the state of the other trait, suggesting that LIGD and PET are unlikely to be functionally linked (i.e., they are not interdependent traits-sensu [7]).

Discussion
Although genetic correlations due to pleiotropy and/or linkage disequilibrium can cause multiple traits to evolve jointly, among-population covariance can also be generated between genetically independent traits [13]. Significant within-population additive genetic variance was detected for the two focal traits studied but no significant additive genetic correlation was found between them and thus, there is no evidence to suggest that these traits are co-dependent. Therefore, the additive genetic covariance between LIGD and PET is unlikely to have contributed to their strong (co)variation at the population level. A caveat here is that this result assumes a single, common additive genetic correlation between the two studied traits for all populations. While additive genetic (co)variances may themselves evolve and be heterogeneous across populations [58,59], limited sample sizes (i.e., number of families per provenance) precluded obtaining accurate estimates for each provenance in the present study. Nevertheless, apart from a suggestion that early lignotuber development is a juvenile characteristic [60], as is the absence of petiole development (i.e., low PET values), we have no physiological and developmental reason to suspect a pleiotropic relationship between LIGD and PET. This conclusion is consistent with the absence of co-located QTL reported for lignotuber size [61] and vegetative juvenility traits [62] in Eucalyptus globulus. Thus, rather than being a covariance generated by divergence of genetically correlated traits, the strong among-population covariance between LIGD and PET is likely a direct consequence of selective covariance [13,14].
Armbruster and Schwaegerle [13] proposed the following possible scenarios leading to selective covariance: (i) traits being functionally linked, whereby the fitness consequences of one trait are affected by the state of the other trait (i.e., interdependent traits-sensu [7]); (ii) geographical variation in an agent of selection that affects both traits similarly; and (iii) geographic covariation in different agents of selection that independently affect the traits. As our fitness surrogate did not appear to be significantly determined by the interaction involving LIGD and PET (Table 2), there is no evidence for correlational selection between these traits (i above). This result suggests that the two focal traits independently affect fitness and thus, are not likely to be functionally linked. Therefore, under the apparent absence of pleiotropic and functional relationships between LIGD and PET (i.e., they are co-specialised traits-sensu [7]), their covariation at the among-population level could potentially be explained by response to the same (ii above) or different (iii above) agents of selection.
While the strong covariation detected at the among-provenance level for LIGD and PET may reflect a classic case of selective covariance driven by adaptation to one or more environmental factors, the possibility of the potential contributions to among-population covariance from genetic drift, mutation and/or gene flow cannot be ignored [7,13,14]. However, for the studied E. pauciflora provenances, there are multiple lines of evidence indicating that the contributions from these latter processes may not be important when compared with the portion of the among-population covariance that may have been generated by the direct effects of natural selection (i.e., selective covariance). First, as noted previously, the Q ST for both traits significantly exceeded neutral expectations, using the maximum F ST estimated from various microsatellite (SSR) markers [26], signalling that provenance differentiation is primarily the result of divergent selection acting across the gene pool of E. pauciflora in Tasmania. Secondly, the present study showed a significant effect of variation in climate on provenance variation in both traits, particularly a major gradient in aridity involving increasing temperatures and decreasing rainfall. The role of these climate variables as agents of selection is supported by reports of similar climate-trait associations in other eucalypt species for both lignotuber size (e.g., E. viminalis [63,64]; E. obliqua [65]) and vegetative juvenility (E. risdonii/tenuiramis [25]). Finally, our selection analysis provided evidence for direct directional selection, at least for LIGD as indicated by a significant favourable effect of increased lignotuber size on the fitness surrogate measured in the arid test environment. While we did not find evidence for direct directional selection for PET, this was consistent with its weaker association with climate when compared with LIGD. It was also consistent with the growing period of our trial site being less deviant, in terms of the home-site of the translocated provenances, for the climate variable (i.e., see precipitation of the driest quarter in Figure 6) that best explained the provenance differentiation in PET.
The regression framework for estimating selection involved studying the relationship between the focal traits measured under greenhouse conditions and a fitness surrogate evaluated on related individuals grown in hot and dry, low-altitude field conditions. Estimation of genetic parameters based on experiments using individuals measured under laboratory/greenhouse and nature/field conditions has also been pursued in other studies for focal traits with unimportant genotype by environment interaction (e.g., [13,66]). In the present case, this assumption is consistent with previous work with other eucalypt species which indicates that while environmental factors may affect the development of lignotubers [67,68] and retention of vegetative juvenility [69,70], the expression of provenance or additive genetic differences is relatively stable across field environments [64,[71][72][73] or experimental treatments [63,70,74,75]. Accordingly, it is reasonable to assume that genotype by environment interactions are not important for the two focal traits studied, and that (apart from scale effects) trait assessments undertaken in the greenhouse and field common-garden trials can be regarded as effectively the same character ( [76], Chapter 19). In addition, while selection estimates are usually obtained for individual populations (e.g., [15,77]), we pooled data from the sampled provenances to adequately represent the range of phenotypic variation found for the species in Tasmania and to enhance statistical power for detecting selection. This approach is similar to the "site-specific" selection analysis pursued by Etterson [78], providing a picture of the overall selection by not restricting the phenotypic variation to the distributional range within each population (see also [79,80]).
A key finding of the present study is the evidence from climate-trait modelling that different, but geographically covarying, causal agents of selection may have shaped the patterns of variation in provenance means for LIGD and PET. Provenance variation in both traits was reasonably well explained by an altitude-related gradient in aridity (i.e., PC1, Table S1; Figures 2 and 3) that occurs across the range of E. pauciflora in Tasmania, which would superficially suggest that selective covariance is due to both traits responding to the same agent of selection. However, this does not appear to be the case. As often occurs with such major gradients, the environmental change is complex and may involve parallel changes in multiple biotic and abiotic factors [12,81,82]. We found that the major climate gradient explaining the provenance differentiation was mainly characterized by a decrease in rainfall and an increase in temperature, a trend often reported in studies of plant adaptation to gradients in altitude [83,84], latitude [85] or aridity/continentality [25,86]. However, when examining the effect of specific climate variables, the variation among provenance means for a given trait appeared to be explained by different facets of the climatic variation. LIGD was better explained by temperature than rainfall variables, and specifically by maximum temperature of the warmest week (TMXWW). In contrast, the provenance differentiation in PET was poorly explained by temperature variables and was better related to rainfall variables, specifically to precipitation of the driest quarter (RDRYQ). Such differential partitioning of the adaptive response of plants to temperature and rainfall variables has been similarly noted at the genomic level [87][88][89], and is likely a common feature of climate adaptation. In our study, the significant negative among-population correlation between LIGD and PET appears to be mainly due to the increase in LIGD associated with the rise in home-site TMXWW from 21 to 24 • C, and the decrease in PET at the dry extreme of the RDRYQ gradient as home-sites for this climate variable were reduced from 150 to 110 mm (Figure 5a,d and Figure S1).
Lignotubers occur in many eucalypt species [90,91], are often evident early in seedling development [92], and have been shown to increase in frequency in provenances with high home-site temperatures [26,93] or low rainfall [65]. These basal swellings comprise protected vegetative buds, vascular tissue and carbohydrate reserves, and allow plants to regenerate after stem death or damage [65,94]. Non-structural carbohydrates, mainly in the form of starch, are mobilised during resprouting when photosynthates are otherwise limited [65,95]. However, while a selective advantage associated with resprouting could lead to large lignotubers being favoured in hot and dry environments, this may not hold in the present case as few plants were recorded as having stem damage after planting in the field trial (< 4%). While we cannot dismiss the possibility that lignotubers may impact plant-water relationships [96], in the absence of stem damage, the positive consequences of variation in LIGD on the fitness surrogate (i.e., above-ground growth) ( Table 2) could be related to the remobilisation of stored non-structure carbohydrates to enhance recovery and rebuilding of photosynthetic capacity following stress events [97], such as defoliation [98][99][100], short-periods of drought or high temperatures [98,101,102]. However, selection analysis using a multiple-trait regression approach does not preclude the possibility that selection may occur for an observed trait via indirect selection acting on unmeasured correlated traits ( [21,57]; see also Methods S3).
The extent to which climate explained the provenance differentiation in PET was less than that observed for LIGD, with unexplained spatial and environmental factors (MEMs) more important for PET. In addition, we did not find evidence for direct directional selection for PET. Gene flow is one possible cause of covariance of traits among populations [13,14] and, in this sense, we cannot dismiss the possibility that the climate signal in PET has been reduced by interspecific gene flow, with introgression from a species with delayed vegetative phase change contributing to population differentiation. In the present case, the lowland related species E. tenuiramis is a candidate, as it does hybridise with E. pauciflora when they co-occur in south-eastern Tasmania [103]. However, if such introgression involved both of the focal traits and contributed to their covariance at the among-population level, significant within-population covariance would likely to have been also detected, at least in the early phase of the introgression [104]. PET increased non-linearly with rainfall of the driest quarter in E. pauciflora. Although sessile and opposite juvenile leaves are only a transient developmental feature of E. pauciflora, the trend for a greater retention of the juvenile leaf form in provenances from drier sites was consistent with that observed by Costa e Silva et al. [25] in a related species complex which has more persistent juvenile foliage. These authors suggested that the delayed development of adult leaves on drier sites could be due to the susceptibility of the petiole to cavitation [105,106], and that early juvenile-adult leaf transition favoured on more mesic sites would allow leaves to be more readily orientated to maximise light interception for photosynthesis [107].

Conclusions
Interpreting the origin of covarying patterns of trait variation in nature is complex and a well recognised challenge in evolutionary biology (e.g., [2,7,13]). However, this is important as such patterns can provide insights into the evolutionary processes and selection agents likely to shape the adaptive response of populations in the future. In the case of the Tasmanian populations of E. pauciflora, we used multiple lines of evidence to argue that provenance covariation for seedling lignotuber size and retention of vegetative juvenility is due to selective covariance arising from adaptation to different facets of the climate, viz. temperature and rainfall variables, respectively, which geographically covary. Our evidence is derived from the comparison of patterns of trait covariation within and among populations, climate-trait modelling and selection analysis, each of which has assumptions and alternative interpretations which need to be considered. Nevertheless, we highlight the multi-trait nature of adaptation likely to occur when tree species experience spatial and temporal changes in climate, which may involve different climate variables acting as agents of selection causing a correlated change among populations in a set of plant traits. This scenario is likely to be broadly applicable to other forest tree species which occur over complex environmental gradients, and we emphasise the need to dissect the environmental drivers and trade-offs among functional traits involved in contemporary adaptation to better predict species and population responses to climate change.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/12/1145/s1, Methods S1: Variance-covariance matrices for random terms defined under the linear mixed model. Methods S2: Verification of model assumptions ((Co)variance components and genetic parameters; Selection analysis). Methods S3: Notes on the interpretation of selection differentials and gradients. Figure S1: Importance of the independent variables in the extended forests model for (a) lignotuber diameter at the cotyledonary node (LIGD), and (b) percentage of nodes expanded at six months of age which had petiolate leaves (PET). Table S1: Bioclimatic variables used in a principal component (PC) analysis to summarise the variation in home-site contemporary climate for the provenances translocated to the studied trial site, as well as to derive the contemporary and growing period climates for this site. Table S2: Summary statistics for the focal functional traits and fitness surrogate measured in E. pauciflora.
Author Contributions: The reported study was conceived by J.C.S., B.P. and P.A.H., T.B., B.P. and P.A.H. were responsible for the design of the field trial. T.B. was responsible for the trial establishment and data collection. Data analyses were undertaken by J.C.S. and P.A.H., with assistance from B.P. The manuscript was written by J.C.S., B.P. and P.A.H., and reviewed by all co-authors. . We are thus grateful to ARC and FCT for their financial support, which provided the opportunity to complete this study.