A new approach to prediction of the age-age correlation for use in tree breeding

Abstract

Key message

Early selection in tree breeding requires a credible age-age correlation. Modelling height growth in provenance and progeny trials, we can predict age-age correlations suitable for use in operational breeding as described in this article.

Context

Tree breeding involves early selection, which is an indirect selection using a genetic correlation. This study describes a procedure of predicting an age-age phenotypic correlation as a surrogate for a genetic correlation. Although the predicted correlations are based on white spruce (Picea glauca) and lodgepole pine (Pinus contorta) data, they can be used in other coniferous species with similar mode of height growths.

Aims

The aim of the study is to predict a correlation coefficient used to adjust breeding values at a measurement age to breeding values at a rotation age. This correlation is derived from the observed height growth trajectories of trees in progeny and provenance trials.

Methods

Correlation prediction equations were developed using modelled height growth in provenance and progeny trials of lodgepole pine and white spruce. The time lag between successive tree ages was used as a correlation predictor variable.

Results

Correlations differed between spruce and pine but the differences narrowed as trees grew older. For example, a correlation between 20 and 100 years was 0.607 for spruce and 0.470 for pine, whereas that of 30 and 100 was 0.826 for spruce and 0.832 for pine. Based on the age-age correlation, the optimum selection age for a 100-year rotation age is 40–50 years. Parameters of the tree height growth function exhibited significant genetic variance and genotype × environment interaction.

Conclusion

After the age of 40 years, age-age correlation for height may be less important for selection and genetic gain prediction than the correlation between height and diameter, which is declining with tree age.

1 Introduction

In commercial forestry, optimum rotation age (ORA) for wood production is the age at which a forest plantation yields the maximum profit (Chang 1984). Genetics, climate, soil properties, silviculture, and other factors that affect tree growth determine ORA. When wood production (yield) is the goal of forest management, tree breeders select genotypes (parent trees or clones) that will maximize yield at ORA. The question is how to identify such genotypes at an early age? How to estimate genetic gain at an early age without overestimating genetic gain at ORA?

While efforts are underway to use DNA marker-aided and genomic selection (e.g., Grattapaglia and Resende 2011; Resende et al. 2012; Isik 2014), field progeny trials remain the primary tool for estimating genetic parameters, predicting breeding values, selecting genotypes and predicting genetic gain. When trees are young, ranks of genotypes for growth traits (height, diameter, and volume) change over time (e.g., Mullin and Park 1994; de Sousa et al. 2005). A correlation among measurements of the same trait at different ages is called an age-age correlation. This correlation is high when genotypic ranks are stable over time and low when they substantially fluctuate. A low age-age correlation suggests that the genetic gain predicted at an earlier age will overestimate genetic gain at ORA, because not all genotypes selected for their superior growth at a young age would have been selected at ORA.

Tree breeders frequently use the term optimum selection age (OSA) to mean the age beyond which changes in the age-age correlation are minor. Consequently, selected genotypes and percentage genetic gain predicted at OSA and ORA ought to be similar. Because OSA is usually unknown, Zobel and Talbert (1984) recommended selecting genotypes when the age of trees in the progeny trials is at least half the ORA. At this age and size, trees are old enough for breeders to be confident that ranking of genotypes will remain relatively stable. White et al. (2007) lists other selection criteria linked to ORA. It suffices to say that these selection criteria and rules of thumb are not always feasible in all species. For example, mid-ORA is certainly feasible in some tropical and subtropical species where ORA is 20–30 years such that 8–10 years of field testing is adequate (e.g., Gill 1987; Cotterill and Dean 1988). In contrast, in the interior northern boreal conifers in Canada, mid-ORA of 40–60 years is certainly too long to delay selection and prediction of genetic gain. Therefore, breeding boreal conifers requires a different approach that allows for selection at much younger ages while avoiding overestimating the genetic gain at ORA.

The need to undertake early selection without overestimating expected genetic gain at ORA is of particular importance in Canada where the public owns 94 % of forested land (Natural Resources Canada 2014). Private forest companies manage these forests through forest management agreements with governments (Beckley 1989). When a company plant trees with a specified expected genetic gain, it receives an equivalent increase in allowable cut from the existing forests. This increase in today’s allowable cut in exchange for an expected yield increase in future forests is called allowable cut effect (Luckert and Haley 1995). Although the allowable cut effect (ACE) provides an immediate return on investment for companies, governments bear the risk by offering genetic gain that might not be realized at ORA. Therefore, for governments, the age-age correlation has both technical tree breeding and public policy implications. For example, to mitigate the risk, Alberta and British Columbia, Canada, mandate the use of the age-age correlation to convert genetic gain at the measurement age to genetic gain at ORA.

The use of the age-age correlation to predict genetic gain at ORA from the genetic gain predicted at an earlier measurement age is consistent with a concept of correlated response to selection (Falconer and Mackay 1996). The challenge is how to obtain a correlation that meaningfully and convincingly relates observed values of a trait at a measurement age and expected values at ORA. To be used in tree breeding, such a correlation must be estimated in a way that takes into consideration the biological nature of the way trees grow and be based on meaningful predictors. Currently, Alberta and British Columbia use the correlation from the equation developed by Lambeth (1980). There are legitimate concerns about the operational use of correlations from this equation, which are addressed in this article.

In this article, I review methods that have been used to obtain age-age correlations for use in tree breeding and present a new method whereby correlations are linked to height growth trajectories in provenance and progeny trials. The present work used data from lodgepole pine (Pinus contorta Dougl.) and white spruce (Picea glauca [Moench] Voss) trials in Alberta. The two species make up more than 80 % of reforestation in Alberta. Nevertheless, the method and correlations developed in this study can be used in other coniferous species with similar mode of height growth. For simplicity, correlations for ORA of at least 50 years are included in the tables. Correlations for short ORA can be obtained by substituting an appropriate predictor in the presented equations.

1.1 Theory and current practices

The use of an age-age correlation (r _t,T ) in tree breeding is based on the quantitative genetics concept of correlated response to selection (Falconer and MacKay 1996). If we treat height at a younger (H _t ) and older (H _T ) age as different traits, we can predict how selection for H _t will change H _T as a correlated response (Eq. 1).

$$ C{R}_{H_T}={r}_{t,T}i{h}_t{h}_T{\sigma}_{H_T} $$

(1)

where $C{R}_{H_T}$ = correlated response (change) in H _T due to selection for H _t ; i = selection intensity at age t; ${\sigma}_{H_T}$ = phenotypic standard deviation for H _T ; h _t and h _T = square root of the heritability for H _t and H _T , respectively; r _t,T = genetic correlation between H _t and H _T . Without knowing the variance and heritability for H _T , which will be observed far in the future, it is prudent and for practical reasons to assume that h _t = h _T ; i will be the same whether selection is done on H _t or H _T , and the phenotypic standard deviation for H _t approximates ${\sigma}_{H_T}$. With these assumptions, Eq. 1 simplifies to,

$$ C{R}_{H_T}={r}_{t,T}i{\sigma}_{H_t}{h}_t^2 $$

(2)

where $i{\sigma}_{H_t}{h}_t^2$ = expected genetic gain for height at a measurement age (H _t ). Therefore, expected genetic gain at ORA can be estimated by a simple multiplication of expected genetic gain at a measurement age by an age-age correlation (r _t,T ). This is why r _t,T is an important statistic for tree breeders working with long ORA species. To be used operationally, r _t,T must be as realistic as possible. Although this article deals with r _t,T for height growth, the same concepts apply to other traits.

Some form of a correlation between H _t and H _T should be expected because of the cumulative nature of perennial height growth (Eq. 3).

$$ {H}_T={H}_t+{H}_i $$

(3)

where H _i is the growth increment accrued since the last time H _t was measured.

Thus, r _t,T is a correlation between H _t and H _t  + H _i (Eq. 4).

$$ {r}_{t,T}=\frac{\operatorname{cov}\left({H}_t,{H}_t\right)+\operatorname{cov}\left({H}_t,{H}_i\right)}{\sqrt{\operatorname{var}\left({H}_t\right)\operatorname{var}\left({H}_T\right)}}=\frac{\operatorname{var}\left({H}_t\right)+\operatorname{cov}\left({H}_t,{H}_i\right)}{\sqrt{\operatorname{var}\left({H}_t\right)\operatorname{var}\left({H}_T\right)}} $$

(4)

Consequently, even if H _t and H _i were not correlated, H _t will be correlated with H _T to the extent that H _t is a component of H _T . The amount of H _i added annually will determine the rate at which H _t ceases to be a significant component of H _T and the rate at which r _t,T declines with tree age or size. The greater the H _i relative to H _t the lower the r _t,T and vice versa. In the juvenile phase when trees have high annual height growth increments (AHI), H _i will accrue faster with an increase in tree age than in a mature phase when AHI is low. Thus, attempts to predict age-age correlations must recognize that changes in AHI in the life of trees will affect r _t,T and OSA.

There are plenty of examples of observed r _t,T in the forestry literature (e.g., Xie and Ying 1996; Hodge and White 1992; Tauer and McNew 1985; Kung 1973). Attempts have been made to use these correlations to develop equations for predicting ${\widehat{r}}_{t,T}$ beyond the observed field trial periods, which are often much shorter than ORA. Using reported r _t,T , Lambeth (1980) developed ${\widehat{r}}_{t,T}=1.02+0.308\times LAR$, where LAR = ln (t/T). Since then, many equations involving LAR in various forms have been fitted to height, diameter, and volume (e.g., Lambeth and Dill 2001; Gwaze et al. 1997; Jansson et al. 2003; Ye and Jayawickrama 2012). The LAR-based models are by far the most attempted and cited ${\widehat{r}}_{t,T}$ prediction equations. As mentioned earlier, the Lambeth (1980) equation is used for operational breeding in Alberta and British Columbia, Canada (Xie and Yanchuk 2003). It is also widely used in forest biometrics when the age-age correlation is involved to incorporate genetic gain into yield models (e.g., Newton 2015).

The prevalence of LAR-based models is likely due to the simplicity of getting t/T that fits every situation. As a ratio, the range of t/T is the same (0.0–1.0) regardless of the age of trees, length, and number of serial measurements. In principle, we can fit a LAR equation based on t/T from any material and use it to obtain ${\widehat{r}}_{t,T}$ that is similar to those from Lambeth (1980) and many other LAR-based models. For example, Fig. 1 shows a regression of r _t,T on LAR for white spruce seedlings. These seedlings were raised in the greenhouse and their heights measured every 2 weeks for 36 weeks distributed equally over two growing seasons (Rweyongeza et al. 2004). With height growth from germination to the end of the first growing season (18 weeks), the equation is ${\widehat{r}}_{t,T}=1.0686+0.3047\times LAR$ with r ² of 0.919. At the end of the second growing season (36 weeks), the equation is ${\widehat{r}}_{t,T}=0.9888+0.2804\times LAR$ with r ² of 0.876 (Fig. 1). Both equations are very similar to Lambeth (1980) and other reported LAR-based models. An age ratio can be substituted into these equations to obtain a correlation for practical use in a real breeding program.

Fig. 1

Age-age correlation in white spruce seedlings measured at 2-week intervals for 36 weeks (two growing seasons) in the greenhouse. Because of the age ratio as a predictor variable, the prediction equation could be used to predict correlations for trees of any age

In LAR-based models, t/T is just scaling the predictor variable to make it appear as if the observed r _t,T used to fit the equation were observed over the entire rotation age (ORA). Hence, substituting t/T in a LAR-based equation is simply obtaining correlations for intermediate ages within the range of the correlation matrix of observed r _t,T with which the equation was fitted. Therefore, ${\widehat{r}}_{t,T}$ from LAR-based models are not predictions of future ${\widehat{r}}_{t,T}$ as intended and are therefore misleading.

Contending that perennial height growth is a cumulative trait, which is a function of annual increments (AHI), and r _t,T is determined by $LAG=T-t$ among AHI, Kremer (1992) fitted ${\widehat{r}}_{t,T}=1.079-0.132\times LAG-0.0039\times (LAG)$ ². He then did simulations to see how ${\widehat{r}}_{t,T}$from a LAG-based model is affected by the (i) input correlation matrix of AHI, (ii) change in the additive coefficient of genetic variation for AHI, and (iii) modelled AHI to 50 years. He concluded that the structure of the input correlation matrix of AHI and the changes in the genetic variance of AHI (${\sigma}_{AHI}^2$) were the major determinants of variation in ${\widehat{r}}_{t,T}$. He further observed that, while r _t,T declined with an increase in LAG as expected, the change in ${\sigma}_{AHI}^2$ did not follow any age pattern. Hence, the randomness of changes in ${\sigma}_{AHI}^2$ questions the value of individual AHI in predicting ${\widehat{r}}_{t,T}$. In a separate study of the same species (Pinus pinaster Ait), the change of r _t,T with LAG for AHI was completely random (Costa and Durel 1996).

There is a large body of literature linking the rate of height and diameter growth to existing tree size (e.g., Bond et al. 2007; Vanderklein et al. 2007; Niklas 2007). Therefore, AHI formed in a specific year is intrinsically dependent on the total height that existed prior to its formation rather than individual AHIs formed in previous years. Kremer (1992) found that r _t,T for AHIs separated by 13 years was almost zero. This, together with the randomness of ${\sigma}_{AHI}^2$ suggests that r _t,T of AHIs may not be useful in predicting ${\widehat{r}}_{t,T}$. Individually, AHI reflects the variation in weather and other temporal environmental factors that affect tree growth. On the long term, these temporal variations in tree growth are averaged out making total height a better predictor of ${\widehat{r}}_{t,T}$ than AHI.

Kung (1993) argued that, due to its symmetry, the correlation matrix can be viewed as a symmetrical response surface model with a ridge on the diagonal (r _t,T = 1.0) and slopes inclining toward the two corners with r _t,T decreasing with an increase in LAG. After fitting ${\widehat{r}}_{t,T}={b}_0+{b}_1t+{b}_2T+{b}_3{t}^2+{b}_4{T}^2+{b}_5tT$ (where b ₁ = b ₂ and b ₃ = b ₄) and finding that first order and interaction terms were not statistically significant, the model was reduced to ${\widehat{r}}_{t,T}={\beta}_0+{\beta}_1\times {LAG}^2$ (Model1). In addition, Kung (1993) argued that, in the same way r _t,T depended on LAG, the degree of non-determination (DON = 1 – r ², where r = r _t,T ) would depend on LAG. Hence, he fitted DON = β ₀ + β ₁ × LAG (Model2) and ${\widehat{r}}_{t,T}=\sqrt{\left({\beta}_0+{\beta}_1\times LAG\right)}$ (Model3). Comparing ${\widehat{r}}_{t,T}$ estimated from Model1, Model2, Model3, and Lambeth (1980) against observed r _t,T , he concluded that Model2 and Model3 overestimated, whereas Model1 and Lambeth (1980) underestimated the correlation.

Gwaze et al. (1997) fitted both the LAR- and LAG-based equations and compared ${\widehat{r}}_{t,T}$ from both equations with corresponding values from the Lambeth (1980) equation. It was observed that the LAG equations fitted the data better than LAR equation and produced ${\widehat{r}}_{t,T}$ that were consistent with observed values than did the Lambeth (1980) equation. It suffices to say that, unlike LAR, the LAG-based models do not have the problem of scaling the predictor variable to erroneously imply that the observed r _t,T values used to fit the model spanned the entire range of ORA.

An obvious problem one would encounter in a LAG-based model is its inability to predict ${\widehat{r}}_{t,T}$ beyond the range of the data used to fit the model, except when fitting a simple linear and second degree polynomial equation. My initial attempts showed that fitting any other model that introduces a curvilinear relationship between r _t,T and LAG produces an equation, which when used to predict ${\widehat{r}}_{t,T}$ for much larger LAG yields values that are very close to ${\widehat{r}}_{t,T}$ corresponding to the largest LAG of the input correlation matrix. This points to the inability to extrapolate this model beyond the range of LAG in the observed data. The correlation prediction equations developed in the present work have addressed this problem.

The other observed feature of the Lambeth (1980) equation is that all t/T yielding the same ratio have the same predicted correlation (${\widehat{r}}_{t,T}$). For example, when t/T is 0.5, ${\widehat{r}}_{t,T}$ is 0.81 even though this may be a correlation between ages 2 and 4, 25 and 50, 50 and 100 years or any other mid-ORA selection ages. This is an unlikely expectation arising from the way Lambeth (1980) developed his equation. The same feature would be encountered in all LAR- and LAG-based equations reviewed in this article. This feature is a result of fitting a single equation through an entire correlation matrix.

Plotting r _t,T from any source against LAG will show that a correlation matrix is a collection of many scatter plots, each corresponding to correlations between a measurement at a specified age (X _t , X _t + 1, X _t + 2, X _t + 3, …, X _t + n ) with measurements at subsequent ages (Y _t + 1, Y _t + 2, Y _t + 3, …, Y _t + n ). If we fit an equation for each of these scatter plots, the same LAG will have different predicted correlation (${\widehat{r}}_{t,T}$) depending on the equation from which it was predicted. This is illustrated in Fig. 2 using data from Rweyongeza et al. (2004) for correlations involving seedling heights at weeks 2, 4, 6, and 8 and heights up to 36 weeks. At the same value of LAG, older trees should have higher ${\widehat{r}}_{t,T}$ than young trees. Fitting one equation for the entire age-age correlation matrix is equivalent to averaging all ${\widehat{r}}_{t,T}$ from individual age-specific equations thereby assigning the same ${\widehat{r}}_{t,T}$ to all cases where LAG or LAR is the same regardless of the age of trees. The correlation prediction equations developed in the present work have addressed this problem.

Fig. 2

Age-age correlation for height growth in white spruce seedlings grown in the greenhouse for 36 weeks illustrating the importance of fitting age-specific correlation prediction equations. The number in brackets indicate the age of seedlings and the equations correspond to regression lines from bottom to top, respectively

2 Materials and methods

2.1 Data description

The data used in this study came from a series of white spruce and lodgepole-jack pine “complex” provenance and progeny trials scattered across Alberta. The term “complex” is used here to imply that white spruce trials may contain hybrids of P. glauca (Moench) Voss and Picea engelmanii (Parry ex Engelm). Likewise, trials of lodgepole pine (P. contorta Doug var. latifolia [Engelm.]) may contain jack pine (Pinus banksina Lamb) and lodgepole-jack pine hybrids. The hybridization between white and Engelmann spruces and lodgepole and jack pines occurs naturally in Alberta. Details of the data are summarized in Table 1.

Table 1 Description of Alberta provenance and progeny trials used in the study

2.2 A new method of predicting age-age correlations

The method of predicting ${\widehat{r}}_{t,T}$ presented here (i) avoids using LAR because this variable does not predict realistic ${\widehat{r}}_{t,T}$ as intended; (ii) recognizes that trees to do not grow indefinitely at the same rate (AHI); (iii) uses LAG as a predictor variable while avoiding a simple linear and polynomial regressions; and (iv) enables LAG to predict ${\widehat{r}}_{t,T}$ meaningfully beyond the range of the data from which the equation was developed.

Studies show that forest trees exhibit growth phases with different AHI. The growth rate is high and exponential during the juvenile phase. As the juvenile phase ends, trees attain vegetative and morphological complexity and reproduction begins, AHI declines (Kramer and Kozlowski 1979). Tree height growth follows a sigmoid growth function whereby AHI is lower, higher and lower in the early, middle and mature phase, respectively (Kramer and Kozlowski 1979). Conifers spend hundreds to thousands of years in the mature phase (Kramer and Kozlowski 1979) with few centimeters of AHI while expanding in diameter. The culmination of AHI and OSA lies in the lower portion of the mature phase. The time trend in the genetic variance for height growth appears to follow these height growth phases (Namkoong et al. 1972; Namkoong and Conkle 1976). Therefore, a realistic tree height growth model must permit for a declining growth rate in the mature phase. Only with such a model will age-age correlation attain an optimum value at OSA.

A simple linear equation is to be avoided because it implies that ${\widehat{r}}_{t,T}$ declines with an increase in LAG at the same rate throughout ORA. This is inconsistent with the sigmoid height growth pattern observed in perennial plants. A quadratic equation implies that after attaining OSA, ${\widehat{r}}_{t,T}$ would decline with an increase in tree age (t). Kremer (1992) advanced this idea by assuming that during ORA, the correlation would increase with t in the first 1/3 phase because fast growing genotypes have high AHI than slow growing ones. In the second 1/3 phase, the correlation would remain constant because AHI of fast growing genotypes has peaked. In the last 1/3 phase, the correlation would decline reversing the first phase trend because AHI of slow growing genotypes has surpassed that of fast growing genotypes. It is possible that slow growing genotypes may have a relatively higher AHI toward the end of ORA than fast growing ones. However, this does not translate into rank reversals for total height, which is the basis for genotypic selection. Thus, the quadratic and other polynomial equations are not attempted in the present article.

The new method of predicting age-age correlations involved (a) predicting height at different ages between age 5 and 120 years, (b) generating the age-age correlation matrix from the predicted heights, and (c) developing age-age correlation prediction equations using this correlation matrix.

A review of the Alberta forest inventory data showed that in the boreal forest, the average height of 120 years old white spruce and lodgepole-jack pine complex is 20 and 23 m. The data in Meng and Huang (2010) and Lotan and Critchfield (1990) support this generalization for lodgepole pine. Table 1 shows that the average height for trees in the oldest white spruce trials were close to half the size expected at ORA in the natural stands. Likewise, height of lodgepole-jack pine complex trials between age 25 and 30 years could be expected to be half the height expected at the ORA in the natural stands.

The initial attempt was to fit a sigmoid growth curve to observed data as follows,

$$ {H}_t=\frac{k}{1+b{e}^{-rt}}+\varepsilon $$

(5)

where H _t = total height at age t (years); k = upper asymptotic height; b = a constant with no biological interpretation (Richards 1959); e = the base of the natural logarithm; r = the growth rate; and ε = residual. The time at the point of inflection (t _0.5) occurs mid-way between the upper and lower asymptote (${\scriptscriptstyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}k$) and is calculated as ${t}_{0.5}={\scriptscriptstyle \frac{1}{r}} \ln (b)$. This logistic equation is further described by Nair (1954). Review of the results showed that the growth rate would begin to decline early such that the predicted height at ORA would be lower than the actual height observed in natural stands. This is because the latest H _t measurement greatly influenced the point of inflection (Meng and Huang 2010). Therefore, Eq. 6 was used to predict height at 3–5-year intervals beyond the latest H _t measurement until predicted heights were in the range of 20–25 m similar to heights expected at ORA in natural stands.

$$ {H}_{kjin}=a{X}^b+{\varepsilon}_{kjin} $$

(6)

where H _kjin = total height at age t (years) of nth tree in ith family or provenance in jth replication (block) at kth test site; ε _kjin = the residual; a and b are regression coefficients.

Therefore, these additional predicted height points were combined with observed heights to create “hybrid datasets” for fitting a sigmoid growth curve on individual tree basis (Eq. 5) using PROC NLIN (SAS Institute 2004).

Pearson’s correlation coefficients between observed height and height predicted by Eq. 5 were greater than 0.95 at all sites. Individual-tree height growth functions (Eq. 5) were used to predict total height (Ĥ _t , that is, height predicted by age) from age 5 to 120 years. Pearson’s correlation coefficients for Ĥ _t were calculated on individual species and single-site basis using PROC CORR (SAS Institute 2004). Earlier analysis showed that Pearson’s correlation coefficients between observed total height and tree age in calendar years were greater than 0.95. Therefore, the age difference LAG in calendar years among successive Ĥ _t was used as a predictor variable for ${\widehat{r}}_{t,T}$ (Eq. 7).

$$ {\widehat{r}}_{t,T}={\beta}_0{e}^{\beta_1{d}_1} $$

(7)

where d ₁ = T − t = $LAG$ in calendar years; β ₀ and β ₁ = regression coefficients; e = the base of the natural logarithm; and all other terms are as previously defined. To maintain consistence with previous notations, LAG is used in place of d ₁ for the rest of this section.

To allow for ${\widehat{r}}_{t,T}$ to differ among measurements with the same LAG, separate equations were fitted for each t instead of fitting a single equation for the entire correlation matrix. For example, the equation for predicting ${\widehat{r}}_{t,T}$ involving height at age 5 years (Ĥ ₅) was developed using correlations (r _5,T ) involving Ĥ ₅ with heights at subsequent ages (Ĥ _T ), where T = 10, 15, 20, …, 120. The equation for predicting ${\widehat{r}}_{t,T}$ involving height at age 6 years (Ĥ ₆) was developed using correlations (r _6,T ) involving Ĥ ₆ with heights at subsequent ages (Ĥ _T ), where T = 10, 15, 20, …, 120. The equation for predicting ${\widehat{r}}_{t,T}$ involving height at age 10 years (Ĥ ₁₀) was developed using correlations (r _10,T ) involving Ĥ ₁₀ with heights at subsequent ages (Ĥ _T ), where T = 15, 20, 25,…, 120. The equation for predicting ${\widehat{r}}_{t,T}$ involving height at age 11 years (Ĥ ₁₁) was developed using correlations (r _11,T ) involving Ĥ ₁₁ with heights at subsequent ages (Ĥ _T ), where T = 15, 20, 25…, 120. Equations for predicting ${\widehat{r}}_{t,T}$ of all other ages were developed using the same sequence as those illustrated above. These equations were fitted using PROC NLIN (SAS Institute 2004).

In addition to the age-age correlation prediction equations, the intraclass correlation was calculated to measure variation among families and provenances for parameters of the logistic growth functions k, b, r _, and t _0.5. Where progeny trials included provenances from bulk seedlots, they were dropped before performing the analysis of variance. All analyses of variances were implemented in PROC MIXED (SAS Institute 2004) as described below (Eqs. 8 and 9).

$$ {y}_{ijn}=\mu +{\alpha}_i+{\beta}_j+\alpha {}_i\beta_j+{\varepsilon}_{ijn} $$

(8)

where y _ijn = observed value of the nth tree in jth provenance (or family) in the ith replication; μ = site mean; α _i = effect of the ith replication; β _j = effect of the jth provenance (or family); α _i β _j = provenance (or family) × replication interaction (experimental error); and ε _ijn = residual. Except μ, all effects were considered random effects with σ ² _α , σ ² _β , σ ² _αβ , and σ ² _ε variance components, respectively.

$$ {y}_{lijn}=\mu +{\tau}_l+{\alpha}_i\left({\tau}_l\right)+{\beta}_j+{\tau}_l{\beta}_j+{\varepsilon}_{ljin} $$

(9)

where y _lijn = observed value of the nth tree in the jth provenance (or family) in the ith replication within the lth test site; μ = general mean; τ _l = effect of the lth test site; α _i (τ _l ) = effect of the ith replication within the lth test site; β _j = effect of the jth provenance (or family); τ _l β _j = provenance (or family) × site interaction; and ε _lijn = residual. Except μ and τ _l , all effects were considered random effects with σ ² _α , σ ² _β , σ ² _τβ , and σ ² _ε variance components, respectively. Intraclass correlations were calculated as in Eq. 10 (for provenance or family on individual sites) and Eq. 11 (for provenance or family across sites) using respective variance components on individual sites and across sites.

$$ {g}_i=\frac{\sigma_{\beta}^2}{\sigma_{\beta}^2+{\sigma}_{\alpha \beta}^2+{\sigma}_{\varepsilon}^2} $$

(10)

$$ g{}_{ac}=\frac{\sigma_{\beta}^2}{\sigma_{\beta}^2+{\sigma}_{\tau \beta}^2+{\sigma}_{\varepsilon}^2} $$

(11)

where g _i and g _ac = intraclass correlation in individual sites and across sites, respectively.

3 Results

In this study, the power function (Eq. 6) was fitted to observed data only to provide few additional data points for fitting the logistic growth function (Eq. 5). The pseudo r ² for Eq. 6 functions were greater than 0.95 showing a near-perfect fit to the model. Detailed results from this stage of the analyses are not presented in this article.

Summary statistics for parameters of the growth functions appear in Table 2. Also included in Table 2 are intraclass correlations (g _i ) on individual sites at provenance and family levels, which measure the genetic variability for parameters of the growth function among populations and families. It can be seen that the level of genetic variation (provenance or family) for parameters of the logistic growth functions differed considerably among sites within the same series of trials. Table 3 summarizes variance components as percentages of the total variance and the intraclass correlations across sites (g _ac ). Lower values of g _ac (Table 3) compared to g _i values (Table 2) are indicative of a substantial genotype × environment (GE) interaction in the parameters of the logistic growth function at both the provenance and family level. Based on the Wald Z statistic, the GE interaction was statistically significant (P < 0.05), which is evident in the relative value of σ ² _β and σ ² _τβ when expressed as percentages of the total variance (Table 3).

Table 2 Summary statistics and family or provenance intraclass correlations for parameters of the logistic growth functions for spruce and pines in Alberta

Table 3 Variance components as percentages of the total variance for cross site analyses of parameters of the logistic growth function (Eq. 5)

Equations for predicting ${\widehat{r}}_{t,T}$ are summarized in Table 4 for t of 10 to 50 years. For white spruce, the correlation between observed and predicted values were 0.30–0.88 (pseudo r ² = 0.09–0.77). Corresponding values for lodgepole pine were 0.28–0.86 (pseudo r ² = 0.08–0.74). Low values are due to the fact that, in developing ${\widehat{r}}_{t,T}$ prediction equations, the input correlations (r _t,T ) were not averaged across sites. This variation in r _t,T across sites lowers the correlation between observed and predicted values. In contrast, if r _t,T is averaged across sites prior to developing ${\widehat{r}}_{t,T}$ prediction equations, the correlation between observed and predicted values is 0.84–0.94 (pseudo r ² = 0.71–0.88) for white spruce and 0.63–0.97 (pseudo r ² = 0.40–0.94) for lodgepole pine. Whether or not input ${\widehat{r}}_{t,T}$ values are averaged across sites prior to fitting, ${\widehat{r}}_{t,T}$ prediction equations does not change the resulting ${\widehat{r}}_{t,T}$ prediction equations. Otherwise, all equations in Table 4 were statistically significant (P < 0.0001). Table 5 contains ${\widehat{r}}_{t,T}$ values for selected long ORA that would normally be encountered in northern temperate and boreal countries such as Canada. Correlations for other ORA can be obtained by substituting LAG in respective equations (Table 4).

Table 4 Age-age correlation prediction equations

Table 5 Predicted age-age correlations for long rotation ages (ORA) in white spruce and lodgepole pine

4 Discussion

This study showed that families and provenances varied significantly for the parameters of the logistic growth function and this variation was greatest for k (Table 2). Because k is the prediction for height the trees can potentially attain, its variability reflects the general extent of genetic variation in height growth as previously reported (Rweyongeza et al. 2007, 2010). Variation for r and t _0.5 were greater in lodgepole pine than in white spruce (Table 2). Lodgepole pine is more shade intolerant (Lotan and Critchfield 1990) than white spruce. This study used height growth data from trials that have closed canopy. If competition in closed canopy trials affects the trajectory of tree height growth among families or provenances, its effects on g _i and g _ac would likely be more visible in lodgepole pine (shade intolerant) than white spruce (shade tolerant). This underlines the need for having species-specific or genus-specific age-age prediction equations (where data exist) instead of a single equation for all conifers as in Lambeth (1980).

The GE interaction in the parameters of the growth function is equally the result of GE interaction in height growth as previously reported for white spruce (Rweyongeza 2011) and Pinus pinaster Ait (Danjon 1994). The GE interaction for parameters of the growth functions will not affect the application of the age-age correlation developed in this study (Table 4). This is because, (i) in principle, these correlations represent cross-site averages of ${\widehat{r}}_{t,T}$ that could be generated on individual sites in the same way as the cross-site breeding values do, and (ii) the correlations will be applied beyond the narrow environment of the test site. Therefore, these correlations take into consideration the variability in the response of families and provenances to the environment encountered in actual reforestation programs.

The consequence of fitting separate correlation prediction equations for t instead of a single equation for the entire correlation matrix is clearly demonstrated in this study. For example, for white spruce (Table 4), ${\widehat{r}}_{t,T}$ is 0.859 (ages 25 and 50); 0.896 (ages 30 and 60); 0.963 (ages 40 and 80); and 0.990 (ages 50 and 100 years). For lodgepole pine (Table 4), ${\widehat{r}}_{t,T}$ is 0.811 (ages 25 and 50); 0.889 (ages 30 and 60); 0.980 (ages 40 and 80); and 0.998 (ages 50 and 100 years). Under Lambeth (1980), these mid-ORA selection ages would have the same correlation of 0.81.

This study used phenotypic correlations generated from height measurements of individual trees to develop age-age correlation prediction equations. Lambeth and Dill (2001) and a review by White et al. (2007) suggest that the phenotypic correlations predicted by the Lambeth (1980) equation are always lower than corresponding genetic correlations. White et al. (2007) observed that the reported genetic correlations are usually 0.05 to 0.2 greater than the phenotypic correlations from Lambeth (1980). However, as explained earlier in this article, the way the Lambeth (1980) equation was developed precludes any realistic comparison with observed genetic correlation from any trial and species. In addition, published genetic correlations are often associated with high standard errors and some exceed the permissible range of −1.0 to 1.0 (e.g., Lambeth et al. 1983; Tauer and McNew 1985). This makes unbiased Pearson’s correlations the more appropriate substitutes to guard against using erroneously high correlations to develop age-age correlation prediction equations for practical use. According to Namkoong and Kang (1990), the presence of a phenotypic correlation does not guarantee presence of a genetic correlation. However, for an age-age correlation that is partly due to autocorrelation (Eq. 4), a high level of similarity between a genetic and phenotypic correlation should be expected. Hence, for practical purposes, phenotypic correlations developed in the present study are considered good substitutes for genetic correlations.

Examination of the correlations in Table 4 shows that at ORA of 100 years, OSA likely lies between 40 and 50 years. Coincidentally, this is close to half-ORA Zobel and Talbert (1984) advocated. Indeed, between 35 and 50 years, ${\widehat{r}}_{t,T}$ changes only slightly. This suggests that delaying selection or even continuing to monitor the height ${\widehat{r}}_{t,T}$ after age 35 years is unnecessary. Height, diameter at breast height (DBH), and taper determine volume, which is a measure of wood production in forestry. Figure 3 shows a decline in the correlation (r _h,d ) between height and DBH as trials grow older. The range of r _h,d in Fig. 3 is 0.76–0.94 (mean = 0.85) for white spruce and 0.58–0.95 (mean = 0.76) for lodgepole pine. The decline is obviously greater in lodgepole pine than in white spruce.

Fig. 3

Age-related decline for the correlation between height and diameter at breast height

According to Lotan and Critchfield (1990), stand density greatly affects diameter and yield per hectare in lodgepole pine. Stand density is much lower in field trials than in fire-origin natural stands. Nevertheless, any shading due to crown closure in field trials will affect lodgepole pine more than white spruce. In Scots pine (Pinus sylvestris L.), Kroon et al. (2008) showed that genetic and phenotypic correlations between height and diameter on three sites were less than 0.80. Volume production was genetically and phenotypically better correlated with diameter (r > 0.95) than height (r < 0.85). Huang et al. (1992) showed that the relationship between height and diameter of all major forest tree species in Alberta was nonlinear. Therefore, evidence does not support estimation of volume genetic gain based on selection for height growth only. One way to overcome this problem is to estimate genetic gain directly from volume when trees are old enough to provide meaningful DBH measurements. For young trees where height measurements are the only reliable data, prediction of volume genetic gain based merely on height breeding values carries some risk of selecting wrong genotypes and overestimating genetic gain.

5 Conclusions

This study developed a method of obtaining age-age correlations for converting height genetic gain at a measurement age to genetic gain at a rotation age. The method is based on modelling height growth trajectories of white spruce and lodgepole pine in provenance and progeny trials in Alberta. The age-age correlation prediction equations developed in this study were based on height predicted by the logistic growth function. Parameters (b, k, and r) of the growth function have impact on ${\widehat{r}}_{t,T}$, because they determine the predicted heights. According to Richards (1959), b depends on the timing of the first measurements, whereas r determines the shape of the growth curve. As the asymptotic parameter, k represents height expected to be reached by a tree during the prediction period, whereas t _0.5 is the point at mid k. Because all parameter are estimated from tree measurement data, the number of measurements, timing, and intervals between serial measurements will definitely affect the logistic growth functions, the predicted heights, and consequently ${\widehat{r}}_{t,T}$. It is expected that measuring trees long enough at short regular intervals beginning at early ages will provide the best prediction of tree height growth trajectories and consequently more reliable age-age prediction equations. The earliest and latest height measurements for the data used in the present study are 5 and 32 years, respectively. Although measurement intervals differ among trial series (Table 1), they are short enough to provided adequate tracking of the tree growth trajectory. This attest to the strength and reliability of the age-age correlation prediction equations developed in the present study.

Based on ${\widehat{r}}_{t,T}$ and the assumptions behind the methodology, it can be concluded that selection at 40–50 years is sufficient for estimating height genetic gain at 100 years with no need for adjustment with the age-age correlation. Although the correlations and correlation prediction equations presented in this article were developed using data from white spruce and lodgepole pine in Alberta, they may be used for other conifers with similar mode of height growth. For example, white spruce correlations may be used in other spruces, whereas lodgepole pine correlations may be used in other pine species. Moreover, the assumptions on the conifer height growth model employed in this study are considered realistic enough to allow the correlations to be used in conifers other than Pinus and Picea species. Use of these correlations for deciduous species is left for tree breeder’s discretion, because the mode of height growth of coniferous and deciduous species may be very different.