Genetic Parameters of Diameter Growth Dynamics in Norway Spruce Clones

: The breeding of Norway spruce in northern Europe has substantially contributed to the production of high-quality wood. The vegetative propagation of robust elite clones could help to sustain the provision of high-quality timber in the face of changing climates. For the adequate evaluation of genetic gains, the altered tree growth dynamics of the clones need to be understood, yet essential information about the long-term growth dynamics of improvedboreal trees is still lacking. We examined a 50-year-old clonal plantation in Latvia to distinguish the clonal effects on diameter growth function parameters and estimate the genetic parameters. A mixed-effect modelling approach was used, in which the clones were applied as random effects on the parameters of the Chapman–Richard equation. All model parameters showed signiﬁcant variance in the genotypic coefﬁcients of variation CV g which ranged between 11.0 and 17.1%, with the highest being for the growth rate. The heritability ( H 2 ) of the diameter at breast height (DBH) reached 0.35 at the age of 40, while CV g decreased from 12.9% to 7.8% between the ages of 20 and 45. Age–age genotypic correlations were positive and were strong or very strong (>0.76). The realised genetic gain varied from − 6.3 to +24.0% around the trial mean. A substantial improvement in DBH was indicated when elite clones were selected for vegetative propagation based not only on early measurements, but also considering the genetic variance in the model parameters.


Introduction
In the Baltic Sea region, the planting of forests has a substantial economic impact [1]. The breeding of Norway spruce (Picea abies (L.) Karst.) has taken place since the middle of the 20th century due to the financial importance of coniferous trees [2]. The profitability of long-term investment in genetic selection has been proven by increased productivity, quality and resistance to the risks of a changing environment, which has resulted in an improved monetary value [1][2][3][4][5]. The growing demand for wood-based biomass continues to maintain the significance of tree breeding in the future [1]. In addition, the financial performance of genetically improved material has been enhanced by the inclusion of carbon pricing [4], which indicates the additional benefits of breeding in terms of carbon sequestration.
When appropriately managed (e.g., in a sparsely planted plantation), Norway spruce can rapidly reach its target diameter [6], thereby shortening rotation time, reducing establishment costs and increasing financial outcomes [7,8]. In turn, a reduced rotation time can mitigate the risk of biotic and abiotic damage, such as that from wind, drought, root rot, bark beetle attacks, etc., which are commonly acknowledged threats for successful Norway spruce management in a changing climate [9][10][11][12][13].
Tree breeding has been estimated to advance production by 10-35% compared to unimproved material [14]. Substantial variations in genetic gains exist among the different improvement levels, including the combination of additive and non-additive genetic site index was 36.0 m. The mean annual temperature in the study area was circa +6.0 °C, with the mean monthly temperature ranging from −6.4 °C in February to +17.1 °C in July. The mean annual precipitation was circa 700 mm [43].
The plantation was established in 1964 using vegetatively propagated (grafted) planting material from 20 selected fast-growing plus trees of local origin, at a spacing of 5 × 5 m (400 trees ha −1 ). The rootstock material comprised seedlings of local origin. In total, 421 trees were planted in randomly distributed single-tree plots  replications (ramets) per clone). Weed control was carried out in the planting year and the first year after planting. No thinning was conducted prior to the sampling. No measurements had been performed before during the trial.
Both DBH and height were measured for all mature trees (i.e., 50 years old or more) in the plantation. The plantation growth and yield data for trees at the age of 50 years were published previously by Katrevičs et al. [6].
Cores from pith-to-bark at breast (1.3 m) height were collected at 5-mm increments from 221 50-year-old trees, which represented 19 of initially planted 20 clones (7-19 ramets per clone). Only straight trees with no visible crown asymmetry or other stem defects, such as visible root rot, double tops or severe browsing damage, were selected for core sampling. Annual ring width data were obtained using high-frequency densitometry with a LignoStation [44]. These were cross-dated and verified by graphical inspection and using the COFECHA software [45].
The tree ring data allowed for the study of the age-diameter relationship at an annual resolution [46]. The DBH for each year was reconstructed using annual diameter increments, which were based on the measurements that were taken at the age of 50 years. For the diameter reconstruction, the tree stems were assumed to be circular so that the DBH could be calculated from the incremental cores. As a result, time-series data for the DBH of each tree were obtained for the ages of 20-50 years ( Figure 1).

Modelling Approach
The clone effect on the parameters of the diameter-age relationship was investigated by applying the Chapman-Richard base equation: where DBH iA is the DBH for the i-th clone at age A, β 1 is the asymptotic diameter parameter, β 2 is the rate parameter, β 3 is the shape parameter and ε iA is a normally distributed zero-expectation random error due to the DBH that was observed at age A [47]. We followed the approach of using models that had been previously fitted successfully by other researchers [48]. The preliminary analyses and previous studies have suggested that the Chapman-Richard function is a biologically reasonable selection for modelling tree growth curves [34,46,49,50]. The clone effect on the diameter-age relationship was modelled as random effects on the parameters in Equation (1) using the function: where DBH iA and ε iA are as defined as for Equation (1), β 10 , β 20 and β 30 are fixed effect parameters and b 1i , b 2i and b 3i are random effect parameters for the i-th with no correlation among these parameters). Heteroscedasticity in the error term ε iA was not detected. The autocorrelations of the errors that were due to ring width measurements on the same trees were modelled by mixed first-order autoregressive and moving average structures (ARMA (1,1)). To remove the effects of the various variances in the different clones, we modelled different variances for each factor (clone) using the variance function [51].
Equation (2) was fitted using nonlinear mixed-effect regression with the nlme package in R software [52]. The effects of the clone on the model parameters were investigated using a common approach: eliminating the random effect parameters from the model one by one by using random effects for all parameters initially. Different combinations of parameters were tested and the best model was selected using the likelihood ratio test, Akaike's information criterion (AIC) and the Bayesian information criterion (BIC). The fitted model was evaluated using mean bias and root mean squared error [53]. The significance of random factors in the final fitted model was evaluated by applying likelihood ratio tests of model reductions for each random effect term, which were provided by the rand function in the lmerTest package [54]. The normality of the final fitted model residuals was tested by applying the Jarque-Bera test.

Genetic Parameters
To estimate the variance components, the following model was used: where y ij is the observation of the ith tree from the jth clone, µ is the overall mean, C i is the random effect of the cloning and ε ij is the random error. To evaluate the dynamics of the genetic parameters, (i) the broad-sense heritability (H 2 ) of DBH for each year and (ii) the age-age genotypic correlations of DBH among the years were estimated as follows: are the estimated variance components of clone and residual, respectively, and: whereσ 2 (Gx) andσ 2 (Gy) are the genotypic (clone) variances at two different ages and Cov (Gxy) is the estimated genotypic covariance between the two measurements.
For the final fitted model, we adapted a formula for calculating the genotypic coefficient of variation (CV g ) [55] in order to describe the extent of genetic variability among the clones for any random effect parameter: where CV g is the genotypic coefficient of variation,σ 2 b x is the estimated variance for the random effect parameter x and β x is the fixed effect parameter x. The CV g of DBH was calculated for each year. For comparison, the H 2 and CV g of tree height at the age of 50 years was estimated.
The realised genetic gain at the age of 50 years (final harvest age) was estimated for each clone using the best linear unbiased predictions (BLUPs) from the final fitted model as a percentage of the trial mean DBH.

The Growth Model
The model form with all three random parameters (random asymptote (b 1i ), random rate (b 2i ) and random shape (b 3i )) had the lowest value of the AIC statistic (−189.152) ( Table 1). Overall, the model fit was significantly improved by addition of the random rate and/or shape into the model that already contained the random effect parameter b 1i compared to the anamorphic random asymptote equation (p < 0.01). Table 1. Parameter estimates (fixed effect parameters: β 10 , β 20 and β 30 ; estimated variance components: and model statistics (AIC, Akaike's information criterion; BIC, Bayesian information criterion; logLik, log likelihood test) of the Chapman-Richard base equation when applying different combinations of random asymptotes (β 1 ), random rates (β 2 ) and random shapes (β 3 ), which were fitted to diameter at breast height from the 19 clones that were used in the study.

Parameter in Model
Parameter Estimates The residuals of the final model were distributed symmetrically around zero with an approximately even variance and did not violate the assumption of normality (p = 0.154). The model provided a good representation of the dataset ( Figure 2). The distribution of the random effects did not indicate any noticeable violations of the assumption of normality.
For the final fitted model, the mean bias was −0.57 cm with a root mean square error of 3.61 cm.  From the data that were used in this study, genetic variety (clone) significantly (p ≤ 0.05) affected the asymptotic DBH and the rate and shape parameters of the Chapman-Richards diameter-age model. Thus, there was a significant polymorphism among the diameter-age trajectories of the different clones.

Dynamics of Genetic Parameters
The estimated broad-sense heritability (H 2 ) ± standard error (SE) for DBH increased from 0.19 ± 0.087 at the age of 20 years to 0.35 ± 0.131 at the age of 40 years (Figure 3). Afterwards, H 2 gradually decreased and was 0.32 ± 0.123 at the age of 50 years. On the contrary, the genotypic coefficient of variation (CV g ) was the highest at the age of 20 years (12.9%) and decreased to 7.8% at the age of 45 years, after which it remained stable for the following 5 years. The H 2 and CV g for tree height at the age of 50 years was 0.41 ± 0.088 and 5.6%, respectively.

Dynamics of Genetic Parameters
The estimated broad-sense heritability (H 2 ) ± standard error (SE) for DBH increased from 0.19 ± 0.087 at the age of 20 years to 0.35 ± 0.131 at the age of 40 years (Figure 3). Afterwards, H 2 gradually decreased and was 0.32 ± 0.123 at the age of 50 years. On the contrary, the genotypic coefficient of variation (CVg) was the highest at the age of 20 years (12.9%) and decreased to 7.8% at the age of 45 years, after which it remained stable for the following 5 years. The H 2 and CVg for tree height at the age of 50 years was 0.41 ± 0.088 and 5.6%, respectively. The estimated CVg for the DBH growth model parameters β1, β2 and β3 were 11.0, 17.1 and 11.9%, respectively. Overall, variance in the asymptotic DBH due to the clonal effects was slightly greater ( = 26.69 cm) than the within-group (within-clone) error variance (σ ԑ = 23.05 cm) ( Table 1). The realised genetic gains of the clones at the final harvest age varied from −6.3 to +24.0% around the trial mean ( Figure 4).
The estimated age-age genotypic correlations (rG) were positive and mainly very strong (>0.80), although they were slightly weaker, yet still strong (>0.76), for the older trees (age 20-22 versus age 46-50 years). There was a trend of slightly stronger correlations between similar ages with increasing age. For instance, DBH between age 45 and 50 was 100% genetically correlated, while rG between ages 20 and 15 was 0.97. The estimated CV g for the DBH growth model parameters β 1 , β 2 and β 3 were 11.0, 17.1 and 11.9%, respectively. Overall, variance in the asymptotic DBH due to the clonal effects was slightly greater (σ 2 b 1 = 26.69 cm) than the within-group (within-clone) error variance (σ 2 ε = 23.05 cm) ( Table 1). The realised genetic gains of the clones at the final harvest age varied from −6.3 to +24.0% around the trial mean ( Figure 4).
The estimated age-age genotypic correlations (r G ) were positive and mainly very strong (>0.80), although they were slightly weaker, yet still strong (>0.76), for the older trees (age 20-22 versus age 46-50 years). There was a trend of slightly stronger correlations between similar ages with increasing age. For instance, DBH between age 45 and 50 was 100% genetically correlated, while r G between ages 20 and 15 was 0.97.

Discussion
Growth curve fitting is a common procedure that is used to understand general biological growth trends [56] and incorporating factors that affect this curve can improve the accuracy of the function under different conditions [33,57]. These factors include any silvicultural treatment over time [58], among which clonal effects are perceived. In our study, the Chapman-Richards equation, with its sigmoid form, inflection point and asymptote, fitted the data well ( Figure 2) and the function parameters were meaningful for the analysis of the response of DBH to clonal effects over time [58,59].

Discussion
Growth curve fitting is a common procedure that is used to understand general biological growth trends [56] and incorporating factors that affect this curve can improve the accuracy of the function under different conditions [33,57]. These factors include any silvicultural treatment over time [58], among which clonal effects are perceived. In our study, the Chapman-Richards equation, with its sigmoid form, inflection point and asymptote, fitted the data well ( Figure 2) and the function parameters were meaningful for the analysis of the response of DBH to clonal effects over time [58,59]. Applying a mixed model, we estimated the population mean response that was common to the entire clonal plantation as fixed effect parameters [34], while random effect parameters that varied around the fixed mean [60] represented the specific responses of the clone grouping variable (Table 1, Figure 4). The possibility of estimating the variance within and among the clones for the random effect parameters [51,61] was an important advantage in terms of analysing the growth modelling results from a tree breeding perspective, since this field of study widely utilises genetically determined variations in traits of interest (in our case, the model parameters) to estimate quantitative genetic parameters and breeding values [55].
Our final model fitting results showed that all three model parameters (asymptote, growth rate and shape) were significantly affected by clone. There are not many previous studies on this issue, but those that do exist have mainly explored the tree height growth response to genetic effects. For loblolly pine (Pinus taeda L.), Sabatia and Burkhart [62] found that clone significantly affected the asymptotic height and shape parameters of the Chapman-Richards height-age model, while Knowe and Foster [63] concluded that half-sib families had an effect on the asymptote and rate. Half-sib genetic varieties and provenance have been reported to only affect the asymptote of the Korf function [29,64]. In contrast, Sprinz et al. [65] only reported the effects of half-sib families on the shape parameter. In Chinese fir (Cunninghamia lanceolata (Lamb.) Hook.), height growth trajectories in stands with different provenances only differed in the asymptote [66]. A strong genetic influence on the asymptote has also been reported for Konishii fir (Cunninghamia konishii Hay.) using an application of the Weibull-based function of DBH [67].

Dynamics of Clone-Specific Diameter Growth and Its Genetic Parameters over Time
Although variance component analyses are of primary interest in quantitative genetic studies, the random effect variances in the model parameters have been rarely reported or interpreted within ecological modelling [68]. Our results of the random effects showed significant variances in all three model parameters that were determined by the genetic effects. The overall magnitude of variation in DBH among the clones was greater than that within the clones (Table 1), which indicated the great importance of random effect variance, i.e., clonal differences, for selection [61,68,69]. Moreover, the genotypic coefficients of variation for the model parameters (asymptote, growth rate and shape) ranged between 11.0 and 17.1%, which exceeded the variation in DBH at the time of final felling (7.8%) and suggested the potential for genetic improvement in the asymptotic DBH, as well as a more rapid radial growth trajectory. For DBH growth function of Konishii fir, provenance and family have been reported to account for 26% of the variation in the asymptote, but only 4% in both the rate and shape [67]. For the loblolly pine height curve, the coefficients of variation among the clones in the asymptote and shape of the Chapman-Richards function have been estimated to be 4.2 and 2.8%, respectively [70].
By comparison, the CV g for DBH in trials involving 19-year-old Norway spruce clones has been reported to vary between 13.6 and 15.9% [71]. Still, our results showed that CV g and H 2 were not constant and changed over time. The estimated H 2 increased from 0.19 at the age of 20 years to a peak of 0.35 at the age of 40 years. A similar trend has been observed for the narrow-sense heritability h 2 of DBH in open-pollinated progenies of Norway spruce in southern Sweden, for which h 2 stabilised at circa 0.22 at the cambial age of 10 years after increasing from close to zero near the pith [72]. For Scots pine full-sib families in Sweden, a trend of increasing h 2 has been observed for wood quality traits and cumulative ring width, for which this parameter increased from very low to 0.25 at the age of 20 years [73]. Hannrup et al. [71] reported a H 2 of 0.34-0.50 for 19-year-old Norway spruce clones, while a low H 2 (< 0.14) was estimated for juvenile white spruce (Picea glauca (Moench) Voss) somatic clones 4 years after outplanting [74].
In contrast to H 2 , CV g decreased by almost half from 12.9% at the age of 20 years to 7.8% at the age of 45, after which it remained stable (Figure 3). Similar to descendent trend that was observed for CV g in our study, spruce clones in series of experiments in northern Germany demonstrated a steady decrease in genetic variance in height from circa 20% at the age of 3 years to a plateau of 7% after the age of 8 years [75]. Joint site data from Sweden showed that the coefficient of additive genotypic variation CV a for DBH in Norway spruce decreased from 15.45% to 11.91% at the ages of 12 and 21 years, respectively [76]. The estimated CV a showed a marked decline with age for both H and DBH (from 25-33% to 7-14%) in Silver fir (Abies alba Mill.) up to the ages of 10-15, after which they became stable [77].
In our study, H 2 and CV g stabilised after the age of 40-45 years (Figure 3), which was likely due to reaching the moment of canopy closure and the intensification of inter-tree competition [78] as diameter increment is considered to be the trait that is most (first) affected by competition [79]. Canopy closure depends on the initial planting density and the increment of the trees. The onset of inter-tree competition in loblolly pine has been estimated to be after 5 years in the most densely planted plots (1.2 × 1.2 m) and after 8.6 years in the most sparsely planted plots (3.7 × 3.7 m) [80]. The latter age corresponds to circa 35% of the rotation age for managed P. taeda [81], which indicates a longer period of competition-free early growth for the wider spacing. Although canopy closure in Norway spruce progeny plots in northern Europe with conventional spacing is typically observed at the age of 10-20 years [79], the low planting density (5 × 5 m) in the study site could have delayed it substantially [82,83].
Nevertheless, we acknowledge that the exceptionally sparsely planted trial is an unusual competitive environment for plantations. Some of the commonly known risks for fast growing sparely planted Norway spruce include reduced wood density and stronger branching, which reduce timber quality [71,72,84]. Still, despite the unfavourable negative phenotypic correlations between fast growth and wood quality traits, selection for both may be possible at a clonal level [12,85]. Our study site also differed from common Norway spruce plantations in terms of its reproductive material (vegetatively propagated grafted clones), which is not commonly used as planting stock. Liziniewicz et al. [39] concluded that growth differences among different genetic varieties can be assumed to be independent of plant type, although they compared seedlings to rooted cuttings. The rootstock × scion effect has been reported as being inessential to the growth of loblolly pine, which is contrary to the genetic and site impacts [86]. The fast growth of the studied plantation indicated no potentially negative effects of cyclophysis [87][88][89][90], which has occasionally been detected in clonal trials [91]. However, we stress that our interest was in fitting the model as an analytical tool to investigate the differences in the response curve trajectories of different clones of the same propagation type with a wide range of realised genetic gains for DBH ( Figure 4). Still, the clones characterise a local population and thus, had to be generalised with caution.

Age-Age Genotypic Correlations
The growth pattern of each individual tree is shaped by a combination of both its genetics and environment [92]. Still, the DBHs of individual trees are unlikely to be affected by competition at a genetic level over time from the indirect additive effect of neighbouring trees [79], which was indicated by the strong to very strong positive genotypic age-age correlations that were observed in this study (r g > 0.76).
Somewhat stronger r G at older ages have also been reported previously for various conifer tree species, which could be associated with the cumulative nature of the growth traits [93,94]. Our results corresponded well to earlier studies of growth traits in progeny trials, which showed very strong genotypic relations among similar ages, with a slight decline as the age differences increased [72,75,77,78,94].
Overall, the strong age-age correlations and clone-specific DBH growth trajectories, which were mostly without pronounced rank shifts during the studied period (Figure 1), justify the early selection for DBH. We observed that the clone-specific DBH growth curves are primarily relatively parallel to the population mean curve ( Figure 4) and that there were no trends of superior clones, in terms of asymptotic DBH, having a more rapid growth rate.
Nevertheless, all DBH growth curve parameters showed substantial genotypic variance (CVg ≥ 11.0%, Table 1) and separate clones (e.g., clone numbers 19, 24, 53) had obviously different growth shapes and rates (Figures 1 and 5), which were likely associated with the different genetically determined growth responses of certain clones to climatic factors [95].
Forests 2022, 13, x FOR PEER REVIEW 11 of 15 growth rate. Nevertheless, all DBH growth curve parameters showed substantial genotypic variance (CVg ≥ 11.0%, Table 1) and separate clones (e.g., clone numbers 19,24,53) had obviously different growth shapes and rates (Figures 1 and 5), which were likely associated with the different genetically determined growth responses of certain clones to climatic factors [95]. Thus, we argue that the selection of superior clones using early measurements of DBH could be combined and improved with information about the genetic variations in growth curve parameters when such data (e.g., increment cores) become available from long-term trials. In particular, intensive tree breeding that applies vegetative propagation would benefit from the knowledge of genotype-specific growth trajectories when a set of only a few superior clones is selected. Precise growth functions can improve the effectiveness of tree breeding, which supports the early selection of superior genetic entries using more reliable information on their expected future performance and the subsequent economic outcome of their selection [62,63]. Even small differences in the predicted growth traits can have a substantial impact on the ranking of clones [63]. The knowledge of genetic variances in function parameters can be utilised to develop and test more practically desirable and dynamic base-age-invariant functions [96] by incorporating genetic effects as modifiers for the model parameters [97], which is just as important in practical applications in forestry.

Conclusions
We found different growth patterns among the studied Norway spruce clones. There were significant differences between clones in terms of all three of the parameters of the Chapman-Richards function: asymptote, rate and shape. All diameter growth curve parameters demonstrated genotypic variations between the clones that were sufficient for selection with a wide range of realised genetic gains in DBH (up to +24%) that had the potential for genetic improvement. The studied clones possessed considerable heritability Thus, we argue that the selection of superior clones using early measurements of DBH could be combined and improved with information about the genetic variations in growth curve parameters when such data (e.g., increment cores) become available from long-term trials. In particular, intensive tree breeding that applies vegetative propagation would benefit from the knowledge of genotype-specific growth trajectories when a set of only a few superior clones is selected. Precise growth functions can improve the effectiveness of tree breeding, which supports the early selection of superior genetic entries using more reliable information on their expected future performance and the subsequent economic outcome of their selection [62,63]. Even small differences in the predicted growth traits can have a substantial impact on the ranking of clones [63]. The knowledge of genetic variances in function parameters can be utilised to develop and test more practically desirable and dynamic base-age-invariant functions [96] by incorporating genetic effects as modifiers for the model parameters [97], which is just as important in practical applications in forestry.

Conclusions
We found different growth patterns among the studied Norway spruce clones. There were significant differences between clones in terms of all three of the parameters of the Chapman-Richards function: asymptote, rate and shape. All diameter growth curve parameters demonstrated genotypic variations between the clones that were sufficient for selection with a wide range of realised genetic gains in DBH (up to +24%) that had the potential for genetic improvement. The studied clones possessed considerable heritability (H 2 = 0.32) and a significant genotypic coefficient of variation (CV g = 7.8%), which tended to stabilise at the age of 40-45 years. Although the mainly very strong age-age genotypic correlations would justify selection for DBH at earlier age, we suggest that growth curve parameters are just as important for tree breeding in terms of the selection of elite clones, once the data from long-term clonal tests are available. The substantial genetic variation in growth rate, shape and asymptote suggests the potential for more precise selection using predictions for not only final dimensions, but also desirable patterns of growth trajectories. Based on this information, clone-specific genetic modifiers should be tested using dynamic base-age-invariant functions for future growth predictions in practical forestry in order to improve the prediction accuracy for genetic entries with various genetic gains.