Directional Selection on Tree Seedling Traits Driven by Experimental Drought Differs Between Mesic and Dry Populations

We evaluated population differences and drought-induced phenotypic selection on four seedling traits of the Australian forest tree Eucalyptus pauciflora using a glasshouse dry-down experiment. We compared dry and mesic populations and tested for directional selection on lamina length (reflecting leaf size), leaf shape, the node of ontogenetic transition to the petiolate leaf (reflecting the loss of vegetative juvenility), and lignotuber size (reflecting a recovery trait). On average, the dry population had smaller and broader leaves, greater retention of the juvenile leaf state and larger lignotubers than the mesic population, but the populations did not differ in seedling survival. While there was statistical support for directional selection acting on the focal traits in one or other population, and for differences between populations in selection gradient estimates for two traits, only one trait—lamina length—exhibited a pattern of directional selection consistent with the observed population differences being a result of past adaptation to reduce seedling susceptibility to acute drought. The observed directional selection for lamina length in the mesic population suggests that future increases in drought risk in the wild will shift the mean of the mesic population toward that of the dry population. Further, we provide evidence suggesting an early age trade-off between drought damage and recovery traits, with phenotypes which develop larger lignotubers early being more susceptible to drought death. Such trade-offs could have contributed to the absence of population mean differences in survival, despite marked differentiation in seedling traits.


INTRODUCTION
Global warming is leading to major shifts in the distribution of plant and animal species globally (Chen et al., 2011). Sessile species, such as forest trees, are among the most vulnerable owing to their often poor dispersal abilities and long life histories (Petit and Hampe, 2006). Already, extreme climate events causing heat and drought stress have been linked directly (Allen et al., 2010;Matusick et al., 2013) and indirectly (Buotte et al., 2016) to forest dieback in many parts of the world. Ongoing climate change is expected to exacerbate such stresses (Allen et al., 2010;Schwalm et al., 2017), and it has long been argued that the rate of climate change may exceed the adaptive capacity of tree populations (Davis and Shaw, 2001;Jump and Peñuelas, 2005;Aitken et al., 2008).
Tree adaptation to major environmental gradients such as aridity is polygenic and genome-wide (Kramer et al., 2015;Jordan et al., 2017;Steane et al., 2017), involves diverse functional traits (Alberto et al., 2013;Kremer et al., 2014), and potentially also genetic constraints (Bourne et al., 2017;Costa e Silva et al., 2020). Detecting adaptive changes in local populations of longlived organisms, such as forest trees, is difficult due to numerous factors (Alberto et al., 2013), including long-life cycles (Petit and Hampe, 2006), ontogenetic (Brunner et al., 2016) and plastic (Nicotra et al., 2010;McLean et al., 2014) changes in phenotype, superimposed on the normal selective filtering of mal-adapted inbred progeny during stand development (Koelewijn et al., 1999;Costa e Silva et al., 2011;Griffin et al., 2019). Molecular markers are often used to detect signatures of adaptive genetic variation in tree species (Keller et al., 2012;Alberto et al., 2013;Dillon et al., 2014;Pluess et al., 2016;Frachon et al., 2018) and, as the effects of putative adaptive markers/SNPs become validated, there will be increasing opportunity for their use to monitor adaptive changes in forest tree populations following selection. However, traditional analysis of phenotypic selection and quantitative genetic studies of functional traits continue to provide valuable insights into the effects of natural selection, and have been applied in many animal and plant species (Kingsolver et al., 2001(Kingsolver et al., , 2012Shaw and Etterson, 2012;Merilä and Hendry, 2014). Despite this, there is a paucity of empirical studies in forest tree populations that address the adaptive value of functional traits through their impact on fitness components (Gómez, 2004;Alberto et al., 2013;Ramírez-Valiente et al., 2014Costa e Silva et al., 2018;Shaw, 2018, 2019), as is the case for most long-lived species (Monica and Lauren, 2003). Lande and Arnold (1983) developed a least-squares multiple regression framework for the analysis of phenotypic selection based on the relationship between relative fitness (i.e., the fitness of an individual divided by the mean fitness of the population) and phenotypic measures for a set of traits. This approach assumes a causal fitness-trait relationship, such that a trait (or a trait combination) is expected to be under selection when significantly influencing fitness. The partial regression coefficients thus represent selection gradients that quantify the direct effects of trait variables on relative fitness, and the complete model with linear, quadratic and cross-product terms provides a quadratic approximation to the individual fitness surface (Phillips and Arnold, 1989). In particular, by characterizing directional selection, the linear selection gradient has a direct application in the prediction of evolutionary change between generations, under the assumption of multivariate normal distribution of breeding values (Lande, 1979). Fitness components used to evaluate selection commonly involve discrete variables (e.g., viability and fecundity measured by survival and the number of offspring), and generalized linear models have been described for the estimation of selection gradients with non-normal distributions of fitness responses (Janzen and Stern, 1998;Morrissey and Sakrejda, 2013).
Numerous studies have shown significant directional selection in wild populations (Kingsolver and Diamond, 2011;Kingsolver et al., 2012;Siepielski et al., 2013). Such selection can vary with time/seasons (Siepielski et al., 2009;Steele et al., 2011;Ramírez-Valiente et al., 2014Shaw, 2018, 2019) and environment (Dudley, 1996;Etterson, 2004;Heschel and Riginos, 2005;Donovan et al., 2009;Lau et al., 2014;Warwell and Shaw, 2018), and this variation drives population divergence and local adaptation within a species (Siepielski et al., 2013;Hendry, 2017). However, there are often multiple phenotypic solutions to the same adaptive requirement (e.g., drought survival-McDowell et al., 2011;Choat et al., 2018;Klein et al., 2018) and increasing evidence from common-garden trials that the focal traits under selection can be population specific (Dudley, 1996;Etterson, 2004;Heschel and Riginos, 2005;Donovan et al., 2009). Understanding population differences in the focal traits under selection following the same selection event can allow insights into functional traits that may be underlying local adaptation (Dudley, 1996;Monica and Lauren, 2003) and how these vary within a species. Such understanding is a key step toward evaluating the effects of selection due to, for example, climate change (Hoffmann and Sgro, 2011;Alberto et al., 2013;Merilä and Hendry, 2014), including outcomes for different populations or parts of the species distribution.
The present study focuses on the Australian tree species Eucalyptus pauciflora Sieber ex Spreng.. Eucalyptus pauciflora has one of the broadest altitudinal ranges of any eucalypt species (i.e., growing from sea-level to near the alpine tree line on both continental Australia and the island of Tasmania), and is one of the main tree species used for ecological restoration in the dry midlands of Tasmania (Prober et al., 2016;Bailey et al., 2021). Climate change is predicted to result in most Tasmanian E. pauciflora populations being outside the contemporary climate envelope of the species and at risk of maladaptation by 2080 (Harrison, 2017). This risk of maladaptation mainly reflects the predicted increasing temperatures and evapotranspiration across Tasmania (Corney et al., 2010;Harrison et al., 2017), and thus increased risk of drought stress. Both water deficit and high temperatures contribute to drought stress (Park Williams et al., 2013;Duan et al., 2014;Mitchell et al., 2014) and could, independently or combined, drive selection on plant traits (Costa e Silva et al., 2019). The challenge in understanding adaptation is to unravel these multiple effects. Here we focus on just a single factor-water deficit, which has both plastic (Cano et al., 2014;Duan et al., 2014;McKiernan et al., 2016) and genetic-based (Ammitzboll et al., 2020) selective impacts on eucalypt seedlings.
As with other forest tree studies on drought-induced selection (e.g., Warwell and Shaw, 2019), we focused on the seedling life stage as: (i) there is usually significant seedling mortality in wild and planted eucalypt populations, with water deficit being the most common cause (Battaglia and Reid, 1993;Stoneman, 1994;Tozer and Bradstock, 1998); (ii) the phenotypic and genetic variation in morphological traits at this life stage has been wellcharacterized in E. pauciflora (Gauli et al., 2015); and (iii) there are signals that variation in some seedling morphological traits of E. pauciflora is adaptive (Gauli et al., 2015;Costa e Silva et al., 2019). We specifically aimed to: (1) assess the difference between dry and mesic populations in four seedling traits that potentially may influence the plant's response to acute drought; (2) evaluate the consequences of the variation in the focal seedling traits on fitness in each population; and (3) estimate whether the relationship between fitness and the focal traits differs between populations. We hypothesized that, if selection acting on the studied traits via seedling survival was driven by acute drought stress in wild populations, then an imposed artificial drought would elicit a shift in character states of the mesic population toward those of the dry population.

Populations Studied
The current study used seedlings from two populations of E. pauciflora, which were grown in a common-garden glasshouse trial subject to an artificial drought. These two lowland Tasmanian populations were from Brushy Lagoon and Ross, and were from mesic and dry parts of the climatic range of E. pauciflora on the island, respectively. These populations were selected to accentuate rainfall differences and minimize the differences in other environmental factors. The mean annual precipitation at Ross is nearly 50% lower than at Brushy Lagoon, but Ross also has a slightly higher mean maximum temperature of the warmest period than Brushy Lagoon (Table 1). These populations are located over 90 km apart, well-exceeding the distance of significant genetic spatial autocorrelation in the species (i.e., 27 km- Gauli et al., 2013). The populations have similar high outcrossing rates (Gauli et al., 2013) and inbreeding levels in their open-pollinated progeny (Gauli et al., 2014), but differ markedly in numerous seedling traits (Gauli et al., 2015).

Drought Experiment
Each population was represented by 12 families corresponding to 12 mother-trees randomly selected within native stands and situated by more than two canopy heights apart. The seedlings of each family were obtained from open-pollinated seeds collected randomly from near the middle of the canopy of the mothertree. Seedlings from these families were grown in soil-filled pots (9 cm × 9 cm × 18 cm) and, when the plants were 14-week-old, the pots were arranged into an experimental design for imposing drought stress under glasshouse conditions. The experiment was conducted in two glasshouses, though a large portion of window panels were removed so that seedlings were exposed to ambient air conditions. Temperature and relative humidity were similar between the two glasshouses, and temperature reflected external ambient air temperature over the duration of the experiment (Supplementary Figure 1). Over the course of the dry-down period, the average daily maximum temperature inside the glasshouses was 16.9 and 16.0 • C, and the average daily minimum temperatures were 6.5 and 6.8 • C. Daily maxima did not exceed 27 • C and daily minima did not go below −1 • C. Humidity ranged from an average daily maximum of 88.5 and 90.9% to an average daily minimum of 51.6 and 58.5% in the two glasshouses, respectively. The experimental layout was effectively a split-plot design comprising four replicates (blocks), two in each glasshouse. Each replicate included two main plots (hereafter denoted as population plots), one of the Brushy Lagoon families and one of the Ross families, and families were randomized as single-tree plots within these population plots.
Overall, there were 387 and 386 seedlings in the Brushy Lagoon and Ross populations, respectively, available for analysis. Across the experiment, the sizes of the population plots varied between 92 and 101 seedlings in Brushy Lagoon, and between 93 and 100 seedlings in Ross. The family sizes across all four replicates ranged from 22 to 38 seedlings in Brushy Lagoon, and from 19 to 38 seedlings in Ross, with an average of 32 seedlings per family in both populations. The degree of unbalance in family sizes was moderate: the variation of family sizes as measured by the coefficient of variation was 16% in both populations.
The water deficit treatment was applied using floral foam, in a manner adapted from Fernández and Reynolds (2000) and Marchin et al. (2020). Each population plot consisted of a water basin containing a layer of floral foam (Oasis Ideal Standard Super 60; Apack Pty Ltd., Dandenong South, Australia). Potted seedlings were placed on the foam, with a small amount of fine washed sand in drainage holes to ensure effective contact with foam. Water basins were topped up twice daily to within approximately 10 cm from the top of the foam, enabling selfwatering by wicking, though with reduced water availability for plants (see Marchin et al., 2020). Exposed foam, above the top of the water basins, was wrapped in plastic to minimize water loss through evaporation.
After placing seedlings in the experimental design, they were left to acclimate for 10 days, over which time morphometric measurements were taken from the seedlings (see below). After 10 days, no further water was added to the water basins, initiating the drought component of the experiment. At the end of the experiment 67 days later, seedling survival was scored as a binary outcome, based on the following assessments: a seedling was classified as dead (scored as 0) when all leaves, including laterals, and the apical meristem had dried out; otherwise it was classified as alive (scored as 1). As scored, seedling survival refers to observable "canopy" mortality immediately following termination of the dry-down experiment, and does not account for the possibility of differential recovery if rewatering had occurred. This measure of seedling survival was used as a fitness component for individuals in the analysis of phenotypic selection.

Trait Measurements
Phenotypic trait measurements were taken from the seedlings during the 10-day acclimation period. At this time, a single leaf was sampled from the third pair of leaves up the main stem (cotyledonary node = 0), and individually digitally photographed with a constant size scale for the measurement of leaf dimensions (see below). Following the removal of the third-node leaf, these seedlings were also measured for the lignotuber diameteri.e., the diameter (mm) across the lignotubers developed in the axils of the cotyledons-and for the stem diameter (mm) perpendicular to the axis of maximum lignotuber diameter. The seedling height (cm) and the node number of the first petiolate leaf (counted from the cotyledonary node) were also recorded. The node number of the first petiolate leaf is a trait that reflects the ontogenetic transition from sessile (juvenile) to petiolate (mature) leaves in Eucalyptus spp., and is hereafter referred to as "vegetative juvenility" (i.e., the higher the node number, the later the petiole development is in nodes, and thus the more the seedling is vegetatively juvenile). Leaf dimensions were manually measured from the digital images of the seedling leaves, and included the lamina length (mm) and lamina width (mm). In the present study, seedling leaf morphology, lignotuber size and vegetative juvenility were represented by four functional traits-lamina length, leaf shape index, adjusted lignotuber size and vegetative juvenility (see below). Photos of lignotubers and vegetative juvenility in seedlings of E. pauciflora are provided in Figure 3 of

Data Analyses
The data analyses entailed two main components. Firstly, population mean differentiation and variance estimation were addressed for the focal functional traits and seedling survival. Then, phenotypic selection analysis was pursued within populations to measure directional selection acting on the focal traits through seedling survival, and differences between populations in estimated selection gradients for a trait were subsequently evaluated. To account for the inherent positive allometric relationship between lamina length and lamina width, the lamina width measurements were adjusted for seedling differences in lamina length in the data analyses. The adjusted lamina width thus provides a measure that reflects leaf shape, and is hereafter denoted as "leaf shape index" (with increasing values meaning broader leaves). The lignotuber diameter measurement is also inherently dependent on stem diameter and, following Costa e Silva et al. (2019), was also adjusted in the data analyses to provide a measure independent of stem diameter, hereafter referred to as "adjusted lignotuber size." Finally, in addition to the focal traits, seedling height and stem diameter were modeled as "nuisance variables" in the phenotypic selection analysis. This was done to account for the effects of initial plant size on drought-induced mortality, as the seedlings were growing in pots with a defined soil volume, and bigger plants are expected to inherently use more water (Ammitzboll et al., 2020). The SAS 9.4 software (SAS, 2017) was used in the data analyses. In particular, linear mixed models were estimated with the procedure MIXED, and generalized linear mixed models were estimated with the procedures GLIMMIX and NLMIXED.

Linear Mixed Model for the Continuous Variables (Lamina Length, Leaf Shape Index and Adjusted Lignotuber Size)
The following linear mixed model (LMM) was fitted for the continuous variables: where y is a n x 1 vector of observations (n = total number of seedlings in the experiment); β is a p × 1 vector of fixed-effect parameters; b 1 and b 2 are q × 1 vectors of random-effect parameters for families within the Brushy Lagoon and Ross populations, respectively; X is a design matrix relating y to the parameters in β; Z 1 and Z 2 are design matrices relating y to the random effects in b 1 and b 2 , respectively; e is a n × 1 vector of residuals of the observations. The vector β included the intercept, and effects for population, block and populationby-block interaction for all traits. The analyses of leaf shape index and adjusted lignotuber size involved the adjustment of the lamina width and lignotuber diameter measurements using lamina length and stem diameter as covariates, respectively. In these cases, β included also the corresponding coefficients for the regression parameter and its interaction with population. The family effects and residuals were modeled by considering a heterogeneous variance structure at the population level. In this sense, it was assumed that are family variances in Brushy Lagoon and Ross, respectively; I q 1 and I q 2 are identity matrices of dimension equal to the number of families in Brushy Lagoon (q 1 ) and Ross (q 2 ). Conditional on the family random effects, the residuals were assumed e 1 e 2 ∼ MVN(0, R), , where σ 2 e 1 and σ 2 e 2 are residual variances in Brushy Lagoon and Ross, respectively; I n 1 and I n 2 are identity matrices of dimension equal to the number of seedlings in Brushy Lagoon (n 1 ) and Ross (n 2 ); it was further assumed that the residuals and random effects were independent. Visual inspection of histograms and quantile-quantile plots of conditional studentized residuals showed residual distributions conforming reasonably well to a normal distribution for all traits. In addition, the residuals did not appear to co-vary with random effects, as generally indicated by small and statistically nonsignificant correlations estimated within either population. The (co)variance structure of y was defined by V = ZGZ + R, where Z = Z 1 0 0 Z 2 ; through G, the V-matrix captures the covariance of the observations within families, which will be reflected in the estimation of the fixed-effect parameters and their standard errors (Raudenbush and Bryk, 2002; see next paragraph). The variances in G and R were estimated by restricted maximum likelihood (REML; Patterson and Thompson, 1971). Following the estimation of variance components, fixed-effect parameters were estimated by generalized leastsquares, i.e.β = (X V −1 X) −1 X V −1 y, with the associated model-based sampling (co)variance matrix being given by Var(β) = (X V −1 X) −1 . This matrix was corrected by using the Kenward-Roger (KR) adjustment, to provide estimates of standard errors ofβ that account for potential finite-sample bias and the uncertainty in the estimation of G and R (Kenward and Roger, 2009). Based on a corrected Var(β) matrix, the Roger (1997, 2009) approximation to compute denominator degrees of freedom was applied to improve statistical inference about the fixed effects. For sample sizes comparable to the number of families and the within-family sizes in our study, simulation work with LMMs using REML indicated that applying the KR method to fixed-effect parameters produced accurate estimates of standard errors, and maintained Type-I error rates and confidence interval coverage close to the specified (α = 0.05) nominal level (Ferron et al., 2009;Baldwin and Fellingham, 2013;McNeish and Stapleton, 2016a; see also Supplementary Methods 1).
Statistical inference about the fixed-effect model terms was based on F-tests, which indicated that the interactions of population with blocks, and with either lamina length or stem diameter used as covariates (in the LMM for leaf shape index or adjusted lignotuber size, respectively), were not statistically significant (p > 0.05). Thus, a final LMM without the interaction terms was then fitted for all traits. One-sided likelihood ratio tests (Self and Liang, 1987) were applied to assess whether estimated family variances were significantly greater than zero. In preliminary analyses, a one-sided likelihood ratio test was also used to test the variance estimate associated with a random term modeling a block-by-family interaction: as no significant effects were detected (suggesting an efficient randomization of seedlings within blocks), the interaction term was excluded in the LMM defined in Equation (1). The population least-squares mean (LS-Mean) for a trait was calculated with a linear function ofβ estimates, i.e. Lβ, where L is the coefficient matrix pertaining to the LS-Mean. The standard error of a LS-Mean estimate was approximated by the square-root of LVar(β)L . The KR approximation for computation of degrees of freedom was also applied to obtain the critical values for the student's t-distribution used for estimating the 95% confidence interval of a LS-Mean.

Generalized Linear Mixed Model for the Discrete Variables (Seedling Survival and Vegetative Juvenility)
The following generalized linear mixed model (GLMM) was fitted for the discrete variables: where g(.) is a link function relating the conditional expectation [i.e. E(y | X, b 1 , b 2 ) = µ, where µ is the mean vector] of the response variable to the linear predictor η = Xβ + Z 1 b 1 + Z 2 b 2 ; the vectors and matrices in η were defined as in Equation (1) (except that η does not contain a definite error term, nor continuous covariates in X as in the LMMs for leaf shape index and adjusted lignotuber size); for b 1 and b 2 , the model specification with heterogeneous variances and the distributional assumptions were the same as in the LMM. Conditional on the random effects, the observations were assumed to be independently distributed, with Bernoulli and Poisson distributions for the binary (seedling survival) and count (vegetative juvenility) variables, respectively. For the binary variable, E(y | X, b 1 , b 2 ) is also the conditional probability of seedling survival [i.e. P(y = 1 | X, b 1 , b 2 )]. The probit link function was used to model survival, such that and (.) is the standard normal cumulative distribution function. For vegetative juvenility, g(.) was the logarithmic link function, and thus where exp denotes the exponential function. The choice of the probit link function for seedling survival was based on the conclusions obtained by Hahn and Soyer (2005) with a randomintercept GLMM [i.e., akin to the model described in Equation (2)]. Previous GLMM analyses conducted separately in each population showed that the variance of the Pearson residuals was close to one for seedling survival, but was clearly smaller than one for vegetative juvenility, thus indicating that the GLMM for the latter trait was under-dispersed relative to the Poisson distribution. Consequently, as an alternative to the standard Poisson model, we used a quasi-Poisson GLMM for vegetative juvenility, under which the equi-dispersion assumption can be relaxed (Wedderburn, 1974;Gardner et al., 1995). The quasi-Poisson GLMM considered heterogeneity of dispersion effects at the population level, and further details are provided in Supplementary Methods 2.
Parameter estimation under the GLMMs fitted for the discrete variables was based on the residual pseudo-likelihood (RSPL) method, which enables extension of the concept of REML variance component estimation to the non-normal case (Wolfinger and O'Connell, 1993). Estimates of β were obtained by generalized least-squares, and the KR approximations described above were used to improve the model-based standard errors ofβ and the statistical inference concerning fixed-effect model terms. For conditions covering those in the current study (i.e., sample sizes, survival prevalence, and mean vegetative juvenility), simulation work with cross-sectional clustered data indicated that RSPL together with the KR adjustments performed better than maximum likelihood (ML) methods, regarding the accuracy of estimates of the random-intercept variance (e.g., family variance) and the coverage of fixed-effect confidence intervals for binary and count variables (McNeish, 2016(McNeish, , 2019; see also Supplementary Methods 1).
Results from F-tests indicated that the fixed-effect term for the population-by-block interaction was significant (p < 0.05) for seedling survival, but not for vegetative juvenility; thus, the interaction term was excluded in a final GLMM for the latter trait only. The population LS-Mean for a variable was estimated on the linked scale by Lβ and then converted to the original scale via the inverse link transformation. The standard error of a LS-Mean on the original scale was approximated by φ(Lβ) . stde(Lβ) and exp(Lβ) . stde(Lβ) for seedling survival and vegetative juvenility, respectively, where: φ(Lβ) is the probability density function of the standard normal distribution evaluated at Lβ; and stde(Lβ) is the estimated standard error of the LS-Mean on the linked scale. The 95% confidence interval of a LS-Mean on the original scale was obtained by applying the inverse link transformation to the confidence limits on the linked scale.
Under the RSPL method, the absence of a true log-likelihood function for the overall optimization process precludes the comparison of nested models using likelihood ratio tests (e.g., Hedeker and Gibbons, 2006;McNeish, 2016). Nevertheless, for comparison with RSPL estimates, family variances were also estimated for seedling survival and vegetative juvenility using ML, which enabled their statistical significance to be assessed by one-sided likelihood ratio tests. In this context, the Conway-Maxwell Poisson (CMP) distribution was assumed in the GLMM fitted for vegetative juvenility, since the estimation of a Rside multiplicative dispersion parameter as in the quasi-Poisson model (see Supplementary Methods 2) is not permitted with ML methods (e.g., Capanu et al., 2013). The CMP distribution is a two-parameter generalized form of the Poisson distribution that allows the modeling of empirical distributions of count variables that are over-or under-dispersed relative to the Poisson distribution (Shmueli et al., 2005;Guikema and Goffelt, 2008). Details on the parameterization of the CMP distribution, and on the GLMM for joint modeling the mean response and dispersion of vegetative juvenility, are given in Supplementary Methods 3. The GLMM fitted for seedling survival was defined as in Equation (2). ML estimation of the model parameters was based on an adaptive Gaussian quadrature, which entails a numerical integration method to approximate the marginal log-likelihood of the data (Pinheiro and Bates, 1995); 10 quadrature points were used to integrate the likelihood function over the randomeffects distribution.
For both the LMMs and GLMMs, the intra-class correlation (ICC) coefficient was calculated as the proportion of the total variance attributable to the family variance, and thus was used to measure the relative importance of the family variance for each variable. In particular, for vegetative juvenility and seedling survival, the ICCs were estimated as described by Nakagawa et al. (2017)

Analysis of Phenotypic Selection
In the present study, we focused on the estimation of linear selection gradients to measure the direction and strength of directional selection acting on traits in each population, and then we evaluated the difference between populations in selection gradient estimates for a given trait. The analysis of phenotypic selection assessed the fitness consequences of variation in four traits-lamina length, leaf shape index, adjusted lignotuber size and vegetative juvenility-using seedling survival as a fitness component. In addition to the target traits, seedling height and stem diameter were modeled as "nuisance variables." Prior to analyses, the data were pre-processed as follows: (i) linear leastsquares regression was used to pre-adjust lamina width and lignotuber diameter in a population for seedling differences in lamina length and stem diameter, respectively, and the adjusted observations (obtained from the regression residuals) were then mean standardized-these data thus refer to leaf shape index and adjusted lignotuber size; and (ii) the remaining variables were centered at their means and then mean-standardized in a population. Mean standardization has been suggested for comparing the strengths of selection estimates on a similar scale (Hereford et al., 2004), and scaling in relation to the mean is both permissible and meaningful for the independent variables used in the current analysis (Houle et al., 2011).
The analysis of phenotypic selection was undertaken separately in each population, with the GLMM fitted for seedling survival including: a vector β with the intercept, block effects, and the linear regression coefficients for the focal and nuisance variables; and a vector b of random effects for families within the population, assuming b ∼ MVN(0, G), with G = I q σ 2 b . Akin to the GLMM described in Equation (2), g(.) was the probit link function, and the model parameters were estimated by the RSPL algorithm.
As opposed to least-squares multiple regression, the marginal effect of a predictor variable (i.e., the partial effect of a unit change in the predictor on the probability of seedling survival) cannot be directly obtained from the estimates of β in the GLMM. In this context, for a focal variable, the directional selection gradient can be obtained from the average gradient of the estimated fitness surface, which can be approximated by averaging the slope of the fitness surface over the distribution of the phenotypic observations (Janzen and Stern, 1998;Morrissey and Sakrejda, 2013). This procedure was undertaken as follows: (i) the marginal effect at each observation point in the data was estimated for a given focal variable; (ii) the estimated marginal effects were averaged over all observations in a population; and (iii) the average of the marginal effects was standardized to the relative fitness scale (i.e., division by the population mean of seedling survival), yielding a measure that can be quantitatively interpretable as a selection gradient. Following these sequential steps, Supplementary Methods 5 provides the methodological details used to obtain directional selection gradients, based on the model parameters estimated from the random-intercept GLMM fitted to the data of a given population. The resulting selection gradient estimates have a population-averaged interpretation (see Supplementary Methods 5), although this does not apply to the estimates of β from the GLMM, which have a conditional interpretation (e.g., Molenberghs et al., 2002;Gardiner et al., 2009;McNeish, 2016; see also Supplementary Methods 7).
Non-parametric bootstrap was used to obtain standard errors and 95% confidence intervals for the directional selection gradient estimates in a population, and for their differences between populations for a trait. Bootstrap samples were generated while preserving the conditional independence of seedlings within families. Further details on this procedure, and on the estimation of standard errors and 95% confidence intervals, are given in Supplementary Methods 6. Statistical support against a null hypothesis being true (i.e., against the absence of an effect for either a selection gradient in a population or its difference between populations) was provided by a 95% confidence interval not overlapping with zero. Although directional selection was the main focus of the current study, we also computed 95% confidence intervals from bootstrapping for quadratic and correlational selection gradients (with their estimation based on second derivatives approximated at each observation by secondorder central finite differences; e.g., see Supplemental File, point 8.4, in Franklin and Morrissey, 2017). As these confidence intervals overlapped with zero, the effects of non-linear selection were judged to be unimportant, and they will not be further considered here.

RESULTS
A diagrammatic summary of the overall results is provided in Figure 1. Table 2 presents the estimates of LS-Means, and related results from statistical tests of the population effect, for four functional traits and the fitness component, measured in the mesic (Brushy Lagoon) and dry (Ross) populations. Results from statistical inference about all the fixed-effect model terms included in the final LMMs and GLMMs fitted to the continuous and discrete variables are shown in Supplementary Tables 1, 2, respectively. When compared to the mesic population, the LS-Means in the dry population were lower for lamina length, and higher for leaf shape index, adjusted lignotuber size and vegetative juvenility. Population-mean differences were statistically significant at the 0.1% (lamina length, adjusted lignotuber size, and vegetative juvenility) and 1% (leaf shape index) significance levels for the functional traits. The LS-Means for the probability of seedling survival were similar in the mesic and dry populations, and not statistically different (p > 0.05). For discrete variables (vegetative juvenility and seedling survival), LS-Means based on the conditional fixedeffect coefficients from the GLMM (as presented in Table 2) were similar to corresponding estimates calculated from fixedeffect coefficients with a marginal (i.e., population-averaged) interpretation (see Supplementary Methods 7).
REML estimates of family and residual variances are shown in Table 3 for the continuous variables. Based on one-sided likelihood ratio tests, statistically significant (p < 0.05) family variance estimates were found for lamina length, leaf shape index and adjusted lignotuber size in both populations, suggesting significant genetic variation within each population for these traits. The proportion of the total variation explained by the family variance (i.e., ICC) was ≈ 0.2 for lamina length and leaf shape index, but lower for adjusted lignotuber size (i.e., 0.03 and 0.07 in Brushy Lagoon and Ross, respectively; Table 3). Two-sided likelihood ratio tests undertaken to evaluate whether residual variances differed between populations indicated significant variance heterogeneity for lamina length and adjusted lignotuber size (p < 0.001), but not for leaf shape index (p > 0.05). In addition, for the difference between populations in family variances, two-sided likelihood ratio tests did not indicate significant variance heterogeneity at the 5% significance level for any of the traits. However, the latter results may reflect the lack of statistical power to detect variance heterogeneity (see also below) for some traits, given that family variance estimates differed between populations by a factor of two for lamina length and three for adjusted lignotuber size ( Table 3). Table 4 provides the family variances estimated by RSPL and ML for the discrete variables. The ML estimates of family variances were noticeably smaller than their counterparts estimated by RSPL: for Brushy Lagoon and Ross, the differences of the ML estimates relative to the RSPL estimates were 33 and 18% for vegetative juvenility, and 16 and 11% for seedling survival (Table 4). Although the true parameter values are not known, these results are consistent with the downward finitesample bias that may emerge for ML estimates of variance components in GLMMs with discrete variables when the number of clustering units is not large (see references in Supplementary Methods 1). RSPL estimates of family variances were thus considered to be more accurate than their ML counterparts, given that simulation studies (e.g., McNeish, 2016, 2019) using sample sizes comparable to the number of families and their sizes in the current study indicated substantial downward bias in ML estimates of variances for the random-intercept term, and no evidence of noticeable bias for corresponding RSPL estimates. Underestimation of variances will also lead to underestimation of standard errors for fixed-effect parameters, which may result in less conservative hypothesis tests for fixedeffect model terms. Nevertheless, for the examined vegetative juvenility and seedling survival data, the method to estimate model parameters was inconsequential regarding conclusions on the statistical significance of the population effect (although the F-value increased under ML estimation, particularly for vegetative juvenility; Supplementary Table 2).
Based on RSPL estimates of family variances, the ICC ranged from 0.25 to 0.28 for vegetative juvenility, and from 0.02 to 0.06 for seedling survival ( Table 4). As RSPL applies linearization methods to approximate the model rather than the likelihood function, ML estimates were used to conduct onesided likelihood ratio tests for family variances of vegetative juvenility and seedling survival. These tests detected statistically significant (p < 0.05) family variance estimates (and thus genetic variation) for vegetative juvenility in both populations, but only in Ross for seedling survival (Table 4). However, for seedling survival in Brushy Lagoon, the likelihood ratio test is likely to have a low statistical power: with our family sizes, the number of families would need to be larger for detecting a non-null family variance, as suggested by simulation work with binary variables and weak ICCs (Austin and Leckie, 2018). Lack of statistical power may also have affected two-sided likelihood ratio tests pursued to evaluate whether family variances differed between populations, as variance heterogeneity was not found to be statistically significant (p > 0.05), despite the ≈ three times difference between populations in family variances for seedling survival (Table 4). Conversely, under a GLMM using a CMP distribution for vegetative juvenility, two-sided likelihood ratio tests indicated statistically significant (p < 0.001) heterogeneity between populations in dispersion of the observations (see footnotes of Supplementary Table 2). Table 5 shows the estimates of directional selection gradients and their population differences for the focal functional traits. For Brushy Lagoon, statistical support against the null hypothesis of no effect on relative fitness (i.e., as indicated by a 95% confidence interval not overlapping with zero) was found for the selection gradients of lamina length and vegetative juvenility, with directional selection favoring seedlings exhibiting lower values in these traits (i.e., a decrease of 1% from the mean of lamina length and vegetative juvenility was expected to increase relative fitness by 1.5 and 0.3%, respectively). For Ross, statistical FIGURE 1 | Infographic summarizing the differences between the mesic (Brushy Lagoon) and dry (Ross) populations of Eucalyptus pauciflora for the least-squares means and directional selection gradients estimated for four seedling functional traits. Population-mean differences were statistically significantly (p < 0.01) for all focal traits assessed ( Table 2). The directional selection gradient estimates for each population are symbolized by the arrows, which refer to the partial linear effect of a focal trait on seedling canopy (relative) survival evaluated after 67 days of dry-down. The increasing length and width of the arrows is indicative of increasing strength of the directional selection gradient estimates, and their directionality is shown by the +/− signs. In addition, a solid arrow indicates that the bootstrapped 95% confidence interval (CI) for the selection gradient estimate did not overlap with zero, providing statistical support against the null hypothesis of no effect on relative seedling survival ( Table 5). A dashed arrow signifies that the bootstrapped 95% CI for the selection gradient estimate overlapped with zero, being indicative of no statistical support for the existence of fitness consequences of variation in the trait. The strong negative selection gradient estimated for lamina length in the mesic population suggests that a future occurrence of severe drought in the wild will shift the mean of the trait in the mesic population toward that of the dry population. The negative selection gradient estimated for adjusted lignotuber size in both populations suggests an early age trade-off between drought damage and recovery traits, with the canopies of phenotypes which develop larger lignotubers at the seedling stage being more susceptible to drought death. support for the existence of fitness consequences of variation in the traits was provided for the selection gradients of leaf shape index and adjusted lignotuber size, indicating directional selection favoring phenotypes with lower values in these traits (i.e., an increase in 1.2 and 0.9% in relative fitness was expected from decreasing leaf shape index and adjusted lignotuber size by 1% from their means, respectively). Statistical support for the difference between populations in directional selection was observed in the selection gradient estimates for lamina length and leaf shape index, but not for adjusted lignotuber size and vegetative juvenility. Similar conclusions regarding the potential importance of directional selection acting on the traits within populations, and its difference between populations for a trait, were indicated from formal statistical inference about parameter estimates (on the probit scale) obtained from a GLMM fitted for data of the two populations combined (Supplementary Table 3).

DISCUSSION
This study identified significant functional trait differences and differing acute drought-induced phenotypic selection gradients between two populations from contrasting environments of the forest tree E. pauciflora (summarized in Figure 1). Significant genetic variation within populations, and significant mean differences between populations in all four focal traits studied suggest adaptive differences between the mesic and dry populations. While seedling survival did not differ between the two populations, when subject to acute drought stress, the traits under phenotypic selection varied between populations, suggesting different response strategies. Statistical support for the presence of fitness consequences of variation in the focal traits was detected for all traits in one or other population, but only one trait-lamina length-exhibited a pattern of directional selection consistent with a shift in the mean of the mesic population toward The results obtained from the significance test (F-statistic and p-value) undertaken for the population effect in the model are also given. Lamina length (mm); leaf shape index = adjusted lamina width (mm); adjusted lignotuber size = adjusted lignotuber diameter (mm); vegetative juvenility = node number of the first petiolate leaf; seedling survival = survival measured as a binary outcome (where 1 = alive and 0 = dead). Estimation of model parameters used restricted maximum likelihood in the linear mixed models (lamina length, leaf shape index and adjusted lignotuber size), and residual pseudo-likelihood in the generalized linear mixed models (vegetative juvenility and seedling survival). The analyses of leaf shape index and adjusted lignotuber size involved the adjustment of the lamina width and lignotuber diameter measurements using lamina length and stem diameter as covariates, respectively; thus their least-squares means are adjusted estimates for seedling differences in lamina length and radial stem growth, respectively. For vegetative juvenility and seedling survival, the tabulated values in the "Least-squares means" columns pertain to the inverse linked scale (i.e., the observed data scale). For family variances, a one-sided likelihood ratio test (Self and Liang, 1987) was applied to test whether a given estimate was significantly greater than zero (p-values are given within parentheses). Lamina length (mm); leaf shape index = adjusted lamina width (mm); adjusted lignotuber size = adjusted lignotuber diameter (mm) (see the footnotes of Table 2 for further details). Standard error estimates for measuring the uncertainty of the estimated family variances are not presented. As the estimation of variance components usually constrains them to be positive, and their sampling distribution may thus be asymmetric, corresponding model-based estimates of standard errors may not provide accurate measures of uncertainty. The inaccuracy in standard error estimates is likely to occur for family variances, as the number of families sampled within each population was not large (Maas and Hox, 2004;McNeish and Stapleton, 2016b).
that of the dry population following acute drought stress. Further, we provide evidence suggesting an early age trade-off between drought susceptibility and a recovery trait, with seedlings which develop larger lignotubers being more susceptible to drought death. Such trade-offs may have contributed to the observed similar population survival, despite the marked population differentiation in seedling traits. Genetic variation within, and marked mean differences between, the mesic and dry populations for all four focal traits suggest adaptive differences between these populations, with seedlings of the dry population having on average smaller and broader leaves, longer vegetative juvenility, and larger lignotuber size than the mesic population. The mesic and dry populations compared were selected to accentuate rainfall differences and minimize the differences in other environmental factors. Accordingly, we expect that the observed population differences in the four functional traits reflect historic natural selection and adaptation to differences in water availability. Lamina length was strongly positively related to leaf area (r = 0.88, n = 380; unpublished data) and, as a surrogate of leaf size, its observed reduction in the dry population, for example, accords with the trait-environment associations reported in other Eucalyptus species (McLean et al., 2014). It is also consistent with species-level trends in south-eastern Australia, where species occurring in low rainfall areas tend to have smaller leaves through various changes in leaf length and width (McDonald et al., 2003). The longer vegetative juvenility observed in seedlings from the dry population also accords with the broader trends observed in a previous glasshouse trial with 37 E. pauciflora populations, where seedling vegetative 4 | Family variances estimated by residual-pseudo likelihood (RSPL) and maximum likelihood (ML), and estimates of intra-class correlation (ICC) coefficients, for discrete variables (vegetative juvenility and seedling survival) measured in the mesic (Brushy Lagoon) and dry (Ross) populations of E. pauciflora.

Family variance (RSPL)
Family variance (ML) Intra-class correlation For family variances estimated by ML, a one-sided likelihood ratio test (Self and Liang, 1987) was applied to test whether a given estimate was significantly greater than zero (p-values are given within parentheses). Vegetative juvenility = node number of the first petiolate leaf; seedling survival = survival measured as a binary outcome (where 1 = alive and 0 = dead). Estimates of family variances refer to the linked scale (i.e., logarithmic and probit scales for vegetative juvenility and seedling survival, respectively). Likelihood ratio tests were not pursued to test family variances under RSPL estimation, since RSPL uses linearization methods to approximate the model rather than the likelihood function, which precludes the comparison of nested models using likelihood ratio tests (e.g., Hedeker and Gibbons, 2006;McNeish, 2016). Residual variances are not provided: the generalized linear mixed models (GLMMs) fitted for vegetative juvenility and seedling survival do not contain a definite error term, and thus the residual variance is not defined as in the linear mixed model for continuous traits. The ICCs presented above correspond to GLMMs estimated with RSPL, and were calculated as described by Nakagawa et al. (2017) (see also Supplementary Methods 4). As in Table 3, standard error estimates for measuring the uncertainty of the estimated family variances are not given (see the footnotes in Table 3).  Non-parametric bootstrapping was pursued to provide standard errors, and 95% confidence intervals (within parentheses), for the selection gradient estimates and their population differences. Lamina length (mm); leaf shape index = adjusted lamina width (mm); adjusted lignotuber size = adjusted lignotuber diameter (mm); vegetative juvenility = node number of the first petiolate leaf. As traits were mean-standardized, a given selection gradient pertains to the effect on relative fitness-i.e., percentage change in relative seedling survival-expected from increasing the focal trait by 1% from its mean, while holding all the other traits constant. Statistical support against a null hypothesis being true (i.e., against the absence of an effect for either a selection gradient in a population or its difference between populations) is provided by a 95% confidence interval not overlapping with zero (indicated in Italics).
juvenility was retained longer in populations from drier homesites (Gauli et al., 2015;Costa e Silva et al., 2019). In this range-wide trial, the trait assessed was the percentage of nodes with petiolate leaves, which is negatively related to our measure of juvenility used in the present study. The difference between the dry and mesic populations is also consistent with trends in the E. risdonii/tenuiramis complex, where the retention of sessile juvenile leaves is extreme (into the reproductive stage) in populations from hot and dry environments (Costa e Silva et al., 2018). Populations of eucalypt species from hot/dry environments also have larger lignotubers than populations from mesic environments (Ladiges, 1974;Walters et al., 2005; Costa e Silva et al., 2019), as evident in our comparison. Despite mean differences between the mesic and dry populations in putatively adaptive traits, the two populations did not differ significantly in survival following acute drought stress as induced in our dry-down experiment ( Table 2). However, the two populations varied in the traits detected to be under directional selection, indicating that differences between the populations reside in the manner they selectively respond rather than in overall fitness under acute drought stress. Significant genetic-based family variance in survival was detected within the dry population, together with directional selection toward narrower leaves (smaller leaf shape index) and smaller lignotubers under acute drought stress. In contrast, no statistically significant family variance in survival was detected in the mesic population, and acute drought stress favored smaller leaf area (through shorter leaves) and earlier vegetative maturation (earlier development of petiolate leaves). However, statistical support for an important difference between the dry and mesic populations in directional selection gradients was provided for lamina length and leaf shape index only. For a given trait, the population difference in the directional selection gradients could in part be due to a difference between populations in the phenotypic variance (Steele et al., 2011;De Lisle and Svensson, 2017). This would particularly apply to lamina length, where the statistical support for a difference between populations in directional selection ( Table 5) could, to some extent, reflect a stronger selection in the more variable population (i.e., Brushy Lagoon; Table 3). Alternatively, the strength and direction of directional selection can be a function of the deviation of the population mean phenotype from the optimal phenotype for the environment (Steele et al., 2011;Hendry, 2017, p. 88). Regardless of the cause, given the significant genetic variation detected for all functional traits within each population (family variances in Tables 3, 4), evolutionary change would be expected in these traits following the observed directional selection acting on them (Lande, 1979;Lande and Arnold, 1983).
The directional selection observed in lamina length, reflecting leaf area, is consistent with the hypothesis that the initial trait difference between the two populations could reflect, at least in part, the effects of past natural selection on the trait to reduce seedling susceptibility to events of acute drought stress. However, with drought risk likely to increase in the wild with climate change (and thus with more frequent episodes of severe drought), the direction and strength of the selection gradient estimates in the dry-down experiment suggest that, in the future, directional selection would shift the mean lamina length of the mesic population toward that of the dry population. A reduction in leaf area may also underlie directional selection favoring narrower leaves within the dry population (i.e., smaller leaf shape index) in the dry-down experiment. As noted above, reducing leaf area is well-recognized as an adaptation to drier conditions. In dry environments, selection has been shown to favor smaller leaves (Dudley, 1996;Ramírez-Valiente and Cavender-Bares, 2017). Smaller leaves may facilitate adaptation to drier conditions via greater main vein density reducing hydraulic vulnerability (Scoffoni et al., 2011) or through having a thinner boundary layer of still air allowing better heat shedding than larger leaves, and thus reducing the need for transpiration and associated water loss to shed heat (McDonald et al., 2003). As also mentioned above, the difference between the two populations in directional selection for a trait could reflect adaptation to different trait optima. For example, selection for narrower leaves, but not shorter leaves, in the dry population may be due to the already shorter leaves of this population being close to a phenotypic optimum, and thus only the variation in leaf shape index would determine selection mediated by seedling survival. In Quercus, narrow leaves were favored in a dry year but not in a wet year, although a statistically significant year-bytrait interaction could not be detected (Ramírez-Valiente et al., 2015). However, there was evidence in the current study for an important difference between the mesic and dry populations in directional selection gradient estimates for leaf shape index. In a species-level study in south-eastern Australia, combinations of reductions in leaf width in addition to leaf length were responsible for a reduction in leaf area in low rainfall areas (McDonald et al., 2003).
The selection gradient estimates for lignotuber size showed an opposite trend to what would be expected if the differences between the mesic and dry population means were the result of differences in adaptation arising from past natural selection to reduce susceptibility to acute drought stress. The dry population had larger lignotubers than the mesic population, but experienced directional selection in favor of seedlings with smaller lignotubers, shifting its mean toward the more mesic population. The same trend in directional selection, toward smaller lignotubers, was evident in the mesic population, but without statistical support for an important effect on seedling survival (as evaluated based on seedling canopy mortality). Susceptibility to drought damage and recovery can be considered as genetically different "traits" in eucalypts (Ammitzboll et al., 2020). While our unexpected finding is consistent with lignotubers having a key role in drought recovery, it signals a cost to their development. Such cost suggests a resistance vs. tolerance trade-off, common in the evolution of adaptations to periodic strong selection events arising from biotic (Best et al., 2008;Moreira et al., 2014) and abiotic (Vesk and Westoby, 2003;Agrawal et al., 2004;Ramírez-Valiente and Cavender-Bares, 2017) stress. In the present case, the observed cost may be transient and only evident at this early seedling stage (Orians et al., 2010;Hoque and Avila-Sakar, 2015), or persist through the life-cycle but out-weighed by the benefits of the better development of a recovery organ and gains in tolerance.
Lignotubers are generally considered a store of non-structural carbohydrates and dormant buds to enable basal re-sprouting following the death of the main stem (Walters et al., 2005;Borzak et al., 2016;Ammitzboll et al., 2020). If the main stem is not severely damaged, it is also possible that reallocation of non-structural carbohydrates, such as starch, may help recover hydraulic function or photosynthetic activity (as suggested in Costa e Silva et al., 2019). However, the lignotuber does not appear to be an organ particularly specialized for the storage of non-structural carbohydrates such as starch, and likely serves primarily as a bud bank (Carrodus and Blake, 1970;Smith et al., 2018). The fact that selection favored smaller lignotuber sizes in the arid population under acute drought stress suggests that there is a cost associated with the early allocation of seedling resources to the development of this recovery organ. Non-structural carbohydrates support respiration and growth when photosynthesis is reduced during periods of water stress to avoid transpiration (McDowell et al., 2011;Mitchell et al., 2013). If the ready availability of nonstructural carbohydrates is reduced due to resources being allocated to the lignotuber, or their storage within the lignotuber, then seedlings with larger lignotubers may be more susceptible to "carbon starvation" during prolonged water stress, and thus death of the seedling canopy. Lignotubers themselves have also been reported to have smaller and contorted vessels compared to those in the stems, and their resistance to water flow may be twice that of the stem (Myers, 1995). An alternative cause of drought-induced selection favoring small lignotubers may thus relate to their offering less resistance to water flow between the roots and the stem. Regardless of the cause, with the arid population having significantly larger lignotubers than the mesic population, a trade-off between avoiding susceptibility to drought damage at the seedling stage vs. future seedling recovery could contribute to the absence of population differences in survival following the acute drought stress imposed. Indeed, with the evaluation of survival based on seedling canopy mortality at the termination of the experiment, re-watering and evaluation of recovery (Lu et al., 2010;Zhang et al., 2017;Ruehr et al., 2019) may well have provided the conditions under which population differences become manifest.
Drought-induced directional selection was found to favor early loss of vegetative juvenility in the mesic population, but had no effect on this trait in the dry population. There was thus no evidence from our phenotypic selection analysis that adaptation to acute drought stress would delay the development of the petiolate leaf in seedlings of the mesic population, and thus would shift the mean vegetative juvenility of the mesic population toward that of the dry population. Using growth performance as a fitness proxy, Costa e Silva et al. (2019) found no association between familylevel variation in seedling juvenility and the fitness proxy in a E. pauciflora common-garden field trial established on a dry site, which is consistent with the results currently obtained in the dry population. In contrast, a common-garden trial of populations from the E. risdonii/tenuiramis complex established on a wet site showed that directional selection favored phenotypes with early ontogenetic transition to the petiolate "adult" leaf (Costa e Silva et al., 2018). This finding is consistent with the selection gradient estimate for vegetative juvenility in the mesic population. However, environmental context is important to the interpretation of trait-fitness relationships (MacColl, 2011), and this result from a wet site does not explain why such selection should also occur under acute drought stress.
Translating results of phenotypic selection analysis into mechanistic responses can be challenging (Wade and Kalisz, 1990;MacColl, 2011;Kingsolver et al., 2012). Whilst all four focal traits showed evidence of directional selection associated with acute drought stress in one or other of the tested populations, it is possible that selection on a measured trait may occur via indirect selection acting on unmeasured traits, correlated with the measured focal trait (Lande and Arnold, 1983;Mitchell-Olds and Shaw, 1987). Indeed, of the four traits studied, vegetative juvenility would be the one most likely to reflect indirect selection involving unmeasured ("hidden") traits that could also influence seedling survival. While we assessed vegetative juvenility of the seedlings based on the node at which the first petiolate leaf was produced, this is potentially a surrogate for a wide range of physical, chemical and physiological changes that can occur, as documented in other eucalypt species (James and Bell, 2001;Borzak et al., 2015;Vlasveld et al., 2018;Lucani et al., 2019). It is therefore possible that the directional selection detected against vegetative juvenility in the mesic population could reflect indirect selection on a "hidden" (unmeasured) trait with a negative phenotypic correlation with vegetative juvenility. Nevertheless, it is also possible that there is direct phenotypic selection on the early development of the petiole, which is only manifest in the early maturing mesic population and could be specific to the timing of the onset of acute drought stress with respect to the seedling ontogeny. Such direct selection could, for example, relate to the enhanced potential for changing leaf orientation, a key trait affecting light interception and heat avoidance in eucalypts and facilitated by the petiole (Cameron, 1970;King, 1997;James and Bell, 2000).
Our selection analysis focused on survival of acute drought stress as evaluated based on seedling canopy mortality. However, historic adaptation to long-term dry environments may be expected to involve a complexity of different selection forces driving a suite of adaptive differences associated with water availability more generally, not only acute drought stress. From an evolutionary perspective, growth in water limited environments may require different adaptations to those involved in surviving periodic acute drought stress. Xylem embolism is believed to be a key cause of tree mortality under severe water stress (Urli et al., 2013;Brodribb, 2017;Choat et al., 2018), including eucalypts (Li et al., 2018), although long periods of water stress and a prolonged reduction in photosynthetic activity can also contribute to mortality through "carbon starvation" and susceptibility to biotic attack (Galiano et al., 2011;McDowell et al., 2011;Mitchell et al., 2013). These selective mechanisms are part of a broader suite of traits which control water uptake and loss (Li et al., 2018), the potential coordination of leaf and stem traits (Bourne et al., 2017) and the important role of recovery mechanisms in drought adaptation (Zeppel et al., 2015;Ruehr et al., 2019). The different roles of such factors in adaptation to acute vs. chronic water stress may have contributed to the absence of our population difference in the studied fitness component, despite differences in focal trait mean values. In particular, historic adaptation in the dry population may potentially involve selection favoring seedling traits, as well as adult traits, that enhance survival in water limited environments more so than survival during periodic acute drought stress.
Nevertheless, it is possible that population differences in fitness following acute drought stress do exist, but depend upon their pattern of acclimation (also referred to as preconditioning and priming- McDowell et al., 2008;Wang et al., 2017). Non-lethal water stress has been shown to result in marked metabolic changes in E. pauciflora seedlings with, for example, leaf carbohydrate concentrations increasing and contributing to osmotic adjustment (Warren et al., 2012). In other eucalypts, differential plastic pre-adjustments to water stress have already been described for functional traits at the species (Warren et al., 2011(Warren et al., , 2012, provenance (Li et al., 2000;McLean et al., 2014) and genotype (Pita and Pardos, 2001) levels in Eucalyptus. Further, it is possible that population differences in susceptibility to drought may only be expressed following a specific type of drought, such as one which includes water limitation combined with high temperature stress (Groom et al., 2004;Mitchell et al., 2014). This is particularly relevant as our dry-down experiment was conducted over autumn, and the maximum daily temperature the seedlings experienced averaged only 16-17 • C, and did not exceed 27 • C on any one day. The capacity to adapt to acute drought stress is likely to become more important into the future, with predictions of increasing extreme climate events. The results of this study suggest that: firstly, there may be genetic variation available in the focal traits (i.e., depending on the extent of genetic constraints in determining adaptive evolution; e.g., Costa e Silva et al., 2020) upon which they can respond to selection; and secondly, that the traits under phenotypic selection may vary between populations, reflecting in part differing historic selection pressures and adaptations, but also indicating potential different acute drought response strategies that may occur under future changes in the climate.

CONCLUSION
In conclusion, our dry-down experiment revealed different traits under directional selection between the studied dry and mesic populations, despite no difference in overall seedling survival. These results highlight the complexity of drought adaptation and the potential role of different traits and strategies in responses to acute drought stress. Further work on more populations from a wider range of environments, as well as replicated populations within the same environment, would help elucidate the various response strategies of E. pauciflora to drought stress, including key traits and mechanisms in different environments. Similarly, a better understanding of interacting factors, such as heat and drought pre-conditioning, in determining the drydown tolerance of plants is required. Ultimately, a greater understanding of the variation in traits under selection and potential mechanisms across environments will contribute to predictions of resilience, or vulnerability, of populations to future climate stress.

DATA AVAILABILITY STATEMENT
The data upon which this research work is based is currently being used in another study, but can be requested from the corresponding author for consideration.

AUTHOR CONTRIBUTIONS
JCS conceived and conceptualized the research question and analytical design, performed the data analyses, and prepared the associated document provided in Supplementary Material. Data for this study was collected from a dry-down experiment designed and undertaken by RJ, EP, and SMP. BMP and RJ were responsible for data scoring of the experiment and compiling the datasets used in the analyses. JCS, RJ, and BMP wrote the manuscript. All authors contributed to the manuscript revision, read, and approved the submitted version.

FUNDING
The contribution of JCS to this research work was supported by Fundação para a Ciência e a Tecnologia I.P. (FCT), Portugal, through the Norma Transitória DL 57/2016/CP1382/CT0008 and UID/AGR/00239/2019. This research was also funded by the Centro de Estudos Florestais, Portugal, a research unit that is funded by FCT (Unit Project Reference: UIDB/00239/2020). BMP acknowledges support of the Australian Research Discovery (ARC) grant DP190102053.

ACKNOWLEDGMENTS
We thank and acknowledge Dorothy Steane's contribution to the design and implementation of the dry-down experiment, and Tim Brodribb for helpful discussion. We also thank Paul Tilyard, Fiona Walsh, and Dale Worledge for assistance with running and scoring the dry-down experiment. We are grateful to all institutions mentioned for their financial support, which provided the opportunity to complete this study.