The provenance family trials of B. alnoides were established at Mengla and Pingxiang in 2002, Hua’an in 2003, and Changning in 2007 ( Figure 1 ), where 400, 386, 280, and 250 half-sib families were involved, respectively, all from 25 provenances. Their seeds were collected from the natural forests in Yunnan and Guangxi, China, and their seedlings were raised separately at a nursery near each site. Specimens of each family were made when sampling in the natural forests, and species identification was done by Dr. Jie Zeng. Voucher specimens of 25 provenances were deposited at the Herbarium of the Research Institute of Tropical Forestry, Chinese Academy of Forestry, Guangzhou, China, and their information were shown in Supplementary Table 1 . The information of geographic locations and climatic conditions at the four sites and geographic locations of 25 provenances could be seen detailedly in our previous study (Yin et al., 2019). The seedlings were planted in randomized complete block design with single individual plots and 12 to 19 blocks at a spacing of 2 m × 3 m.

Tree growth was investigated at four sites in late 2017 to early 2018. The blocks with survival rates lower than 60% were discarded at each site so as to keep consistency of survival rate among provenances and sites and thus avoid the impact of mortality on the H-D relationship. The current survival rate ranges from about 63.16% to 66.61% at four sites. All trees in well-reserved blocks were measured for DBH to the nearest 0.1 cm using a diameter tape and for H to the nearest 0.1 m using a Vertex IV Altimeter (Haglöf Sweden AB, Västernorrland, Sverige). A total of 7,724 trees were used for model analysis, including 3,403, 953, 1,438, and 1,930 trees at the ages of 15, 15, 14, and 10 years in Mengla, Pingxiang, Hua’an, and Changning sites, respectively. Variance analyses (ANOVA) and Tukey’s multiple range tests were conducted to estimate the differences in DBH and height among sites and provenances. The scatter plots of H and DBH showed a curvilinear relationship at all four sites ( Figure 2 ).

Ten commonly used H-D equations ( Table 1 ), including two linear and eight nonlinear ones, were used to model the H-D relationships. For evaluating the model performances, the goodness of fit for each model was then assessed with mean absolute error (MAE), Akaike’s information criterion (AIC), adjusted coefficient of determination (R²), largest log-likelihood (LL) and paired t-test (Liu et al., 2017). The equation with a high goodness of fit was selected as the base model for further analysis.

Site and provenance variables usually affected the asymptotes rather than the slopes of the H-D equations (Buford, 1986; Fu et al., 2016; Zhang et al., 2019); they were thus introduced step by step into the asymptote parameter of the base model in the present study. Here, the extended models for Equations (1) and (8) in Table 1 were taken as examples for linear and nonlinear equations, respectively:

where H and DBH are the tree height and stem diameter at breast height, respectively; φ 0 , φ 1 , and φ 2 , are basic parameters; x i and k j are dummy parameters; and ϵ is error term. S i denotes the site dummy variable with i = 1, 2, and 3. Changning site is taken as a reference, and Mengla, Pingxiang, and Hua’an are set as dummy variables: S 1 = 1 denotes Mengla site and 0 for the rest of sites; S 2 = 1 denotes Pingxiang site and 0 for the rest of sites; S 3 = 1 denotes Hua’an site and 0 for the rest of sites; and the Changning site is represented by S 1 = S 2 = S 3 = 0. P j denotes the provenance dummy variable with j = 1, 2, 3, …, 24. Provenance Y is taken as a reference, and other provenances are set as dummy variables: P 1 = 1 denotes provenance A and 0 for the rest of provenances; P 2 =1 denotes provenance B and 0 for the rest of provenances; and so on…. P 24 = 1 denotes provenance X and 0 for the rest of provenances; and the provenance Y is represented by P 1 = P 2 = ⋯ = P 24 = 0. The likelihood-ratio test (LRT) (Fang and Bailey, 2001) was used to estimate whether the site or provenance variables influenced the equation parameters significantly:

where LL i and LL j are largest log-likelihood of extended model i and base model j. Compared with the critical values for chi-squared distribution, if X i − j 2 > X ( n i − n j , 0.05 ) 2 , there exist significant differences between model i and model j at 0.05 level. The model parameters then are influenced significantly by site or provenance variables, where n i and n j are the degrees of freedom of extended model i and base model j, respectively.

Block was then introduced into the dummy model as random effects to analyze the H-D allometry. A method of Davisian and Giltinan (1995) was used to account for the within-tree heteroscedasticity and autocorrelation in the variance–covariance matrix ( R i j l m ) of the error term ( ϵ i j l m ):

where R i j l m is a variance–covariance matrix of the error term ( ϵ i j l m ) in the mth tree nested within lth block for the jth provenance at the ith site; σ 2 is a scaling factor for error dispersion, given by a value of residual variance of the estimated model; and R i j l m and Г i j l m are n i j l m × n i j l m diagonal matrix explaining the variance of within-tree heteroscedasticity and autocorrelation structure of errors. Power variance fuction (Fu et al., 2016) was used to reduce heterogeneity in variance:

where x is fitted values as the selected predictor, γ is the parameter to be estimated, and ϵ ′ is an error term. The all parameters in the nonlinear mixed-effects (NLME) models were estimated with restricted maximum likelihood implemented in R software “nlme” package (Pinheiro et al., 2017).

On the basis of NLME models, the curves of the H-D relationship at four sites and of 25 provenances at each site were simulated to assess variation among sites and provenances. Considering the visualization and simplicity for the H-D relationship in line graphs as well as subsequent provenance selection, the 25 provenances were clustered into several groups through the system clustering method based on Euclidean distance of asymptote parameter k j at each site. The H-D relationship of groups was then showed in a line graphs and used for comparison and selection for provenances.

The H-D relationship was further involved in the selection of excellent germplasm. Because the large-sized timber production is the target of management for valuable tree species and lower asymptote of the H-D relationship at a given DBH is important for wood quality and tree stability, above-average DBH growth and lower asymptote of the H-D curves were thus used as the indicators. The provenances with above-average DBH were ranked by their asymptote parameter in the H-D equations, and excellent provenances with lower values of asymptote parameter were then selected with a rate of 20% at each site, i.e., five provenances of the total 25 were selected. The individual volumes were used to evaluate the effect of selection and were calculated as volume = 0.45/4πDBH² × H (Wang et al., 2013), and selection gains were calculated as selection gain = (selected mean − total mean)/total mean × 100%.

All data analyses and modeling were performed using R software, SPSS 13.0, and Microsoft Excel 2010.

There were significant differences in DBH and H of B. alnoides among the four sites at 0.05 level ( Table 2 ). The DBH and height performed better at Mengla and Hua’an than at the other two sites. At each site, both traits differed significantly among provenances, and well-performed provenances varied with sites, which demonstrated obvious interaction between sites and provenances. Provenance A, G, and B performed the best in both growth traits at Mengla, Hua’an, and Changning, respectively; and provenance V and R performed the best in H and DBH at Pingxiang, respectively.

All parameters estimated in each H-D candidate function were of significant difference from zero at 0.05 level ( Table 1 ). The paired t-test showed that there was no significant difference between observed values and predicted values of each model except model (10). The MAE and AIC of the Weibull model [Equation (8)] were the lowest, and its coefficient of determination (R²) was the greatest among all equations ( Table 3 ). Weibull model also satisfied the biological assumption that H is 1.3 m (breast height) when DBH equals to zero for the H-D relationship (Bi et al., 2012), it was thus selected as the base model.

Site and provenance variable was introduced as dummy variables step by step into the asymptote parameter of the base model [Equation (8)], and the random parameters of block were then added to the asymptote and slope of DBH to develop site-level NLME model [Equation (16)] as follows:

Considering the interaction between sites and provenances, provenance dummy variables were introduced into the base model [Equation (8)] at each site, the provenance-level NLME model was then shown as follows:

where μ 0 l and μ 1 l are the random effects on asymptote and slope caused by the lth block, and μ 0 l ~N(0, σ 0 b l o c k 2 ), μ 1 l ~N(0, σ 1 b l o c k 2 ).

The simulation results and the model performance of the extended models are shown in Tables 4 , 5 , respectively. The performance of site-level NLME model [Equation (16)] and provenance-level NLME model [Equation (17)] were significantly improved with a smaller MAE and AIC values and a larger R² and LL values than the base model [Equation (8)]. The LRT test results also demonstrated that asymptote parameters of the models were influenced significantly by site variables and provenance variables at four sites (P< 0.05). The residual scatter diagram of Equation (8) inferred obvious increase trend of predicted H with increasing DBH. Because the heteroscedasticity was effectively accounted by the power variance function [Equation (15)], this trend disappeared with NLME models [Equations (16) and (17) at four sites] ( Figure 3 ).

All the basic estimated parameters ( φ 0 , φ 1 , and φ 2 ) of the models were significantly different from zero (P< 0.001), and the parameters x 1 and x 2 were also significant (P< 0.001 and P< 0.05) rather than x 3 ( Table 4 ), indicating that the trees at Mengla and Pingxiang sites were significantly taller than those at Hua’an and Changning sites for a given DBH. The H-D curves at four sites by Equation (16) are shown in Figure 4 . The simulation results of provenance-level NLME model [Equation (17)] indicated that the influences of provenance variables on asymptote of the models differed among four sites, and difference of asymptote parameters ( k j ) among provenances was more significant at Mengla than other three sites. To demonstrate the variance of the H-D relationships clearly, the provenances were clustered into three to five groups based on asymptote parameters (k_j ) of Equation (17) at each site ( Supplementary Figure 1 ; Table 6 ). The simulation plots of the H-D curves for these groups were shown in Figure 5 . The H increased rapidly with increasing DBH at Changning site, whereas it gradually increased slowly at the other sites, and it tended to be stable at Mengla site. The differences of asymptote among provenances increased with increasing DBH and were not obvious at Changning site.

The results of elite provenances selection showed that provenances B, W, P, J, and N performed well at Mengla; R, P, U, E, and W did well at Pingxiang; Q, R, U, W, and G did well at Hua’an; and B, V, T, W, and P did well at Changning ( Table 6 ). Their asymptotes parameters values ( k j ) and H-D ratio decreased from 7.17% to 486.05% and from 3.07% to 4.72% at four sites, respectively, and the selection gains of H, DBH, and individual volume ranged from −1.38% to 5.18%, from 1.46% to 9.87%, and from 1.99% to 29.81%, respectively.

The application of sigmoidal nonlinear models with inflexion points and asymptotes in simulation height–diameter (H-D) relationship has a long history in tropical and subtropical zones, and many nonlinear models have been developed (Scaranello et al., 2012). The study of Scaranello et al. (2012) indicated that Weibull and Chapman–Richards nonlinear models could improve the fitting accuracy of the H-D relationship in biomass estimates for tropical Atlantic forests in southeastern Brazil. Linear models are also suitable tools once a linearizing relationship can be established through data transformation and reasonable biological significance can be interpreted (Curtis, 1967; Watt and Kirschbaum, 2011). One typical example was that Buford (1986) used transformed linear models to fit the H-D relationship of loblolly pine and found out that the R² values were consistently larger than 0.88, which provided an appropriate and accurate simulation of height–diameter curves. In the present study, Weibull nonlinear model [Equation (8)] showed the best goodness of fit (R^2 = 0.837) among 10 candidate models when fitted with the dataset of H and DBH and was selected as a base H-D model for B. alnoides.

Whether based on linear or nonlinear models, the way to further improve the goodness of fit is to introduce and explain more effect variables (Zhang et al., 2019). The mixed model and dummy variable model approaches are frequently used to construct the H-D relationship equations for tree species. Mixed-effects models can isolate fixed and random variations and are applied to incorporate complex nested stochastic structure in the H-D curves. For instance, Zhang et al. (2019) introduced climate, site, plot variables, etc., into the NLME model to reveal the influence factors on the H-D relationship of Chinese fir. Dummy variable models are also of adequate methodology, especially for the quantification of categorical variables. For an example, Buford (1986) introduced locations into the linear model as dummy variables and revealed clearly the variation of asymptotes in H-D curves among locations.

For the variation analysis of multi-layer complex variables, both dummy variables and mixed-effects models are all powerful modeling tools. There is a little difference between the two kinds of models when the number of samples in each category is large (Wang et al., 2008; Fu et al., 2012). In the present study, the dataset contained 953–3,403 individuals at four site with average of 38–136 individuals per provenance at each site ( Table 2 ), and both dummy variables model and mixed-effects model were thus applied comprehensively. Sites and provenances were introduced into the base model as dummy variables, and block was introduced as random effect to develop the nonlinear mixed-effects models. Significant variations of the H-D relationship were then observed among sites and provenances at each site, and the clear visualization of variation for the H-D relationship among provenances would provide a great convenience for the following selection breeding for B. alnoides.

A better performance of height growth and higher asymptote of the H-D curves was observed, and the H was infinitely close to the growth limit at Mengla, indicating that better soil nutrients and water conditions might be conducive to realise the genetic potential of H in the present study. A study on Norway spruce indicated that abundant soil water and nutrients were beneficial to height growth and thus had a high asymptote of the H-D curves (Bošeľa et al., 2014). In the study of Trouvé et al. (2015) on Quercus petraea Liebl. and the study of Kempes et al. (2011) on the forests in the continental United States, H growth usually decreased and more resources were assigned into radial growth with increasing water stress, which could also result in a low asymptote of the H-D curves. Zhang et al. (2019) reported that the H of Chinese fir plantations decreased with increasing annual heat moisture index in the subtropical zone. These also verify that sufficient water and nutrients can promote H growth and affect the H-D curves. The study of Homeier et al. (2010) on Ecuadorian montane rain forest and the study of Bošeľa et al. (2014) on Norway spruce in the Western Carpathians demonstrated that H tended to decrease with increasing altitude, and DBH gradually played a more important role than H in the H-D relationship. In the present study, the H-D relationship did not show a clear trend with such a large range variation of altitude (275–1250 m) at four sites. This might be closely dependent on the complex geomorphic features in southern China. Buford (1986) found that the asymptotes of the H-D curves for loblolly pine were different across six sites based on a dummy variable model. The site variables used in the present study are also a general one, and specific variables of site conditions should be included to improve the model accuracy and reveal influencing factors in further studies.

The variation of the H-D relationship among provenances was not consistent at four sites, demonstrating that there were interactions between provenances and sites. This genotype by environment interaction of other growth and quality traits was also observed in our previous study (Yin et al., 2019). The differences among provenances were modeled separately for each site to deal with this interaction in the present study. The specific expression of genetic effect was that different asymptotes and similar slopes of the H-D curves were seen among provenances at each sites. This was following the study of Buford (1986) on loblolly pine at the age of 15 years, in which six of the nine sites showed a provenance effect on the asymptotes of the H-D curves, whereas the study of Egbäck et al. (2015) on loblolly pine at the age of 6 years indicated that there were no significant differences among one seed orchard mix and six families for both the asymptote and slope parameters from the Korf function. This inconsistency of genetic variance may be due to the facts that the differences are larger at provenance level than those at family level, and more genetic differences of the H-D relationship are expressed at the age of 15 years than those at the age of 6 years. In the present study, the variance of asymptotes of the H-D relationship among provenances was relatively smaller in Changning site at the age of 10 years compared with that in other sites at the age of 14–15.

High selective gain is the core purpose of the selective breeding process. In the present study, five excellent provenances were selected on the basis of the lower asymptotes [k_j in Equation (17)] of the H-D curves and DBH above mean values at each site. The selection gains of DBH, height, and individual volume ranged from 1.46% to 9.87%, from −1.38% to 5.18%, and from 1.99% to 29.81%, respectively, for four sites in the present study, whereas they were 7.31%, 4.57%, and 14.84% based on growth and quality traits in our previous study (Yin et al., 2019). Although DBH, height, and individual volume varied among the sites, the maximum selection gains of them in the present study were all higher than those in the previous study. The selection gains of DBH, height, and individual volume at Mengla were lower than those at the other three sites. This could be attributed to the trade-off between lower asymptote and biger DBH of provenances at Mengla site, in which the fertile soil and suitable climate formed biger DBH and also higher tree hight and asymptote of provenances.

Although there is a statistical correlation between DBH and H, the disproportional selection of DBH and H will still occur in the selection process. For example, Yin et al. (2019), using biplot analysis for main genotypic effects and G × E interaction to select the superior provenances of B. alnoides, found that the selection gains of DBH (7.31%) were much higher than that of H (4.57%), whereas Wang et al. (2017) observed that the superior clones’ selection gain of H (14.84%) was much higher than that of DBH (3.95%) during selecting superior clones of B. alnoides based on simple index method. When selecting through volume alone in the present study, the selection gain of DBH (6.0%–10.8%) was higher than that of H (5.3%–8.1%) ( Supplementary Table 2 ). These inferred that excellent germplasms with high selection gains in the H-D relationship could not be obtained through disproportional selection of DBH and H.

The selection results of DBH and H can affect the H-D ratio, which is related to the stem quality and stand stability (Coutts, 1986; Sharma et al., 2016; Zhang et al., 2020). Schmidt and Kändler (2009) once found that the individuals of the higher site index had the lower H-D ratios for Norway spruce; and for a given site index, the lower the H-D ratios, the higher the stem quality. Tree pulling experiments conducted by Peltola et al. (2000) also showed that individuals of birch species (Betula spp.) with the lower H-D ratios had the bigger resistive bending moment for stem breakage. Therefore, excellent germplasms with lower H-D ratios have more production value in plantation management. Moreover, H-D ratio is a factor which should be taken into account during selecting excellent germplasms. Because of the disproportional selection of DBH and H, the excellent germplasms usually had a higher H-D ratio than the unselected germplasm in previous studies. Andersson et al. (2007) and Kroon et al. (2008) all observed that the selection for improved growth, based on H, resulted in trees with a higher H-D ratio. Similarly, Sabatia and Burkhart (2013) also found that the average H of the improved loblolly pine clones increased, whereas the average DBH did not, because the H-D relationship was not considered. Therefore, they also pointed out the high necessity of introducing the H-D relationship into the selection methods system of superior germplasm.

In the present study, the H-D ratios of the superior germplasms reduced from 3.07% to 4.72% at four sites based on selection by asymptote parameter (k_j ) and DBH. While after calculating the selection gains of the H-D ratios in the excellent germplasms of our previous study (Yin et al.2019), we found they were −1.7%, −0.8%, 2.5%, and −0.5% in Mengla, Pingxiang, Hua’an, and Changning, respectively (not published). When using volume for selection alone in the present study ( Supplementary Table 2 ), the modified direction of the H-D ratios was also unstable, which increased at Mengla site with the percentages of 0.84% but reduced from 0.39% to 3.00% at other three sites due to improvement prior to DBH. Their asymptote parameter (k_j ) was all increased significantly from 37.73% to 1023.22% at four sites. The inconsistency between H-D ratio and asymptote parameter (k_j ) was attributable to the fact that the H-D ratio is not constant over time and may include scaling effects (Kroon et al., 2008). These also proved the superiority of asymptote parameter (k_j ) in selecting superior germplasms by the H-D relationship. As a whole, selecting excellent provenances with the lower asymptote parameter (k_j ) of the H-D curves and above-average DBH was a logical and valid genetic improvement strategy. In addition, further studies should be carried out on linking stem form factor with the H-D relationship to predict volume gains more accurately.