Incorporating non-key traits in selecting the Pinus radiata production population

Three of the traits considered of most economic importance in the genetic improvement of Pinus radiata D. Don, termed as ‘key’ traits, are tree diameter (a proxy for stem volume), wood density and wood stiffness. A number of other traits (non-key traits) may be of equal importance to growers depending on where and for what purpose the trees are being grown. Two such non-key traits are: resistance to the needle blight caused by Dothistroma septosporum (Dorog.) M. Morelet, and reduced heartwood content. Data from two trial series (each planted at two sites, from a total of some 330 families (189 were half-sib and 142 were full-sib) with 10 to 30 individuals assessed per family) were analysed to determine the effect forwards selection of key traits had upon genetic gains of these two non-key traits. Multivariate analysis for each trial provided estimates of trait narrow-sense heritabilities (h2) and genetic correlations between traits. Density was the most heritable trait assessed (ĥ2 0.50 to 0.72) with heartwood ring number (0.21 to 0.35), acoustic velocity (0.40), resistance to Dothistroma septosporum (0.20 to 0.34) and stem diameter (0.11 to 0.25) being less heritable. Wood properties were adversely correlated with growth rate to varying degrees. To estimate the impact of differing technical weights on multiple traits a selection index model was used. Strong positive genetic correlations (r g 0.69 to 0.87) between resistance to Dothistroma septosporum and stem diameter means that strong selection for stem diameter after severe Dothistroma attack assures genetic gain in resistance to Dothistroma septosporum. Strong selection for stem diameter compromised wood properties due to adverse correlations between the two. Heartwood ring number was almost uncorrelated with the other key traits, density and stiffness, meaning that zero or slight negative gain would be expected via indirect selection. In such instances, it is advocated that key traits be selected for in the breeding population using a selection index and that non-key traits such as heartwood be selected using independent culling in the production population.


Background
Tree breeding primarily involves selection of superior genotypes from breeding populations that are expected to improve a number of economic traits in the production populations. Selection indices are commonly used in tree breeding to improve multiple traits simultaneously (Cotterill and Dean 1990). To understand the potential trade-offs made during selection for a suite of traits, it is important to understand how such traits are correlated. Selection of measurement traits are constrained by both their cost of measurement, and by their ability to affect the profitability of breeding objectives that are based on defined products and/or industrial processes (e.g., for Pinus radiata D. Don; see Apiolaza andGarrick 2001, Ivković et al. 2006). The three key traits (KTs) for P. radiata now considered to be of most economic importance are diameter (DBH, a proxy for stem volume), density (DEN) and stiffness (PME) (Pers. Com. John Butcher). However, depending on the location of the stand or the intended end-product, one or several other traits referred to as non-key traits (NKTs) may also be of importance to the grower. The expected impacts of forwards selection for these three KTs on two NKTs, resistance to Dothistroma septosporum (Dorog.) M. Morelet (DTR) and heartwood ring number (HWRN), are examined in this paper.
Dothistroma needle blight causes foliage loss and, in moist warm climates, results in growth loss of up to 70 percent (Shaw and Toes 1977) making it an important selection trait for plantings in regions at high risk of infection. Heritability estimates for DTR have been extensively reported. For example, Ivcović et al. (2010) reported narrow-sense heritability estimates ranging from 0.05 to 0.69, with a median of 0.36 in a series of 16 P. radiata trials assessed in Australia for DTR. A positive genetic correlation between DTR and growth at a later age was also reported, with a median value of 0.39, however no correlation was found between density and DTR (Ivcović et al. 2010).
Heartwood in P. radiata is undesirable due to its darker appearance, unreliable durability and difficulty to treat chemically. These properties combine to devalue wood containing heartwood, making increased sapwood content of interest to wood processors. A prior paper, involving the same genetic material, investigated the genetic relationship of heartwood/sapwood with diameter growth (Kennedy et al. 2013) finding sapwood traits to be moderately heritableĥ 2 i ¼ 0:2 to0:4 and sapwood cross-sectional area to be strongly correlated (r g = 0.9) with diameter growth. Heartwood content is strongly correlated phenotypically with diameter growth in conifers (Climent et al. 2002), but evidence of some independent genetic control has been found in P. radiata (Nyakuengama et al. 2000) particularly if assessing HWRN as used in this paper (Kennedy et al. 2013). The practicality of simultaneously selecting for DTR and HWRN, and the three KTs, is dependent largely on the heritability of each trait, and the genetic correlations between them. It may eventuate that some NKTs are favourably correlated with some of the KTs, and in this situation the correlated response of the NKT to selection on the KTs means that the NKT may also benefit indirectly. Conversely, if adverse correlations between the KTs and NKTs exist, simultaneous selection for both categories will prove more difficult.
A multi-trait selection index model was produced to illustrate the generalised benefits and constraints of particular selection strategies, to identify situations where alternative selection approaches might need to be considered and to test a number of selection scenarios among KTs and NKTs. The implications of selecting for KTs and the NKTs of HWRN and DTR on the gain that might be achievable in the production population was examined using individual-tree breeding values generated from mixed-model analysis.

Methods
The field trials Data from two Pinus radiata field-trial series were analysed in this study, the 'Female Tester' series and the 'Dothistroma' series.
The Female Tester series was planted at two sites in New Zealand, one at Esk Forest, Hawkes Bay (39°15' S, 176°42' E) and the other at Woodhill Forest, Northland, (36°46' S, 174°21' E). The trials were of single-tree-plot, sets-within-reps design consisting of 30 replicates. Each trial was made up of five sets, containing 36 seedlots, of which five in each set were controls with Growth and Form (GF, see Vincent and Buck 1998) ratings varying from 6 to 27. The other seedlots were Female Tester crosses; such crosses were composed of a number of pollen parents that had been crossed factorially with five female parents. This was adopted to better compare the top-ranking parents.
Two Dothistroma needle blight screening trials were planted at Kinleith Forest, Central North Island, New Zealand. The trials were planted in 1985 at Phoenix Road (38°23' S, 176°13' E) and Bobcat Road (38°24' S, 175°99' E). The design of these trials was a single-treeplot, sets-within-reps design consisting of 30 replicates. The seedlots were divided into seven sets, each consisting of around 25 seedlots, with all seedlots in set A (not included in the analysis) being self-pollinated. The remainder of the sets were made up of open-pollinated families. Between three and eight controls were included in each set.

Sampling and assessment
Trait abbreviations and units are detailed in Table 1. Summary details of the number of families and individuals assessed for each trait at each trial are given in Table 2. The level of DTR was assessed on each of Dothistroma series trials at ages two, three and six years from planting. Resistance was scored visually as percent of total crown depth not infected by Dothistroma septosporum, in 5% classes. Wood properties are more difficult and costly to assess resulting in only a sub-sample of individuals from each family being assessed for DEN, HWRN and PME ( Table 2). The number of heartwood rings and density were assessed on 5-mm pith-to-bark cores extracted at breast height from trees aged 16 at Esk, 17 at Woodhill, and 23 at Bobcat and Phoenix. The heartwood ring number of each core was used in the trial analysis to select trees that delayed the initiation of heartwood development (i.e. had fewest heartwood rings) to maximise sapwood area. Basic density was assessed, using the maximum moisture content method (Smith 1954) on methanol extracted cores, across rings 1 to 5 from the pith at Esk and rings 6 to 10 from the pith at the two Dothistroma series trials. Standing-tree acoustic velocity was converted to PME using Young's modulus equation (Wang et al. 2001): where ρ is the green density of the material (kg m −3 ) and V is the velocity of the wave through the material (m s −1 ).
The green density of a standing tree was assumed to be a constant, and so a direct relationship between V 2 and PME was assumed for the analyses conducted here.

Statistical analysis
Series were analysed separately, due to limited overlap of genetic entries across the two series. Controls within each trial series were dropped from datasets before analysis. Traits at each site were analysed using an individual-tree linear mixed-effects model using ASReml software (Gilmour et al. 2009). The following base models were used for each series: where y is the vector of observations on a trait, d is a vector of fixed effects (i.e., female, and mean), r is a vector of random replicate effects, a is a vector of random additive genetic effects of individual genotypes, f is a vector of random specific combining effects of families and e is a vector of random residual effects. Set effects were insignificant and not included. Parameters X, Z 1 , Z 2 and Z 3 are known incidence matrices relating the observations in y to effects in d, r, a and f respectively.
where y is the vector of observations on a trait, d is a vector of fixed effects (i.e. mean), r is a vector of random replicate effects, a is a vector of random additive genetic effects of individual genotypes (assuming open-pollinated families to be half-sibs), and e is a vector of random residual effects. Set effects were insignificant and not included. Parameters X, Z 1, and Z 2 are known incidence matrices relating the observations in y to effects in d, r and a respectively.
Base models 1 and 2 were adjusted to include site effects, site was fitted as a fixed effect and all model terms fitted within site. Type B genetic correlation estimates (r B ) between genotypes at different environments for each trait (Burdon 1977) were obtained directly using the CORR command within the ASReml program to structure the G matrix for the additive genetic component estimated by each model; all other model terms were considered independent for each site. Type B genetic correlations for the two Dothistroma series trials were sufficiently strong (~0.9, Table 3) to warrant an across-sites analysis by dropping the site term from the vector of random additive genetic effects of individual genotypes.
Variance components and breeding values used in the multi-trait selection model were estimated for the KTs and NKTs of interest, at each site or series, using where σ 2 a is the additive genetic variance of individual genotypes, σ 2 f the specific combining ability variance of families and σ 2 e the residual effects. As the Dothistroma series trials were open-pollinated it was not possible to estimate σ 2 f for trials within this series. Standard errors of statistics were based on approximations using Taylor series expansion (Lynch 1998) within ASReml.
The additive genetic coefficient of variation (CV a ) was estimated to compare the relative genetic variances of traits across trials: where σ a is the square root of the additive genetic variance and x is the trait mean. And additive genetic correlations between traits were estimated as: where σ a x a y is the covariance between the traits, σ 2 a x is the additive genetic variance for trait x andσ 2 a y is the additive genetic variance for trait y. The estimated additive genetic correlation r a xy between traits was obtained directly using the CORR command to structure the G matrix for each model.

Multi-trait selection model
An index model was constructed in Microsoft Excel software based on the traditional 'Smith-Hazel' index model (Cotterill and Dean 1990) to explore the dynamics of selecting for the KTs and NKTs when different technical weights (level of importance) were imposed. An outline of the model, using notation adapted from Burdon (1989), is presented in the equations below: A selection index (I) can be constructed to account for multiple traits as follows: I ¼ b 1 X 1 þ b 2 X 2 :::::::: where b 1 X 1 is the product of the weight b given to the breeding value X for trait 1, and so forth.
The solution for b that theoretically maximises desired genetic gain is where b is the column vector of index weights 1, 2…n; P is the phenotypic variance-covariance matrix; G is the genotypic variance-covariance matrix; b is the column vector of technical weights.
The expected gain from such an index in terms of a column vector EΔ is given by the formula: Where i is the selection intensity (e.g. 2.42 = a selection rate of 2%).
The selection intensity used in the model was fixed at 2.42 (2%), leaving the only variable as the weight imposed on each trait under selection. The weights applied to the model were altered for each trait to achieve the desired, or close to the desired, gain for the suite of traits under selection (i.e., they are considered technical weights). A number of selection scenarios were investigated by altering the set of technical weights used, in order to illustrate the trade-offs in gain between the traits (Figures 1, 2 and 3). All gains were calculated as a percentage change relative to the population mean for the model in use.
The individual-tree breeding values estimated from the linear mixed-effects multi-trait analysis of the Dothistroma series trials were used to determine the effects of independent culling upon the NKTs after index selection of the KTs. A selection index was constructed for the two KTs assessed (DBH and DEN). The index took the form of Eq. 6, with appropriate weights within the multi-trait selection model applied to DBH and DEN to provide gain estimates in line with breeding objectives 1 for this hypothetical round of the breeding cycle (10% increase in DBH and a 10 kg m −3 increase in DEN; note: this approximately reflected a weight of 2 on DBH and 1 on DEN). Non-key traits were omitted from the selection model, as they would not be considered at this stage of selection. Technical weights were multiplied by the individual-tree breeding values, as calculated by Eq. 1. The top 2% of individuals (252 from 12600) for the index were selected. The 10 best individuals for HWRN were then selected from this top 2%.

Genetic parameters
Type B genetic correlations between the Bobcat and Phoenix sites were strong for all traits assessed (r g~0 .90, Table 3). The estimated heritabilities for traits between these sites were very similar (Table 3) and consistent with the heritability estimated from the across-sites analysis (Table 4). This is unsurprising, as the two trials were planted within very close proximity to one another and established in the same year. Genetic parameters estimated from the multivariate analysis of the trials are presented in Table 4: DEN was negatively correlated with DBH (r g~− 0.4) at the three sites where it was assessed. However, there would appear to be little correlation between DBH and PME (r g range −0.19 to 0.27) with the sign of the correlation varying by site; DTR was strongly correlated between the ages of two, three and six Figure 1 Expected genetic gain for four traits at Esk for three selection scenarios. A key to the weights used in each selection scenario for the respective traits is provided; a selection rate of 2% was used in the model.

Figure 2
Expected genetic gain for three traits at Woodhill for three selection scenarios. A key to the weights used in each selection scenario for the respective traits is provided; a selection rate of 2% was used in the model.  Trait means along with their estimated additive genetic coefficients of variation (CV a %) are also presented along with estimated Type B genetic correlations (r B ) between Esk and Woodhill for traits in common.
years (Table 4), therefore for the remainder of the paper DTR03 will form the basic unit of measurement for this trait. Diameter at breast height was the least heritable trait (ĥ 2 0.11 to 0.25) across all four sites and density the most heritable trait (ĥ 2 0.5 to 0.73) across the three sites it was measured.

Multi-trait selection
For a given set of weights, the predicted gain for each trait calculated using Eq. 8 was plotted, to help build up the picture of the trade-offs when employing multi-trait selection. The graphs in Figures 1, 2 and 3 illustrate the dynamic process that gain is allocated to each trait under multi-trait selection. Across all sites, it is apparent that HWRN cannot be reduced (i.e. a negative gain achieved) without a large reduction in optimal DBH gain. Heartwood ring number is almost uncorrelated to any of the KTs, making it difficult to achieve large gains in both. The pattern of relative gain estimated for the two trial series were similar, but gains were lower for DBH in the Female Tester series trials largely due to the relatively low heritabilities at both the Esk and Woodhill sites (i.e. <0.16). The parameter DTR03 is favourably correlated with DBH06; the bigger the tree the lower the incidence of DTR infection. Owing to this favourable genetic correlation, selection for growth alone will result in an increase in DTR resistance. Strong selection for DBH will, however, compromise wood properties such as DEN and HWRN ( Figure 3). As HWRN is not strongly correlated to the two other KTs, DEN and PME, selection for the three KTs is unlikely to decrease HWRN.

Selection indices and independent culling
The basic approach for selection of NKTs used here was to restrict consideration of these traits to the production population using independent culling to provide targeted seedlots that meet individual growers' needs. An additional advantage of this approach would be that only trees among the prospective production population need be assessed for NKTs of interest, making considerable cost savings. An example of how this approach might work follows.
The top 2% of individuals (252 from 12600) for the index (for the two KTs DBH and DEN) were selected from the Dothistroma test population to form a prospective production population (Figure 4). Selection of these 252 individuals provided an estimated gain relative to the population mean of 21.6% for DBH (24 mm) and 1.6% for DEN (5 kg m −3 ). These gains differed from the predicted gains for DBH and DEN from the general Smith-Hazel index model (18.8% and 6.7% respectively), possibly due to the dependency on precise estimation of variancecovariance matrices to accurately predict breeding values. The predicted breeding values as used in the general index model do not take into account the differing accuracies of the estimation of individual breeding values which may explain the difference in gain (Schneeberger et al. 1992). As such, the gain comparisons among the various scenarios using the technical weights from the general Smith-Hazel index (which is operating as a mass selection index) do not directly apply when applying the same weights to the tree breeding values from mixed-model analysis. Regardless, the relative changes and tradeoffs in gain among traits can be adequately probed with the Smith-Hazel index model, and attaining certain levels of gain from direct best linear unbiased prediction (BLUP) estimates will require additional but more cumbersome sensitivity analysis.
Based on selection of the 252 individuals for DBH and DEN, the associated gain for the NKTs of HWRN and DTR was −6.5 and 53.4% respectively. For growers Figure 4 Bivariate distribution of individual-tree breeding value estimates (% departure from trial mean) for DBH and DEN for the two Dothistroma series trials. Filled diamonds are those individuals selected at intensity of 2% (252 out of 12600 individuals). Index weights were estimated from the multi-trait selection model. The relative gains for all 4 traits of the selected 252 are presented in the top right-hand corner.
where a decrease in the number of heartwood rings is of importance, this slight gain may not be sufficient. However, specific seedlots from the production population could be targeted to minimise the number of heartwood rings using an independent-culling approach. For example, selecting the 10 trees with the lowest HWRN out of the 252 selected for DBH and DEN resulted in a 19.7% decrease in HWRN ( Figure 5). However, these 10 trees come from just three families, and a reduction in selection intensity for independent culling would be required to reduce relatedness in this very small selected population.
This approach of using independent culling to improve NKTs not strongly correlated with any of the KTs is only likely to lead to gains in the current production population, as gene frequencies that improve the NKTs are not increased in the breeding population. The trade-offs in gain selecting for just a few uncorrelated traits, demonstrated earlier, illustrate the importance of concentrating improvement upon only a few traits of most economic importance. As P. radiata is grown on a large number of site types and has such a wide variety of uses, it has been a challenge to keep the breeding objectives down to only a few traits (Shelbourne 1997). Calculation of economic weights for NKTs is likely to prove difficult, as a particular trait may be of importance for some growers and producers and for others not, for example, with widely varying incidence of Dothistroma infection across different regions. Earlier generations of the New Zealand breeding programme partitioned the main breeding population into several breeds that concentrated on improving specific traits (Jayawickrama and Carson 2000). These breeds were later combined to form the current breeding population which had the benefit of reducing costs and increasing selection intensity. Identification of other NKTs uncorrelated to the KTs used for selection may see the return of specific breeds if enough demand for their improvement exists.

Conclusions
If reliable estimates of genetic correlations between the three identified KTs and a NKT of particular interest can be obtained, then it would be possible to increase the weight given to the NKT in question, elevating potential gain from selection. The main drawback to such an approach would lie in the situation whereby the NKT was unfavourably correlated with a number of the KTs. Considerable penalty would be incurred in the unfavourably correlated KTs if any gain was to be achieved through selection for the NKT(s) in question. Under normal circumstances, if such a situation was to arise, one would select for this NKT based on independent culling, but absence of measurement data means this cannot generally be considered at present. More likely, the trees selected for the three KTs would be assessed for the NKTs of interest, reducing operational measurement costs. Non-key traits are best selected from within the production population rather than the breeding population.
Sensitivity analyses using the Smith-Hazel selection index model enabled a general understanding of the levels of gain and trade-offs among key and non-key traits, and the effects that choosing different populations can have on gain. While the desired technical weights, derived from the traditional mass-selection index approach, cannot be directly applied to multivariate BLUP breeding values (and yield the same levels of genetic gains) they do provide an important 'starting point' for further sensitivity evaluations. Both approaches can provide insights and good quantification of what may be coming out of production populations when independent culling is applied to NKT's. As the results presented here indicate, positive gains can still be made in a trait that is not correlated with the key traits under selection. Figure 5 Bivariate-distribution of individual-tree breeding value estimates (% departure from trial mean) for HWRN and DTR for the two Dothistroma series trials. Filled diamonds are the 252 individuals selected on an index for DBH and DEN. Triangles represent the 10 trees of the selected 252 with the lowest HWRN. The relative gain for all four traits of the re-selected 10 is presented in the top right-hand corner.