Combined Analyses of Phenotype, Genotype and Climate Implicate Local Adaptation as a Driver of Diversity in Eucalyptus microcarpa (Grey Box)

: Trees are a keystone species in many ecosystems and a critical component of ecological restoration. Understanding their capacity to respond to climate change is essential for conserving biodiversity and determining appropriate restoration seed sources. Patterns of local adaptation to climate between populations within a species can inform such conservation decisions and are often investigated from either a quantitative trait or molecular genetic basis. Here, we present ﬁndings from a combined analysis of phenotype (quantitative genetic analysis), genotype (single nucleotide polymorphism (SNP) trait associations), and climate associations. We draw on the strength of this combined approach to investigate pre-existing climate adaptation and its genetic basis in Eucalyptus microcarpa (Grey box), an important tree for ecological restoration in south-eastern Australia. Phenotypic data from a 26-year-old provenance trial demonstrated signiﬁcant genetic variation in growth and leaf traits at both the family and provenance levels. Growth traits were only associated with temperature, whilst leaf traits were associated with temperature, precipitation and aridity. Genotyping of 40 putatively adaptive SNPs from previous genome-wide analyses identiﬁed 9 SNPs associated with these traits. Drawing on previous SNP–climate association results, several associations were identiﬁed between all three comparisons of phenotype, genotype and climate. By combining phenotypic with genomic analyses, these results corroborate genomic ﬁndings and enhance understanding of climate adaptation in E. microcarpa . We discuss the implication of these results for conservation management and restoration under climate change. quantitative traits and do these associations support evidence of climate adaptation?


Introduction
Understanding adaptive variation within species is a major goal for forest tree research, past and present [1]. Given their keystone role in many terrestrial systems worldwide [2], tree responses to changes in climate could have significant implications for wider ecosystems. Despite the fact that trees generally have wide climate tolerances, recent environmental changes have had notable impacts on tree populations, including canopy die-back and mortality [3][4][5][6] and shifts in distribution [7]. Understanding the potential of trees to adapt to climate change will help in the management of forests through approaches such as assisting gene flow and provenancing across climatic gradients [8,9].
At present, genomic approaches are becoming more commonly used for assessing adaptive variation in trees ( [10][11][12]; reviewed in [13]) as well as potential genomic vulnerability to projected  For this study, we sampled seven of the 12 planted provenances across three replicates. Provenances from Victoria and New South Wales (NSW) were sampled, avoiding potentially different species recognised after the trial was established (Queensland provenances potentially E. woollsiana and a NSW provenance from Gillgandra marked as E. pilligaensis (Narrow-leaved grey box), now considered E. woollsiana; [42]) and highly disjunct populations where population structure is likely to confound association analyses (South Australian provenances, disjunct from general south-eastern Australian distribution; Figure 1). Table 1. Sample and climate data for the seven Eucalyptus microcarpa (Grey box) provenances measured in the provenance trial, including number of trees sampled for traits, number of trees genotyped, and climate data of original provenance locations. Climate data from [41]. See Table S1 for climate variable definitions.

Number of Samples (n)
Aridity Growth measurements and leaf samples were taken in October-November 2014. Diameter at breast height (DBH) and number of stems at breast height were measured on all surviving trees from the 71 families available across the seven Victorian and NSW provenances. For multi-stem trees, DBH was calculated based on the total area of all stems at breast height. Leaf samples were collected from a subset of trees for both trait measurement and genotyping. For Benalla, only 10 of the 11 families were sampled for leaf traits to remain consistent with other provenances. Approximately two randomly-selected trees per family per replicate were sampled, giving approximately six trees per family in total or 459 trees across 70 families from seven provenances sampled (Table 1). For trait measurements, 10 fully expanded, mature leaves were sampled from a single mid-point in the outer canopy. A single canopy point was used owing to the challenge of sampling from tree canopies >10 m high. To ensure this did not impact results, a second set of 10 leaves was collected from a different side of the canopy for 32 trees across families and replicates. Assessment of trait variation between canopy points for leaf area, length, weight and thickness found little difference in traits between canopy points and greater variance between individual trees than between canopy points (see Supplementary Methods, Tables S2 and S3 for more detail). This suggested most trait variation occurred between rather than within trees and thus the single canopy point was sufficient for the questions of this study. All leaves for trait measurements were stored at 4 • C during fieldwork before later being dried at 40 • C. Several additional leaves were collected and immediately dried on silica gel for later DNA extraction and genotyping.
Dried leaves for trait measurements were weighed and leaf area and length measured using a Licor (USA) LI-3000C leaf area meter and leaf thickness measured using a micrometer (Mitutoyo, Japan). Leaf weight, area, length and thickness measurements were averaged across the 10 leaves to create a single measurement per tree. Specific leaf area (SLA; area/weight) and density ([weight/area]/thickness) were calculated from average dry leaf area, weight and thickness measurements. Height was measured using a clinometer (Suunto, Finland) for the subset of trees with leaf samples and a set of additional trees (total 496; Table 1)

Genetic Variance
To estimate the proportion of phenotypic variance due to genetic effects, family and provenance were fitted as random terms in a bivariate linear mixed model framework for each trait as per [43]: where y is the vector of trait observations, b and u are vectors of fixed and random effect estimates, X and Z are incidence matrices for fixed and random model terms and e is a vector of random residual effects. Models were tested with fixed effects of Replicate and 'Edge'-a binary variable defining a tree's location on the outer row or column of a replicate block. 'Edge' accounted for potential edge effects due to inconsistent boundary tree planting and survival around the trial plots. The number of stems at breast height was included as a fixed effect for growth traits of DBH, height and size ratio. Provenance, as a random effect, was found to significantly improve models for all traits, based on likelihood-ratio tests, and was therefore retained in all models. A Wald test was used to test for significance of fixed effects, with only those fixed effects identified as significant (p < 0.05) retained for that trait model. Final models per trait (Table S4) were checked for linearity, homoscedasticity and normality of residuals by plotting fitted vs. residual values and a histogram of residuals. Narrow-sense heritability (h 2 ) was calculated as: where V A is the additive genetic variance, V P the total phenotypic variance, σ 2 f the family-level variance and σ 2 residual the residual variance of the model. A relatedness coefficient of 1/2.5 was used to account for greater relatedness than the expected 1/4 for halfsibs commonly found within open-pollinated, half-sib families of eucalypts [44,45]. Within family relatedness estimates from SNP data in this study (see Section 2.3.1 Genotyping) suggest this is also appropriate for E. microcarpa (average within family relatedness coefficient = 0.32 ± standard deviation (s.d.) 0.09). This estimate is an approximation of narrow-sense heritability as it includes potential maternal and dominance effects. Initial model testing for heritability estimates was performed in R v 3.2.1 [46], using the R package 'lme4 v 1.1-12 [47], before final calculations in ASReml-R (Release 3; [48]) using the pin.R function to estimate standard errors (SE; [49]).
To estimate genetic differentiation of quantitative traits between provenances, Q ST was calculated from estimated variance components as per [50]: where σ 2 p is the genetic variance between populations and V A and σ 2 f are as defined above. Standard errors for Q ST were calculated in ASReml-R using the pin.R function. Q ST was not estimated for leaf density due to a lack of additive genetic variance for this trait.
To determine whether provenance-level differentiation was potentially due to selection or neutral processes, trait Q ST values were compared to a distribution of neutral F ST values. The expected distribution of Q ST for neutral traits may be approximated by the distribution of F ST for neutral genomic loci [50]. If Q ST is significantly greater than the mean neutral F ST , a trait may be considered significantly more differentiated than neutral expectations and potentially under diversifying selection. The distribution of neutral F ST can be well approximated by a χ 2 distribution; specifically whereF ST = per locus F ST , F ST = mean F ST and n = number of demes [51]. Individual Q ST estimates (as percentages i.e., Q ST × 100) were deemed significantly greater than neutral expectation if the lower SE was significant (p < 0.05) when tested against χ 2 with df = F ST (as a percentage). Mean F ST was calculated from 418 'neutral' SNPs in 157 samples from the seven original provenance populations [12]; see Supplementary Methods). To avoid limitations associated with estimating a distribution of F ST values from a single mean value [50], Q ST estimates were also compared to the empirical F ST distribution of the 418 'neutral' SNPs and similar results were obtained (see Supplementary Methods, Figure S1).

Genetic Trait-Trait Correlations
To explore evidence of a shared genetic basis between traits, trait-trait correlations were performed in ASReml-R to test genetic correlations between traits at both the family and provenance levels, using the mixed model framework previously described. Leaf density was not tested as no additive genetic variance was found. Models included replicate, edge and stem as fixed effects, and covariance between traits was estimated by fitting random family and provenance terms and residuals across traits in an unstructured covariance matrix, assuming heterogeneous variance estimates. Genetic correlations (r) were calculated as: where cov 1,2 is the covariance between trait 1 and trait 2 and σ 2 1 and σ 2 2 are the variance of traits 1 and 2 respectively. A two-tailed likelihood-ratio test, testing against the null model of covariance = 0 ('diag' covariance matrix for random effects), was used to calculate the significance of genetic correlations. To account for multiple testing, q-values were calculated using the R package 'qvalue' v 1.43.0 [52] within each test level separately.

Climate Associations
To investigate climate as a potential selection pressure driving provenance-level trait genetic differentiation, climate associations between quantitative traits and climate of the original provenance location (Table 1) were tested using linear models in R. Ten climate variables previously used in a genomic assessment of climate adaptation in E. microcarpa were tested [12]; (Table S1). Briefly, climate variables were selected based on projected future hotter and drier climatic conditions in south-eastern Australia, increased extreme events, and variable precipitation patterns [53]. Climate variables related to temperature, precipitation, evaporation, moisture and aridity were reduced to a set of 10 key representative variables, minimising redundancy and correlations between variables (Table 1). Climate data for the original provenance locations were downloaded from the Atlas of Living Australia [41] (0.01 degree (~1km) gridded data). Provenance-level Best Linear Unbiased Estimates (BLUEs, equivalent to least-square means) were calculated for each trait using the R package 'lsmeans' v 2.2.3 [54] from trait specific models (Table S4) with Provenance fitted as a fixed effect.
To reduce correlations between traits, BLUEs for growth traits (DBH, height and size ratio) and leaf traits (area, length, thickness, weight, density and SLA) were reduced via separate principal component analyses (PCA) using the R package 'FactoMineR' [55]; see Table S5 for trait correlations and contributions to principal component axes). The first two axes (principal components, PC) for growth (98.4% of variation) and first three axes for leaf trait (97.9% of variation) were tested for associations with climate variables using separate linear models ('lm' function) in R. An analysis of variance (ANOVA) was performed in R ('anova' function) to calculate the p-value for the environmental component of the model. Correction for multiple testing was applied to environmental component p-values using 'qvalue' in R. To ensure associations were not due to neutral population structure, linear models were repeated including latitude and longitude of the original provenance location. This assumes a population structure pattern of isolation-by-distance, which appears to be the case for E. microcarpa [12,56].

Genotyping
To look for links between potentially adaptive genomic variation and quantitative traits, a set of 40 putatively adaptive candidate SNPs were genotyped in 422 trees for which both growth and leaf trait data were collected (Table 1). Using data from genome scans, rather than candidate genes based on known gene function, has the potential to uncover new genomic regions associated with traits, without the need for a priori selection of potential genes of interest [57]. Candidate adaptive SNPs were chosen based on previous genomic analyses of climate adaptation in E. microcarpa, using reduced-representation, DArTseq SNP data [12]. Potential candidate SNPs used in this study were those identified as an F ST outlier in at least one of the four previous analyses-BayeScan [58], hierarchical FDIST2 [59], FDIST2 [60,61] or Bayenv2 X T X [62], or those identified as having a strong association with at least one of the 10 climate variables tested (Bayenv2; [62,63]). In the previous analyses, environmental association analyses were only performed on SNPs identified as F ST outliers. To expand the list of potential candidate SNPs for genotyping in this study, environmental associations were performed in Bayenv2 for all 4218 SNPs, using the same methodology as the previous analysis [12]. The F ST outlier analyses account for population structure using various demographic models to create a 'null' distribution against which individual SNPs are compared, whilst Bayenv2 employs a covariance matrix, generated from genomic data, to account for population structure when performing SNP-climate associations. Candidate SNPs were, therefore, those that were more differentiated than expected (F ST outliers) or were associated with climate even after accounting for population structure (Bayenv2).
To determine a subset of 40 putatively adaptive candidate SNPs for genotyping, potential candidate SNPs were filtered to those with (i) a minor allele frequency of greater than 0.05 in at least six of the seven original provenance sites to ensure a high enough frequency across the dataset to enable effective SNP-trait association analyses, (ii) to SNPs within 2000 basepairs (bp) of a putative E. grandis gene with an Arabidopsis thaliana orthologue (TAIR10) based on E. grandis v 1.1 annotation [64] to enable exploration of potential gene function or physiological processes underlying traits, and (iii) those with an appropriate flanking sequence for DArTmp, the genotyping methodology used (see below).
An additional 25 'non-adaptive' SNPs were included as a proxy for neutral population structure. Based on previous analysis, 'non-adaptive' SNPs were selected from those not significant in all four F ST outlier tests-BayeScan (q > 0.2), hierarchical FDIST2 (p > 0.1), FDIST2 (q > 0.2) and X T X (outside top 10% of 3 runs), not strongly associated with any of the 10 climate variables tested (Bayenv2, Bayes Factor (BF) < 20), and with a minor allele frequency of greater than 0.05 in at least six of seven original provenance sites. Whilst SNPs deemed non-significant in F ST outlier tests may still be under selection [65], this is likely to be a weaker selection and, therefore, not strongly impact the primary aim of these SNPs, being to approximate neutral structure. Despite the significantly smaller dataset, the final dataset of 65 SNPs provided a fair representation of the population structure found in the full dataset of 4218 SNPs [12]-correlation of r = 0.69 for provenance-level pairwise F ST estimates using the 65 genotyped SNPs within the provenance trial compared to the full dataset of 4218 SNPs for the seven original natural sites (pairwise F ST calculated in Arlequin v 3.5.1.2 [66]).
Flanking sequences for each SNP were extracted from genomic read data from the previous analysis [12]. Firstly, a E. microcarpa consensus sequence was created from genomic read data flanking the SNPs of interest by modifying the E. grandis v 1.1 sequence using samtools mpileup v 1.2-10 and bcftools v 1.3 [67]. E. microcarpa sequences for each SNP locus were then created by extracting flanking regions covered by known E. microcarpa reads from the "whole genome" consensus.
To genotype the SNPs, approximately 25 mm 2 of dried leaf material per tree and fasta sequences of SNP loci were sent to Diversity Arrays Technology Pty Ltd. for DNA extraction and subsequent SNP genotyping via DArTmp multiplex PCR. Genotypes were called from allele sequencing read counts using the maximum likelihood method described by [68] in R, using an error rate of 0.05 and only scoring genotypes with a minimum total read depth of 20. Data were filtered to SNPs and individuals with <50% missing data. Little linkage was found between SNPs (average r 2 = 0.003 ± s.d. 0.005), nor within adaptive or neutral SNP subsets (r 2 = 0.004 ± 0.006 and r 2 = 0.004 ± 0.005, respectively; linkage (as measured by correlation, r 2 ) was calculated using the 'r' function in PLINK v 1.90b3p [69]).

Genotype-Phenotype Associations
SNP-trait associations between the 40 putatively adaptive SNPs and nine scored traits were performed on 422 individuals in TASSEL v.5 [70]. To account for population structure, a PCA based on covariance was performed in TASSEL using all 65 SNPs. Overall F ST , estimated via analysis of molecular variance (AMOVA) in Arlequin using all 65 SNPs, found 3% of variance between populations. Therefore, the first PC accounting for the highest amount of variance and reflecting between-population variance (PC1 = 4.1%) was retained in all analyses.
To account for relatedness between individuals, a kinship matrix was created in Coancestry v. 1.0.0.1 [71]. Coefficient of relatedness (2θ) was calculated using the triadic likelihood estimator [72] and using 100 reference individuals. This estimator can account for potential inbreeding, an important consideration given the mixed mating system of eucalypts [73]. Inbreeding (f ) was calculated using the estimator of [74], which bounds inbreeding coefficient estimates between 0 and 1. Two different kinship matrices were tested. The first used individual kinship estimates, allowing for variation in the degree of relatedness between all individuals, including across families and provenances. However, due to the small number of SNPs used in this study, relatedness estimates may be subject to large error rates. In addition, identity-by-state (genetic similarity) may be mistaken as identity-by-descent (true relatedness), especially between provenances where true sibling or cousin relationships are unlikely. To address these issues, an alternative matrix was created using the average family-level relatedness within and between families. This can reduce the effects of individual spurious estimates, especially between provenances; however, it may reduce genuine unique relationships. The average relatedness (2θ) within and between each family and the average inbreeding (f ) within each family was calculated and an individual matrix created using family-level averages.
SNP-trait associations were tested in TASSEL, using a mixed linear model (MLM) based on the same underlying framework as previously described for variance component estimation (Equation (1)); with trait observations as the dependent variable (y), provenance, replicate, edge and stems as well as adaptive SNP and one PC to account for population structure as fixed effects, and kinship included as a random effect to account for relatedness between individuals. Both the raw individual kinship matrix and family-average kinship matrix were tested.

Genotype-Phenotype-Climate Analysis
To provide additional support for climate adaptation in E. microcarpa, and possible underlying genetic mechanisms, we compared results from the three independent pairwise associations of genotype, phenotype and climate performed in this and previous analyses [12]. Results of SNP-climate associations from the previous study were compared to trait-climate and SNP-trait associations in this study. For trait-climate comparisons, climate association results for the principal component(s) to which individual traits loaded were used (Table S5).

Genetic Variation Within Quantitative Traits
Significant additive genetic variance suggested a genetic basis for quantitative trait variation observed in E. microcarpa. Family-level heritability (h 2 ) ranged from 0.106 to 0.373 (Table 2). Significant heritability was suggested for all traits, except leaf area, with estimates more than one standard error (SE) from 0. Additive genetic variance was not found for leaf density. Specific leaf area (SLA) and diameter at breast height (DBH) had the highest heritabilities (0.373 and 0.318, respectively), followed by leaf weight and tree height (0.230 and 0.210, respectively). Heritability estimates for other traits were below 0.2 ( Table 2). Errors for heritability estimates tended to be large, potentially reflecting low replication per family.
Significant trait differentiation was found between provenances with Q ST greater than one SE from 0 for all traits, excluding leaf density (Table 2). Furthermore, Q ST estimates for tree height, leaf area and leaf length exceeded neutral F ST (F ST 'neutral' SNPs = 0.033), suggesting these traits are potentially under divergent selection (Table 2; Figure S1). Low, neutral F ST here is in line with low differentiation identified in E. microcarpa using greater numbers of loci and samples [12]. Q ST ranged from 0.094 to 0.483, with leaf area and leaf length having the highest Q ST estimates (> 0.4), followed by leaf weight and thickness (0.224 and 0.205, respectively). For all other traits, Q ST was less than 0.2. In addition to significant variation between provenances, Q ST estimates also suggested greater than 50% of trait variation could be attributed to variation within provenances (Table 2, Figure S2).

Genetic Trait-Trait Correlations
Not all traits were genetically independent. Correlations in additive genetic variance suggested variation in some pairs of traits had a common heritable basis (Table 3). At the family level, all three growth traits were significantly correlated (|r| = 0.75 to 0.97, q < 0.1), as were all six leaf traits (|r| = 0.36 to 0.98, q < 0.1; Table 3). Between growth and leaf traits, DBH and leaf traits were independent except for leaf area (r = −0.20), height was correlated with all leaf traits except leaf weight (r = −0.30 to −0.56), and size ratio was correlated with leaf weight and thickness (r = −0.55 and −0.42 respectively; Table 3). In contrast, most traits were genetically independent at the provenance level, with only five pairs of traits having significant genetic correlations (r = 0.79 to -0.95, q < 0.1; Table 3). Genetic correlations between DBH and height (r = 0.93), leaf area and leaf length (r = 0.79), and leaf area and leaf weight (r = 0.95) reflected family-level results. Leaf length and SLA were also significantly correlated at the provenance level; however, contrary to family-level results, greater leaf length was correlated with higher SLA (r = 0.93). Although not significant at the family-level, higher leaf weight was correlated with larger DBH at the provenance level (r = 0.92).

Climate Associations with Quantitative Trait Variation
Significant heritable genetic variation (h 2 ) and genetic differentiation between provenances (Q ST ) suggest genetic variation in phenotypes that selection can act on, and that selection has likely driven trait variation in E. microcarpa. An adaptive hypothesis was further supported by associations between phenotypic variation and climate among provenances, implicating local adaptation to climate as a driver of genetic variation in phenotypes.
Provenance level trait estimates (BLUEs) were reduced to two major principal component axes for growth traits and three major axes for leaf traits (Table S5). Nine of 50 potential associations between growth or leaf trait axes and climate were significant, although only marginally after correcting for multiple testing (p < 0.05, q < 0.15; Table 4a; Table S6 for individual trait-climate associations). Stronger associations were found between climate and leaf trait axes than growth trait axes (Table 4a; Figure 2). For growth traits, the second major axis, predominately driven by increasing size ratio and to a lesser extent decreasing DBH (Table S5), was positively correlated with warmest period maximum temperature and mean annual temperature. For leaf traits, the second major axis, driven primarily by decreasing leaf density and increasing SLA and to a lesser extent increasing leaf length and thickness, was positively correlated with the three temperature variables and was negatively correlated with both aridity variables and winter precipitation. The third leaf trait axis, primarily reflecting increasing leaf weight and thickness, was negatively correlated with summer precipitation. driven trait variation in E. microcarpa. An adaptive hypothesis was further supported by associations between phenotypic variation and climate among provenances, implicating local adaptation to climate as a driver of genetic variation in phenotypes.
Provenance level trait estimates (BLUEs) were reduced to two major principal component axes for growth traits and three major axes for leaf traits (Table S5). Nine of 50 potential associations between growth or leaf trait axes and climate were significant, although only marginally after correcting for multiple testing (p < 0.05, q < 0.15; Table 4a; Table S6 for individual trait-climate associations). Stronger associations were found between climate and leaf trait axes than growth trait axes (Table 4a; Figure 2). For growth traits, the second major axis, predominately driven by increasing size ratio and to a lesser extent decreasing DBH (Table S5), was positively correlated with warmest period maximum temperature and mean annual temperature. For leaf traits, the second major axis, driven primarily by decreasing leaf density and increasing SLA and to a lesser extent increasing leaf length and thickness, was positively correlated with the three temperature variables and was negatively correlated with both aridity variables and winter precipitation. The third leaf trait axis, primarily reflecting increasing leaf weight and thickness, was negatively correlated with summer precipitation.  Table 4). Refer to Table S5 for further details on trait contributions and correlations to PC axes.  Table 4). Refer to Table S5 for further details on trait contributions and correlations to PC axes. After accounting for possible neutral demographic effects, associations with growth traits were no longer significant (Table 4b). However, local adaptation cannot be ruled out due to strong spatial autocorrelation of the two temperature variables and latitude (r = 0.91 and r = 0.90 for mean annual temperature and warmest period maximum temperature respectively). Associations with leaf traits remained significant, although p-values were higher. Additional associations were identified when accounting for population structure-the second leaf trait axis was also associated with summer precipitation, precipitation of the driest period and annual precipitation, although these were not significant after correcting for multiple testing (p < 0.05, q > 0.15; Table 4b).

Genotype-Phenotype Associations
Associations between putatively adaptive SNPs and quantitative trait variation further support a genetic basis to trait variation, validating results of previous genomic adaptation analyses [12] and providing insights into possible genomic regulatory factors. Nine of the 40 putatively adaptive SNPs genotyped were significantly associated with at least one of the nine measured traits (Table 5). Arabidopsis thaliana orthologues (TAIR10) of putative E. grandis genes within 2000 bp of associated SNPs (Table 5; [75]) suggest trait variation in E. microcarpa may be influenced by genes associated with development (e.g., CHC1, [75]) and stress responses (e.g., PRK5, [76]) as well as gene regulation (e.g., MYB78, [77]). For growth traits, three significant SNP-trait associations (p < 0.05) were identified, each explaining approximately 0.9%-1.9% of trait variation. For leaf traits, 11 significant associations were identified ( Table 5). The amount of leaf trait variance explained by an individual SNP varied from 0.9%-3%. No SNP-trait associations were significant after correction for multiple testing across SNP-trait comparisons.

Genotype-Phenotype-Climate Associations
Support for climate adaptation was strengthened by four climate-related associations that were corroborated across the three independent analyses of trait-climate, SNP-trait and, drawing on previous analyses, SNP-climate associations [12] (examples in Figure 3). In these cases, the results lend weight to a hypothesis for climate adaptation of genetic variation underlying traits in E. microcarpa and suggest potential candidate genes that may have been targeted by selection. It should be noted candidate genes are interpreted here as indicative only and warranting further analysis, given low linkage disequilibrium in eucalypts [78] and assumptions of sufficient similarity between E. microcarpa and E. grandis genomes and A. thaliana orthologues reflecting possible gene functions. previous analyses, SNP-climate associations [12] (examples in Figure 3). In these cases, the results lend weight to a hypothesis for climate adaptation of genetic variation underlying traits in E. microcarpa and suggest potential candidate genes that may have been targeted by selection. It should be noted candidate genes are interpreted here as indicative only and warranting further analysis, given low linkage disequilibrium in eucalypts [78] and assumptions of sufficient similarity between E. microcarpa and E. grandis genomes and A. thaliana orthologues reflecting possible gene functions.  Two SNPs associated with leaf length also had several climate associations in common with this trait (SNP 2:63702271, 5:5009176; Table 5; Figure 3). Leaf length correlated primarily with the first, but also the second leaf trait axis. Both leaf length, via its correlation to the second leaf trait axis and both SNPs were associated with mean annual temperature and SNP 2:63702271 was also associated with warmest period maximum temperature. In addition, the second leaf trait axis and SNP 5:5009176 were associated with aridity. Based on A. thaliana orthologues of predicted E. grandis genes, possible genes or functions that may influence climate-related variation in leaf length include an MYB domain protein (MYB78; SNP 2:63702271), and proteins of unknown function or hypothetical proteins (SNP 5:5009176; Table 5).
Results for leaf density suggest a possible link to SNP 6:39617441 (putative A. thaliana orthologue, C2 calcium/lipid-binding plant phosphoribosyltransferase family protein; Table 5) and temperature, with both the SNP and leaf density, via its correlation with the second leaf trait axis, associated with mean annual temperature and warmest period maximum temperature.
Finally, a potential link was found between size ratio and SNP 10:29282238, supported by both the trait, via correlation with the second growth trait axis, and the SNP associated with mean annual temperature and warmest period maximum temperature. Based on A. thaliana orthologues of predicted E. grandis genes, a possible SWIB/MDM2 domain or galactose oxidase/kelch repeat superfamily protein may be involved in temperature-related variation in size ratio (Table 5).

Discussion
Building on previous genomic-only assessments of adaptive variation [12], we found evidence for climate driven, local adaptation of genetic variation in growth and leaf traits of E. microcarpa. The identification of adaptive phenotypic variation provides support for local adaptation to climate within E. microcarpa and demonstrates that the genomic variation previously identified likely reflects adaptive phenotypic variation. Furthermore, the identification of both genetic variance and trait variation within provenances suggest there is variation that selection could act upon, potentially enabling adaptation in the future. Although representing a small sample of the genome as well as trait and population diversity, the combination of both phenotypic and genomic analyses bolstered evidence of climate-related, adaptive variation in this species that warrants further investigation, especially for understanding future adaptive potential and seed sourcing under climate change.

Evidence of Genetic Variance and Climate Adaptation
In line with many ecological traits in eucalypts [19], this study detected a genetic component to variation in several traits of E. microcarpa, upon which selection could potentially act. Heritability estimates for growth traits in this study were similar to, or higher than, estimates seen in other eucalypts (e.g., DBH h 2 = 0.06-0.37, height h 2 = 0.13; [22,35]), although estimates for leaf traits were low relative to other studies (e.g., leaf length h 2 = 0.60-0.77; [36,79]). Furthermore, genetic trait-trait correlations identified in this study suggest a common genetic basis, or pleiotropy, between leaf and growth traits. Genetic variance and correlation estimates in this study are likely influenced by low replication within families, increasing the standard error, as well as potential upward bias in additive variance due to maternal, dominance and epistatic effects unaccounted for [80]. Further analysis using more samples per family, as well as more test families and provenances, would greatly improve genetic variance estimates. Despite these limitations, genetic variance detected for leaf and growth traits in E. microcarpa here suggests potential for future evolution in response to environmental change.
Greater quantitative genetic trait differentiation between provenances than expected under a neutral model suggest selection has driven adaptive differences in a number of traits in E. microcarpa. This result is consistent with evidence of selection on leaf and growth traits commonly found in eucalypts as well as other tree species [19,24,35,36,81]. Whilst significant adaptive divergence in quantitative traits was detected between provenances, a large proportion (greater than 50%) of total genetic variance was attributed to variation within provenances. The presence of such standing variation further suggests a genetic source for future adaptation in these traits.
Associations between trait variation and climate suggest local adaptation in E. microcarpa resulting from divergent selection. This supports the results of genomic adaptation analyses in E. microcarpa [12] suggesting significant, climate-associated genomic differences and evidence of phenotypic variation driven by similar climatic variables in other widespread trees species [25,28,35,36]. Trait associations with maximum temperature identified here and in other studies [20,36], suggest climate extremes may be important drivers of adaptive trait differences in trees. In line with a priori expectations of thicker, denser leaves under higher water stress [82,83], thicker leaves in this study (second leaf trait axis) were associated with warmer, more arid home climates (consistent with [24]). However, contrary to this expectation, higher leaf density and lower SLA (second leaf trait axis) were associated with a cooler, wetter climate of origin. Plastic responses of leaf density and SLA have been found to differ between provenances and species of eucalypts [24,39,82]. Leaf density and SLA climate associations in this study may therefore reflect variation in plastic responses to environment at the trial site, especially in the case of leaf density for which no additive genetic variance was found.
One of the challenges for association studies in trees is disentangling local adaptation from neutral processes, especially where correlations between climate and population structure exist. False positives-differences due to drift rather than adaptation-can arise if population structure is not accounted for. Conversely, correcting for population structure that correlates with environment reduces the environmental effect and can result in false negatives, as shown in simulation studies of genotype-environment association analyses [65]. In this study, higher p-values when latitude and longitude were included in trait-climate associations likely reflect reduced power rather than a lack of association. Even with reduced climate effects due to correlations with spatial variables and reduced power (fewer degrees of freedom) given the small sample size (seven provenances), several associations remained marginally significant (p = 0.034-0.073; Table 4). Furthermore, a lack of significance after correcting for multiple testing potentially reflects a further reduction in power for this small study rather than a true lack of association. Indeed, differentiation (Q ST ) greater than neutral population structure supports local adaptation rather than neutral processes. Future analyses, with a greater number of provenances will help to more clearly define those climate variables driving local adaptation in E. microcarpa.

Linking Genotype and Phenotype
Whilst genomics is providing advances for identifying signatures of adaptation, assessing phenotypic variation in combination with genomic data provides greater support for adaptive variation and increases understanding of potential adaptive mechanisms [16,26,84]. SNP-trait associations in this study linked previous genomic analyses [12] with quantitative analyses of climate adaptation in E. microcarpa, validating earlier results, providing stronger evidence of climate as a driver of local adaption and highlighting potential genes or genic regions underlying quantitative traits.
Combining genomic and phenotypic analyses firstly demonstrated that putatively adaptive genomic markers explained, individually, a small proportion of genetic variation in quantitative traits. This, as well as multiple SNPs associating with individual traits, is in line with the highly polygenic nature of quantitative traits, and characteristic of association studies in trees [18,19,85]. Indeed, variance estimates found here are comparable to other SNP-trait association studies in trees [25,29,81]. The small number of SNP-trait associations identified and the weak significance found is likely due to the small size and power of this study, although they may also reflect the tendency of single marker analyses to miss associations due to small effect sizes [84,86]. Similarly, non-significance after correcting for multiple testing may be due to corrections being too stringent given the small study size, combined with the more subtle associations underlying highly polygenic traits. Future analyses are warranted before concrete conclusions are drawn from these results, including studies which employ more genomic markers (e.g., [18,84]), that use polygenic [85] or "omnigenic" [87] assessments of complex traits instead of single SNP associations (e.g., [88]) and that assess a wider range of traits, including phenology and physiological traits. Furthermore, additional common garden sites, in contrasting environments, would assist in separating plastic from genetic responses as well as identifying gene-by-environment interactions that were not able to be separated in this study. Lastly, multi-species analysis could provide additional evidence for genomic regions and mechanisms underlying climate adaptation (e.g., [18]).
As found in other studies [28,29], our results demonstrated how combining both genomic and phenotypic datasets can validate genomic signatures of adaptation and reveal adaptive variation and potential within a species. Phenotypic variation in this study complement and validate previous genomic variation, together demonstrating genetic variation between provenances (heritability in this study and [12]) and high variation within provenances (Q ST and [56]) and thus suggesting the presence of variation upon which selection may act. Results of this study, therefore, support the combination of traditional common garden trials with modern genomic technologies as a method for investigating local adaptation and the potential for future adaptation in trees [16,26].

Conservation and Restoration under Climate Change
Incorporating evolutionary potential into conservation management will be essential for long-term sustainability in natural systems under environmental change [89]. In trees, high standing genetic variation, fecundity, gene flow, and large effective population size are likely to help enhance evolutionary potential by maintaining variation upon which selection can act [81,90]. Supporting previous findings [12,56], significant additive genetic variance and trait variation found in this study suggests high levels of standing genetic variation and therefore potential for future adaptation in E. microcarpa. Additionally, adaptive differences between provenances provide an opportunity for the movement of 'pre-adapted alleles', via gene flow or assisted migration, to help facilitate climate adaptation [90]. However, whilst additive variance is important for adaptation, it does not necessarily equate to adaptive capacity [91]. If genetic correlations identified in this study represent true genetic dependence between traits, pleiotropy may also constrain adaptive responses in one trait in favour of increased fitness in another [19]. Furthermore, additional stresses such as habitat fragmentation, reducing population size and potentially disrupting gene flow, could slow the rate of adaptation within tree populations relative to climate change [92]. Where adaptation may not occur quickly enough, introduction of pre-adapted genetic variation through assisted migration and restoration planting may be necessary to assist evolutionary potential in natural populations [8].
As restoration seed sourcing strategies move away from 'local' and toward capturing evolutionary potential, especially pre-adapted climate-related diversity [9], understanding climate-related variation across species' distributions will help inform effective seed sourcing for long-term sustainability of populations under climate change. The results of this study, together with genomic results [12] suggest the presence of climate-associated genetic variation within E. microcarpa that could be utilised to enhance diversity and pre-adapted climate variation within restoration plantings and the wider landscape. In particular, temperature, aridity and summer and winter precipitation appear to be important climate drivers of adaptation in E. microcarpa. Sourcing seed along these climatic gradients, towards projected future climates, may therefore enhance the long-term potential of restoration sites as well as support adaptation in the wider landscape through the introduction of pre-adapted genetic variation.

Conclusions
Assessing phenotypic variation provides additional support and confirmation of previous genomic-based assessments of adaptive variation as well as climate as a driver of this variation. The presence of adaptive variation within and between populations suggests potential future adaptation in natural populations and may enable selective seed sourcing to enhance climate resilience in restoration plantings. However, whether the variation found is sufficient for adaptation to future climate change remains unknown. Further quantitative trait and genomic studies, especially whole-genome sequencing and polygenic assessments of adaptation, are needed to clarify the extent of adaptive variation in this species and thus potential adaptability, or vulnerability, to climate change.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/5/495/s1, Supplementary Methods, Table S1: Definitions for climate variables used in climate association analysis of quantitative traits and putatively adaptive SNPs in Eucalyptus microcarpa (Grey box).; Table S2: Results of two sample t-tests assessing differences in leaf trait measurements between two different canopy points within a single tree (10 leaves per canopy point) for Eucalyptus microcarpa (Grey box).; Table S3: Mixed model estimates of variance between trees (Tree) and between two canopy points within a tree (Position) for leaf trait measurements in 32 Eucalyptus microcarpa (Grey box) trees, including fold difference between variance estimates.; Table S4: Mixed models for calculating narrow-sense heritability (h 2 ) and Q ST in nine quantitative traits measured in Eucalyptus microcarpa (Grey box).; Table S5: Trait correlations (a) and contributions (b) to trait principal component (PC) axes in Eucalyptus microcarpa (Grey box).; Table S6: Climate associations for individual leaf and growth traits in Eucalyptus microcarpa (Grey box), calculated (a) without and (b) with adjustment for spatial population structure; Figure S1: Empirical per locus 'neutral' F ST distribution (bar plot) for Eucalyptus microcarpa (Grey box) and Q ST estimates; Figure S2: Provenance level, Best Linear Unbiased Estimates (BLUEs) for nine quantitative traits measured in Eucalyptus microcarpa (Grey box).
Data Availability: Raw trait data for all measured trees (FileS2_Emicrocarpa_trait_data.xlsx) and raw genotypes for all genotyped individual SNP data for the 65 SNPs used in this study (FileS3_Emicrocarpa_SNPgenotypes.vcf; FileS4_README_vcf_file.txt) are available from the CSIRO Data Access Portal: https://doi.org/10.25919/ 5e9e4b5b2241e.