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

Abstract

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 findings 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 significant 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 identified 9 SNPs associated with these traits. Drawing on previous SNP–climate association results, several associations were identified between all three comparisons of phenotype, genotype and climate. By combining phenotypic with genomic analyses, these results corroborate genomic findings and enhance understanding of climate adaptation in E. microcarpa. We discuss the implication of these results for conservation management and restoration under climate change.

1. 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 change [14,15]. Genomics provides the advantage of enabling investigations of adaptation in shorter time frames and without the logistical investments required for common garden experiments [1,16]. Furthermore, genomics can identify genes and genomic regions that have responded to selection [17,18].

Nevertheless, it remains important to understand phenotypic variation to validate genomic results and provide insight into mechanistic drivers of adaptive variation. Common garden trials have been extensively used to understand adaptive variation in trees [19], demonstrating significant genetic trait variation along environmental clines in numerous traits, for example growth [20], leaf morphology [21], stress tolerance [22], physiology [23], and even plasticity [24]. Whilst trait variation can be influenced by multiple selection pressures, climate (and especially temperature, precipitation and aridity) has been identified as a key driver of adaptive variation [23,24,25].

Evidence for adaptation can be built by combining both genomic and phenotypic approaches, corroborating evidence of adaptation across all three components of genotype, phenotype and climate [16,26,27]. Such an approach has been useful for several tree species: for example, relationships have been established through a combined approach between outlier single nucleotide polymorphisms (SNPs), temperature and leaf traits in Alnus glutinosa (Black alder) [28], between aridity-associated adaptive genetic variation and leaf traits in Eucalyptus tricarpa (Red ironbark) [11], and between candidate genes, temperature and bud cold tolerance in Pseudotsuga menziesii (Douglas fir) [29]. Demonstrating phenotypic variation and linking climate-adapted traits to climate-adaptive genotypes enables validation of genomic signatures of selection and provides insights into the potential role and function of genomics in mediating phenotypes under selection [1,26,27].

Eucalypts are the dominant tree in the majority of forests and woodlands in Australia [30]. Furthermore, they are widely used in ecological restoration plantings [31]. For many eucalypt species, climate change is projected to decrease the extent of climate-suitable habitat [32]. Migration is unlikely to match rates of climate change, hindered by the long generation times of trees, altered land use and low topographic relief in Australia, the latter making the physical distances required to track current climate niches greater than in areas of high topographic relief e.g., mountains with strong elevation gradients [33,34]. Understanding current climate adaptation will be important for utilising adaptive variation to enhance evolutionary potential in both natural and restored sites for these keystone species.

Adaptive genetic variation has been identified in eucalypts through both phenotypic and genomic approaches. There is an extensive history of quantitative trait studies in eucalypts which have demonstrated genetic variation in many traits related to both abiotic and biotic factors [19], including growth and morphology [35,36], phenology [37], physiology [38], and plasticity [24,39]. Previous quantitative trait loci (QTL) mapping has identified genomic regions associated with particular traits [19]. More recently genomic analyses have revealed potential genes and genomic regions associated with climate adaptation [11,12,40]. Whilst phenotypic and genomic approaches have been used to investigate climate adaptation in eucalypts, few studies have combined them (however, see [11]).

In this study, we combined phenotypic and genomic data to infer evidence of climate adaptation in Eucalyptus microcarpa (Maiden) Maiden (Grey Box). Eucalyptus microcarpa is an important restoration species in south-eastern Australia, used to mitigate widespread habitat loss across the species’ distribution due to clearing for agriculture. Previous landscape genomic analyses identified genomic signatures of climate adaptation in E. microcarpa [12] including genomic variants associated with temperature, precipitation and aridity. Building on these results, this study aimed to investigate evidence of variation in quantitative traits to provide a greater understanding of adaptive variation in this species. It then aimed to further support genomic results by combining genotypic, phenotypic and climate association analyses. In particular, this study asked the following questions: (1) Is there evidence of quantitative genetic variation in growth and leaf traits of E. microcarpa? (2) Is variation in these quantitative traits associated with climate? (3) Are putatively adaptive genomic variants associated with quantitative traits and do these associations support evidence of climate adaptation?

2. Materials and Methods

2.1. Trial Design and Trait Measurements

A provenance trial for Eucalyptus microcarpa was established by the Western Australian Department of Conservation and Land Management in 1988 near Collie, Western Australia (Figure 1). Though outside the species normal distribution and established for a different purpose, this trial provided an opportunity to explore trait variation in well-established, mature trees in a common environment. Temperature at the trial site was mid-range compared to original provenance locations (

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 Index (Ratio) 1		Precipitation (mm)					Temperature (°C)
Site	State	DBH	Height	Leaf Traits	Genotyped	Mean Annual	Maximum Month	Annual (Bio12)	Summer	Winter	Driest Period (Bio14)	Wettest Period (Bio13)	Annual Mean (Bio01)	Max. Month Abs. Mean Max.	Warmest Period Max. (Bio05)
Avoca	Vic	251	65	62	58	0.740	1.797	529	103	176	6	14	13.7	43	28.8
Benalla	Vic	289	69	67	58	0.698	1.745	550	108	174	7	16	15.4	44	30.9
Bendigo	Vic	240	75	71	66	0.640	1.613	489	88	154	5	13	14.4	44	29.9
Deniliquin	NSW	250	71	64	60	0.394	0.989	374	85	109	5	10	15.6	46	31.5
Forbes	NSW	245	69	63	57	0.479	1.017	556	158	134	8	13	16.8	45	33.1
Wagga Wagga	NSW	243	73	66	61	0.487	1.189	489	110	133	7	12	16.4	45	33
West Wyalong	NSW	226	74	66	62	0.427	1.010	466	124	113	7	12	16.4	45	32.8
Total		1744	496	459	422
Minimum						0.394	0.989	374	85	109	5	10	13.7	43	28.8
Maximum						0.740	1.797	556	158	176	8	16	16.8	46	33.1
Trial site	WA					0.796	2.445	587	42	307	2	24	15.2	45	30.6

¹ Ratio precipitation to potential evaporation (pan, free-water surface). Vic = Victoria; NSW = New South Wales; WA = Western Australia; DBH = diameter at breast height; Max. = maximum; Abs. = absolute.

). The site was wetter than original provenance locations, though intrannual variability differed with wetter winters and drier summers (Table 1). Seedlings were planted using a randomised complete block design with four replicates, with families nested by provenance within a replicate. Each replicate contained 10 open-pollinated progeny from 10 mother trees (families) from each provenance, except for the Benalla provenance which had 11 families.

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).

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)

2.2. Phenotypic Analysis

2.2.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]:

$$y=Xb+Zu+e$$

(1)

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²) was calculated as:

$$h^2=\frac{V_A}{V_P}=\frac{2.5{\sigma }_f^2}{{\sigma }_f^2+{\sigma }_{residual}^2}$$

(2)

where $V_A$ is the additive genetic variance, $V_P$ the total phenotypic variance, ${\sigma }_f^2$ the family-level variance and ${\sigma }_{residual}^2$ 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]:

$$Q_{ST}=\frac{{\sigma }_p^2}{{\sigma }_p^2+2V_A}=\frac{{\sigma }_p^2}{{\sigma }_p^2+2\left({\sigma }_f^2\times 2.5\right)}$$

(3)

where ${\sigma }_p^2$ is the genetic variance between populations and $V_A$ and ${\sigma }_f^2$ 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 χ² distribution; specifically

$$\frac{\left(n-1\right){\widehat{F}}_{ST}}{{\overline{F}}_{ST}} \approx {\chi }^2 \left(n-1\right)$$

(4)

where ${\widehat{F}}_{ST}$ = per locus F_ST, ${\overline{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 χ² with df = ${\overline{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).

2.2.2. 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:

$$r_{1,2}=\frac{co{v}_{1,2}}{\sqrt[2]{{\sigma }_1^2\times {\sigma }_2^2}}$$

(5)

where $co{v}_{1,2}$ is the covariance between trait 1 and trait 2 and ${\sigma }_1^2$ and ${\sigma }_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.

2.2.3. 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].

2.3. Genotype-Phenotype Analysis

2.3.1. 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² 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² = 0.003 ± s.d. 0.005), nor within adaptive or neutral SNP subsets (r² = 0.004 ± 0.006 and r² = 0.004 ± 0.005, respectively; linkage (as measured by correlation, r²) was calculated using the ‘r’ function in PLINK v 1.90b3p [69]).

2.3.2. 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.

2.4. 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).

3. Results

3.1. 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²) ranged from 0.106 to 0.373 (

Table 2. Heritability (h²) and between provenance differentiation (Q_ST) for nine quantitative traits measured in Eucalyptus microcarpa (Grey box). SE = one standard error.

Trait		Abbreviation	Units	Family Level		Provenance Level
				h2	SE	QST	SE
Growth traits
	Diameter at breast height	DBH	cm	0.318	0.226	0.094	0.065
	Height	Height	m	0.210	0.137	0.431	0.190	*
	Size Ratio	Size Ratio	height:DBH	0.106	0.080	0.165	0.162
Leaf traits
	Leaf area	Area	cm2	0.155	0.169	0.467	0.220	*
	Leaf length	Length	cm	0.144	0.100	0.483	0.226	*
	Leaf weight	Weight	g	0.235	0.102	0.224	0.145
	Leaf thickness	Thickness	mm	0.193	0.098	0.205	0.147
	Specific leaf area	SLA	mm2 mg−1	0.373	0.116	0.182	0.115
	Leaf density	Density	mg mm−3	0.000	-	-	-

* = Q_ST > F_ST. Lower bound SE (Q_ST – SE) significant (p < 0.05) against χ² distribution with df = mean neutral F_ST (3.3%).

). 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 (${\overline{F}}_{ST}$ ‘neutral’ SNPs = 0.033), suggesting these traits are potentially under divergent selection (Table 2; Figure S1). Low, neutral ${\overline{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).

3.2. 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. Family-level (upper right) and provenance-level (lower left) genetic correlations between traits in Eucalyptus microcarpa (Grey box). Note leaf density excluded due to no additive variance found.

		Growth Traits			Leaf Traits
		DBH	Height	Size Ratio	Area	Length	Weight	Thickness	SLA
Growth traits
	DBH		0.88 *	−0.97 *	−0.20 *	−0.08	−0.03	−0.18	−0.15
	Height	0.93 *		−0.75 *	−0.30 *	−0.56 *	0.16	−0.34 *	−0.35 *
	Size Ratio	0.35	0.81		−0.35	−0.16	−0.55 *	−0.42 *	−0.15
Leaf traits
	Area	0.76	0.53	0.04		0.72 *	0.97 *	0.46 *	−0.36 *
	Length	0.40	0.32	0.38	0.79 *		0.84 *	0.38 *	−0.36 *
	Weight	0.92 *	0.56	−0.01	0.95 *	0.59		0.70 *	−0.77 *
	Thickness	−0.37	−0.22	0.22	−0.74	−0.50	−0.76		−0.98 *
	SLA	0.16	0.19	0.06	0.70	0.93 *	0.34	−0.38

* q < 0.1 (accounting for multiple testing).

). 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).

3.3. Climate Associations with Quantitative Trait Variation

Significant heritable genetic variation (h²) 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 4. Top climate associations (p_env < 0.05) for leaf and growth trait principal component (PC) axes in Eucalyptus microcarpa (Grey box), calculated (a) without and (b) with adjustment for spatial population structure (latitude and longitude). Note p_env = significance of environmental component in the linear model; q_env = accounts for multiple testing.

Trait		Environment	(a) Trait ~ Environment					(b) Trait ~ Lat. + Long. + Environment
			Adj. r2	Assoc.	Fenv [1,5]	penv	qenv	Adj. r2	Assoc.	Fenv [1,3]	penv	qenv
Growth traits
	PC2	Warmest period max. temp.	0.59	+	9.7	0.026	0.13	0.56	+	0.1	0.831	0.811
	PC2	Mean annual temp.	0.53	+	7.7	0.040	0.15	0.56	−	0.1	0.835	0.811
Leaf traits
	PC2	Max. abs. mean max. temp.	0.76	+	20.1	0.007	0.13	0.69	+	7.4	0.073	0.400
	PC2	Mean annual aridity	0.73	−	17.1	0.009	0.13	0.71	−	8.4	0.063	0.400
	PC2	Winter precipitation	0.69	−	14.6	0.012	0.13	0.81	−	13.8	0.034	0.400
	PC2	Max. aridity	0.63	−	11.1	0.021	0.13	0.58	−	4.7	0.119	0.575
	PC2	Warmest period max. temp.	0.62	+	10.6	0.022	0.13	0.70	+	7.9	0.067	0.400
	PC2	Mean annual temp.	0.59	+	9.5	0.027	0.13	0.74	+	9.4	0.055	0.400
	PC3	Summer precipitation	0.52	−	7.6	0.040	0.15	0.26	−	2.5	0.215	0.811
	PC2	Summer precipitation	−0.17	+	0.1	0.732	0.61	0.88	−	23.1	0.017	0.385
	PC2	Driest period prec.	−0.17	+	0.1	0.742	0.61	0.87	−	22.7	0.018	0.385
	PC2	Annual precipitation	0.13	−	1.9	0.231	0.43	0.77	−	11.4	0.043	0.400

Adj. r² = adjusted r² of linear model; Lat. = Latitude; Long. = Longitude; Assoc. = association (+ = positive; − = negative).

a; 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.

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).

3.4. 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. Significant SNP-trait associations (p < 0.05) in Eucalyptus microcarpa (Grey box)and possible gene function of SNPs based on Arabidopsis thaliana orthologues (TAIR10) of predicted E. grandis genes (information from E. grandis v1.1 genome annotation). Mixed linear models (MLM) accounted for kinship using either a kinship matrix based on individual kinship (indv. kinship) or family-average kinship (avg. kinship).

		MLM (Indv. Kinship)		MLM (Avg. Kinship)		Eucalyptus grandis (v1.1) Gene Information (+/−2000 bp)		Best TAIR10 Gene Orthologue
Trait Marker		p	r2	p	r2	Name	Gene Effect	Name	Symbol	Definition
DBH
	3:59841756	0.043	0.019	0.049	0.019	Eucgr.C03147	synonymous	AT3G13980.1
Size Ratio
	10:29282238	0.023	0.012	0.015	0.014	Eucgr.J02333	synonymous	AT3G27150.1		Galactose oxidase/kelch repeat superfamily protein
	10:29282238	0.023	0.012	0.015	0.014	Eucgr.J02334	upstream	AT5G14170.1	CHC1	SWIB/MDM2 domain superfamily protein
	2:58822368	0.060	0.009	0.044	0.010	Eucgr.B03399	downstream	AT1G27150.1		Tetratricopeptide repeat (TPR)-like superfamily protein
Leaf area
	4:30801453	0.069	0.013	0.037	0.016	Eucgr.D01681	synonymous, intron	AT1G64660.1	ATMGL, MGL	methionine gamma-lyase
	11:4085447	0.052	0.016	0.044	0.017	Eucgr.K00355	downstream	AT5G38280.1	PR5K	PR5-like receptor kinase
Leaf length
	2:63702271	0.036	0.016	0.002	0.030	Eucgr.B03985	missense	AT5G49620.1	AtMYB78, MYB78	myb domain protein 78
	4:30801453	0.012	0.021	0.004	0.026	Eucgr.D01681	synonymous, intron	AT1G64660.1	ATMGL, MGL	methionine gamma-lyase
	5:5009176	0.030	0.016	0.130	0.009	Eucgr.E00527	missense	AT4G23020.1
Leaf weight
	5:4625458	0.028	0.017	0.027	0.017	Eucgr.E00491	upstream	AT3G60890.2	ZPR2	protein binding
Leaf thickness
	5:4625458	0.006	0.023	0.007	0.023	Eucgr.E00491	upstream	AT3G60890.2	ZPR2	protein binding
SLA
	5:4625458	0.024	0.018	0.017	0.020	Eucgr.E00491	upstream	AT3G60890.2	ZPR2	protein binding
	11:4085447	0.020	0.021	0.015	0.022	Eucgr.K00355	downstream	AT5G38280.1	PR5K	PR5-like receptor kinase
Leaf density
	5:4625458	0.102	0.010	0.042	0.014	Eucgr.E00491	upstream	AT3G60890.2	ZPR2	protein binding
	6:39617441	0.086	0.014	0.044	0.018	Eucgr.F02999	synonymous	AT1G22610.1		C2 calcium/lipid-binding plant phosphoribosyltransferase family protein

). 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.

3.5. 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.

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).

4. 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.

4.1. 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² = 0.06–0.37, height h² = 0.13; [22,35]), although estimates for leaf traits were low relative to other studies (e.g., leaf length h² = 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.

4.2. 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].

4.3. 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.

5. 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.