Genetic parameters of growth and adaptive traits in aspen (Populus tremuloides): Implications for tree breeding in a warming world

Aspen (Populus tremuloides Michx) is a widespread commercial forest tree of high economic importance in western Canada and has been subject to tree improvement efforts over the past two decades. Such improvement programs rely on accurate estimates of the genetic gain in growth traits and correlated response in adaptive traits that are important for forest health. Here, we estimated genetic parameters in 10 progeny trials containing >30,000 trees with pedigree structures based on a partial factorial mating design that includes 60 half-sibs, 100 full-sib families and 1,400 clonally replicated genotypes. Estimated narrow-sense and broad-sense heritabilities were low for height and diameter (~0.2), but moderate for the dates of budbreak and leaf senescence (~0.4). Furthermore, estimated genetic correlations between growth and phenology were moderate to strong with tall trees being associated with early budbreak (r = -0.3) and late leaf senescence (r = -0.7). Survival was not compromised, but was positively associated with early budbreak or late leaf senescence, indicating that utilizing the growing season was more important for survival and growth than avoiding early fall or late spring frosts. These result suggests that populations are adapted to colder climate conditions and lag behind environmental conditions to which they are optimally adapted due to substantial climate warming observed over the last several decades for the study area.


Introduction
Trembling aspen (Populus tremuloides Michx) is an ecologically and commercially important tree species with high genetic diversity and a broad natural range, including the boreal forest of North America, the eastern United States, and the western mountain ranges from Mexico to Alaska [1,2]. Aspen can regenerate both via sexual and asexual reproduction [2,3] Competing interests: I have read the journal's policy and the authors of this manuscript have the following competing interests: AH received a research grant that included matching financial contributions from industry partners to a government research grant for this study. JSB, representing the consultancy Isabella Point Forestry Ltd., received financial compensations from industry partners for his contributions to experimental design and analysis. The industry partners provided support in the form of research grants to AH and consulting contracts to JSB, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. All data collected by the authors on the timing of budbreak and leaf abscission will be made available through an online repository upon acceptance. However, our adherence to PLOS policies on sharing data and materials is altered with respect to height and diameter measurements that were not collected by the authors and that were made available to us by an industrial tree breeding cooperative. These data are not owned by the authors and can therefore not be made publicly available." initially tested a large number of clones collected from natural stands, targeting plus trees with good form and without signs of pathogens and diseases. In a similar manner, 122 individuals were selected as male or female parents for a partial factorial mating design (to be detailed below). The offspring were planted in ten progeny trials, established between 2005 and 2008 by the Western Boreal Aspen Corporation (WBAC), an industrial collaborative that includes Norbord Inc. (formerly Ainsworth Engineered Canada LP), Mercer Peace River Pulp (previously Daishowa-Marubeni International Ltd.), and Weyerhaeuser Canada Ltd. Detailed origins of the parental material for the progeny trials were previously documented [22]. There were 64 half-sib families as well as 100 full-sib families generated in a partial factorial mating design for the two breeding regions [22]. Each male and female parent was represented by one or two full-sib families and each female parent was also pollinated with a polymix to generate half-sib families. The pedigree structures of northern and southern breeding regions were constructed separately. Trials 1 and 2, planted in 2005, were seedling trials planted in the southern breeding region, sharing the same half-sib and full-sib families. Trials 3 to 8 were planted in 2007 and utilized families from both breeding regions. Families were clonally replicated prior to planting so that these trials have half-sib, full-sib, and clonal structure (i.e., multiple ramets of the same clone, with clones replicating multiple individuals of full-sib families). Trials 9 and 10 were established in 2008, containing different clonal material but overlap in half-sib and full-sib families with trials 3 to 7. These trials were connected through shared halfsib and full-sib families [22].
Seedlings for Trials 1 and 2 were grown in a greenhouse in April to May 2004, hardened in September to October, packed and cold-stored in a refrigerator in winter before planting in May 2005. Clonal planting stock for trials 3-10 was produced over a two-year cycle using a rootling methodology described by Brouard et al. [23]. First, seedlings were grown in 1-gallon pots (3.76 liters) at the Weyerhaeuser Tree Improvement Center in Drayton Valley, Alberta (53˚13'N, 114˚58'W, 869 m) to generate root mass. The following year, the root masses were washed and root cuttings were propagated in Beaver Plastics Styroblock-512 (60 cavities/block 220 ml plug volume) containers at Woodmere Forest Nursery in Fairview, Alberta (56˚04'N, 118˚24' W, 670 m). The resulting rootling ramets were then hardened, packed and cold-stored in a refrigerator in winter before planting in May 2007 and 2008 [23].

Experimental design and phenotypic measurements
All trials were constructed using an alpha design [24]. The use of alpha designs allowed for the flexible allocation of treatment and block numbers, and advantage over conventional randomized incomplete block designs. The location, exact experimental design, numbers of half-sib and full-sib families and clones for each trial are provided in Table 1. Over 30,000 individual trees were planted in this progeny trial series. Border trees surrounded each trial, and all trials except 7 and 10 were fenced to prevent browsing. Mulch layers or brush blanket mats were used to control competing vegetation, and spacing was 3×3m.
Height measurements were carried out multiple times between 2005 and 2013 with extendable measuring poles, and diameter at breast height (DBH) was assessed with a diameter tape. Budbreak scores were obtained based on a repeated scoring method according to Li et al. [25]. Score 0 was recorded for dormant buds, score 1 indicated a swollen bud, score 2 indicated broken bud scales, score 3 was given for the emergence of green leaves, score 4 indicated leaf extension, score 5 indicated more than two leaves emerged, and score 6 indicated fully unfolded leaves. Scores were recorded for each individual tree on April 12,20,22,24,26,May 1,3,9,11,13,15,17,19  For leaf senescence scores, scoring was based on an eight-level scale according to Fracheboud et al. [26]. Score 0 represented uniformly green leaves, score 1 indicated darker than pale green leaves, score 2 indicated a majority of pale green leaves, score 3 indicated more green than yellow leaves, score 4 indicated a majority of yellow leaves, score 5 indicated only yellow leaves, score 6 indicated 20% brown leaves, score 7 was 50% leaf abscission, and score 8 represented �90% leaf abscission. Observation dates for fall phenology in trial 2 were September 1, 7, 21, 29 October 5, 20 in 2011. Scores in trial 3 were recorded on September 1, 6, 23, 30 and October 6, 21 in 2011. Trial 8 was assessed on August 31, 7, 21, 29 September and October 20, 2011.
The day-of-year of a phenology event (DoY) was subsequently calculated for a critical score that showed the best separation among genotypes: score 3 for spring phenology (emergence of green leaves) and score 7 for fall phenology (50% leaf abscission). The date when the phenology reached the critical score of individual trees was determined by either the first record of the critical score, or by linear regression from the bracketing dates. The phenology data is provided as S1 Dataset in the supporting information section.

Statistical and quantitative genetic analysis
All quantitative genetic analysis was conducted with the ASReml-R package [27]. For the seedling trials (1 and 2) we employed the following mixed linear model: where Y is the vector of observations of traits (tree height, day-of-year of budbreak or leaf senescence etc.); β is a vector of fixed effects; a, f, r, b and p are vectors of additive genetic effects, full-sib family (cross) effects, replicate effects, block-within-replicate effects and plot effects, respectively; e is a vector of random residuals; X is the incidence matrix of the fixed effects relating β to observations Y; and Z 1 to Z 5 are the incidence design matrices relating the random effects a, f, r, b and p to observations Y. We assume that the observation vector follows a normal distribution with the expected value of E(Y) = Xβ and with the covariance matrix of Var(Y) = V, i.e., Y~N(Xβ,V). For our data, the β vector has only one element (the overall mean). The vectors of five random effects a, f, r, b, and p, as well as the vector of random residual e are assumed to follow the normal distributions Nð0; ; I e s 2 e Þ respectively. Here, s 2 A is the additive genetic variance, A is the pedigree kinship matrix for describing the additive genetic relationships among individual trees, s 2 f represents 25% of the dominance genetic variance, s 2 r ; s 2 b ; s 2 plot and s 2 e are the variance components corresponding to the vectors of random effects r, b, plot and residual e respectively, I t is the identity matrix of order t (t = f, r, b, plot, e). Thus, the total variance matrix can be partitioned into components due to the five vectors of random effects described above as well as the residuals, The best linear unbiased estimation (BLUE) of fixed effect (β) and best linear unbiased prediction (BLUP) of random effects (a, f, r, b, p) are solutions to the following mixed model equations,

½3�
For the trials with clonal single tree plot trials, we modified model (1) into the following linear mixed model: where the model remains the same as (1) except the plot effect is removed and the effects of clones within full-sib family (c) are added. The c factor accounts for the epistasis and ¾ of the dominance [28,29]. We also assumed that a, f, r, b, p, and e followed the normal distributions as above respectively; random effect c followed the normal distribution as N(0,I c σ C 2 ). Narrow-sense and broad-sense heritabilities were calculated based on following functions: whereŝ 2 A is the additive genetic variance component;ŝ 2 P is the phenotypic variance component represented by the sum ofŝ 2 A ,ŝ 2 NA andŝ 2 e ;ŝ 2 NA is the variance of non-additive genetic effects; the residual error isŝ 2 e [30]. The broad-sense heritability was estimated as: The standard errors of the heritability were calculated with the delta method [31]. We estimated the additive genetic correlation in seedling trials, genetic correlation in clonal trials (r G ), phenotypic correlation (r P ) based on individual trees observations. The linear model of tree growth and leaf phenology for a single site is: where y nimj is the observation of j-th tree of im-th genotype for the n-th trait; t n represents the n-th trait effect, r in is the i-th replicate effect of n-th trait, g nim is the additive genetic effect of m-th genotype of n-th trait in i-th replicate in the seedling trial, while in the clonal trails, g nim is the genetic effect. In the seedling trial, the genotypes are seedlings nested in replicates, while in the clonal trial, genotypes are clones evenly assigned in each replicate. In seedling trial 2, the genetic correlation is due to the additive genetic effect. In the clonal trials, the genetic correlation is due to the total genetic effect, though the additive genetic effect is more significant, and e nimj is the experimental error for each trial. The fixed effects of the mixed model are similar as the previous model, although the trait effect is added as a fixed factor. The genetic correlation (r G i j ) was of the form whereŝ G ij is the estimated genotypic covariance between trait i and j;ŝ G i is the estimated phenotypic standard deviation of trait i. For seedling trials the genetic correlation is the additive genetic effect. The phenotypic correlation was calculated as follows: whereŝ P ij is the estimated phenotypic covariance between trait i and j;ŝ P i is the estimated phenotypic standard deviation of trait i. Based on BLUPs, the breeding value reliability of half-sib parents and individual clones were calculated as follows: where R i is the reliability of the breeding value of the i-th parent, where PEV is the prediction error variance that equals to the standard error square of the predicted breeding value [32]; andŝ 2 A is the estimated additive genetic variance component. The correlations of breeding values among sites were calculated as the Pearson's correlation coefficients of half-sibs (breeding values) and clones (genetic values) between trials with chart. correlation of the R package PerformanceAnalytics. Bootstrapping of correlation coefficients of survival was carried out with the R boot package [33]. The G×E effect is explored with the Type-B genetic correlation for tree height, where the same trait measured in two or more environments over the same genetic composition can be treated as two different genetically correlated traits. R code for estimating genetic parameters according to the models above are available in the Appendix in [22].

Genetic parameters of growth traits
Genetic parameters for height and diameter were estimated at a relatively young age. Trees were between 5 and 8 years old at the time of evaluation with the oldest seedling trials having an average height of 3-4m and most clonal trials reaching average heights of 1-2m (Table 1). Dominance and epistatic variance components for height and diameter were small, less than 10% of the phenotypic variance component, so that most broad-sense heritabilities were only marginally higher than narrow-sense heritabilities ( Table 2). The highest heritabilities were estimated for a relatively small seedling trial 1 with values around 0.5. For all other trials, heritabilities for height and diameter were quite low (or unreliable) with narrow-sense heritabilities typically ranging from 0.1 to 0.2, and broad-sense heritabilities typically ranging from 0.2 to 0.3.
Type-B genetic correlations based on shared clones and shared full-sib families could only be calculated for sister trials (1-2, 3-4, 5-6, 7-8, and 9-10). Genetic correlations among sister trials yielded r GB values around 0.7 with a standard error of approximate 0.08 for the first three pairs, indicating a relatively low degree of genotype-by-environment interactions (G×E). Pairs 7-8 and 9-10 did not yield reliable estimates. Correlations of parental breeding values could be calculated for a larger number of trial pairs that shared parents through the partial factorial mating design. Parent breeding values between sister trials in trials 1 through 6 were generally well correlated, which can be interpreted as low G×E. Trial 8 showed a negative correlation with all other trials, and it should be noted that this trial was planted in a relatively cold environment with the highest elevation (Table 1).

Genetic parameters for adaptive traits
The phenology traits budbreak and leaf senescence were measured at three trials (Table 1), had moderate broad-and narrow-sense heritabilities (Table 3). Narrow-sense heritabilities for budbreak ranged from 0.4 to 0.5, while heritabilities for leaf senescence were slightly lower Table 2. Estimates of narrow-sense and broad-sense heritabilities at ten aspen progeny trials for tree height and diameter at breast height (DBH). Standard errors of the estimates are given in parentheses. with values between 0.3 and 0.4. Broad-sense heritabilities were only marginally higher than narrow-sense heratibilities, and could not be estimated for the seedling trial 2 ( Table 3). Also here, additive genetic variation was most important and dominance and epistatic variance components ranged from zero to 10% of the phenotypic variation. Moderate to strong genetic correlations were found between growth and phenology (r = -0.3 to -0.2 between height and bud break, and 0.6 to 0.8 between height and leaf senescence) with tall trees being associated with early budbreak and late leaf senescence (Table 4). Survival was not compromised by early budbreak or late leaf senescence. In fact, the reverse appeared to be true with negative correlations between survival and budbreak and positive correlations between survival and leaf senescence, just as for growth traits (Table 4).

Trait variation across the landscape
The correlations between fall phenology and height and survival were generally stronger than those between spring phenology and height and survival, indicating that the adaptive value of fall phenology is high. Neither fall phenology nor budbreak showed strong spatial patterns of breeding values across both breeding regions (Fig 2). However, correlations between height and leaf senescence are visible in these maps: comparing Figs 2 and 3, high breeding values for height (green dots) are often associated with late leaf senescence (pink and purple), with a Pearson's correlation coefficient of 0.65 (p<0.0001).

Discussion
Positive associations between height and survival as well as increased growth and survival in trees that break bud early and abscise leaves late suggest that utilization of the growing season may be more important than the avoidance of early fall frosts or late spring frosts at all three test sites where phenology was assessed. Strong additive genetic correlations between growth Table 3. Estimates of narrow-sense and broad-sense heritabilities at three aspen progeny trials for budbreak and leaf senescence. Broad-sense heritabilities were not estimated for the seedling Trial 02. Standard errors of the estimates are given in parentheses.  and phenology indicate that much of the genetic gain at the early stage of stand development will be due to expanding the growing season, which may increase the risk of frost damage in spring and fall.

Trial Code Age of Measurement
In studies with Populus species [34][35][36] and other species [37], high heritabilities of fall phenology traits, and positive genetic correlations between productivity and fall phenology have previously been observed. Genotypes with a delayed senescence in fall may nevertheless be exposed to increased frost risks. However, climate warming trends that have materialized over the last several decades in Alberta may have decreased the risks of early fall frosts [38]. This might potentially explain our positive association of an extended growing season in fall with high survival. Similar to the expectations of Olson et al. [39], genotypes that utilize a longer growing season may be favored by climate warming at northern latitudes. Inadvertent selection, where individuals with better height and diameter growth are chosen, and the genetically correlated leaf senescence is extended in fall, may therefore be an unplanned but effective climate change adaptation strategy. That said, other climatic factors are likely to change, and while the length of time suitable for growth is likely to continue to increase, water limitations may lead to an overall reduction of boreal forest productivity [40,41] In contrast, budbreak is a highly plastic trait in response to interannual variation and longterm trends in temperature. Populations can generally be expected to respond appropriately to climate change trends as long as daily temperature variances do not change for given baseline values. In other words, the frost risk associated with a certain heatsum that triggers budbreak needs to remain the same. Yet, the day-of-year for budbreak may have shifted, when this heatsum is reached. Late spring frosts are also considered a more severe threat than early fall frosts, as they may destroy buds and juvenile leaves and severely compromise early-season growth [42]. In our study, genetic correlations between growth and budbreak were negative (i.e. early budbreak associated with better growth), but they were not nearly as strong as genetic correlations with leaf senescence. Furthermore, survival was not compromised in genotypes that started the growing season relatively early. We, therefore, conclude that correlated selection for early budbreak poses only a small risk. In fact, the results could be interpreted as indicating an adaptational lag with respect to changing climate. The utilization of a longer growing season appears to increase the growth and survival at least in this sample of three field tests. Our within population genetic analysis corroborates results from provenance research that demonstrated adaptational lag among populations through wide-ranging reciprocal transplant experiments [43].
The previous interspecific hybrid breeding mainly exploits heterosis/specific combining ability for growth and wood quality trait improvement [9,15,44]. In our study, intra-specific crosses and within family selection after field selection provide potential gain for growth and predictable phenology response due to the non-additive genetic effect such as specific combining ability and epistasis effect. Instead of producing hybrids with northern species such as P. davidiana and P. tremula [9], for regional tree improvement, adjacent breeding zones could exchange parentages, crosses, and clones in places like Alberta without severe frost risk concern in spring and fall even in the northern plantation sites [10].
We noted that heritabilities found in this study were generally much lower than those reported in previous studies by Gylander et al. [10] in comparable trial series that investigated wild clonal selections, and Gylander et al. [10] reported broad-sense heritabilities of 0.51-0.58 from clonal trials. Kanaga et al. [45] Calculated broad-sense heritabilities for growth traits ranging from 0.30 to 0.50 in a short-term common garden study using 13 aspen clones. Our low heritability estimates for growth traits are likely due to strong environmental microsite variation at post-harvest planting sites and the juvenile age at which trees were evaluated. Heritabilities were particularly low for Trials 7 and 10, and only Trial 7 could not be fenced and showed evidence of browsing, while Trial 8 was located at the more elevated site in the study area. For growth selection at multiple sites, results of >8 growing seasons are adequate for mediocre sites and the productive site. We conclude that selection for growth traits at the current stage promises only small to moderate genetic gains, but higher heritabilities may emerge at a later date of trial evaluation.
Genetic correlations among fall phenology and spring phenology traits with growth traits were nevertheless already high in this study, and could further increase as the influence of microsite variation of the planting site decreases with the age of the trees. Survival was also not compromised, but was positively associated with early budbreak or late leaf senescence, indicating that the growing season length was more important for survival and growth than avoiding early fall or late spring frosts. We found that a substantial portion of the tested genotypes are adapted to a shorter growing season than they have experienced during the testing period. Selecting genotypes for reforestation that utilize a longer growing season may be unproblematic under continued climate warming in northern latitudes, and may lead to overall increases in forest productivity as long as water availability under increased evapotranspiration does not become a dominant limiting factor.