Xu et al. (2016) reported a genetic linkage mapping study on growth traits in a full-sib family of Populus. This family was generated by crossing an eastern cottonwood (Populus deltoides) clone I-69, introduced to China in the 1970s (Wu et al., 1992), as the female parent, and a Euramerican poplar (P. × euramericana) clone I-45 as the male parent. This cross is equivalent to a backcross at the species level. This family was planted in 1984 for 24 years in a uniform site at Zhangji Forest Farm, Xuzhou, Jiangsu, China. Supplementary Table S1 shows the average of annual temperature and annual precipitation at the experimental site from 1987 to 2010. Stem heights (defined as the length of the main stem from the stem–root connection to the tip) and stem diameters at breast height were measured for each tree at the end of each growing season. A part of this family (64 full-sib members) was genotyped, producing 156,362 good-quality SNPs distributed on 19 chromosomes. Of these SNPs, 94,591 were backcross-like testcross markers and 61,771 were F₂-like intercross markers. For a testcross marker, one of the parents is heterozygous and the other is homozygous, and the intercross marker is derived from two heterozygous parents.

Stem growth in Populus trees experiences multiple phases of development including juvenile, mature, and senescence (Bond, 2000). In the study of Xu et al. (2016), only data from the first 14 years were considered, which include two possible phases. Because it is likely that tree growth experiences a distrain phase from year 14, a joint analysis of stem growth that spans24 years is essential to reveal the genetic architecture of phase change in stem wood growth trajectories. To do so, a more sophisticated model, as we will develop in this study, is needed.

Several classical nonlinear growth equations, such as those of Gompertz (Gompertz, 1825), Richards (Richards, 1959), Logistic (Verhulst, 1838), and Von Bertalanffy (Bertalanffy, 1957), describe the “S” shape of growth approximation. However, the growth of most organisms can actually be described as a composite form of multiple “S”-shaped phases because of seasonal fluctuations and differences in the growth rates of the components of organisms (Piantadosi, 1987). Koops proposed a multiphasic growth model, which is superior to single-phase growth in many studies and has been shown to more accurately estimate growth parameters (Koops, 1987; Grossman and Koops, 1988; Kwakkel et al., 1993). For example, Van der Klein et al. (2020) compared the monophasic, diphasic, and triphasic Gompertz and logistic models to describe the weight-age and gain-age functions of hens, and the results showed that the diphasic and triphasic Gompertz and logistic models yielded better fitting effects than the monophasic models. In addition, Treves et al. (2017) studied the multiphase growth pattern of the green alga Chlorella ohadii and metabolic changes during the growth phases. Multiphase analysis has proven to be a more accurate method for fitting biological data analysis and prediction.

According to the multiphasic growth view proposed by Koops (1987), the overall growth of height and diameter of poplar in the first 24 years is expressed as the sum of the growth function of two phases, juvenile and early adult (Xu et al., 2016), leading to a coupled nonlinear governing equation (NGE):

where H(t) and D(t) represent the growth of quantitative traits of height and diameter of trees at age t; H₁(t) and D₁(t) represent growth during the juvenile phase; and H₂(t) and D₂(t) represent the growth in adulthood.

In nature, the interactive relationship between diameter and height can promote or hinder the growth of another trait. The synergistic effect of stem height and radial growth is conducive to the survival and reproduction of trees. However, when resources are in short supply, height growth and diameter growth achieve a trade-off of resource competition (Hulshof et al., 2015). In particular, it is significant and obvious during the juvenile stage of vegetative growth and nutrient accumulation (Fu et al., 2017). The governing relationship between diameter and height during the juvenile stage can be expressed as:

where H₁ and D₁ represent the stem height and diameter of poplar in the first phase, which corresponds to juveniles. According to the form of the equation, the first phase can be divided into independent growth and interaction growth. α_H and α_D represent the independent growth rate of the first phase; K_H1 and K_D1 represent the asymptotic value of independent growth of two characters in the first phase; and β_H←D and β_D←H are dimensionless parameters used to describe the competitive or cooperative interaction between stem height and diameter.

However, when poplar reaches the maturity period of development, the growth patterns of diameter and height change, which can be described by growth equations without interacting parts (Xu et al., 2016). By evaluating the information criterion and fitting optimum of the equations through numerical experiments in Supplementary Table S2, we introduce the nonlinear governing equation (NGE), with H₂(t) and D₂(t) satisfying the following two equations:

where K_H2 and K_D2 are the growth asymptotic values of adulthood height and diameter; p_H and p_D are related to the initial values of adulthood; and q_H and q_D represent the growth rate of adulthood.

The NGE can not only clarify the dynamic change rule of the growth of each character at different phases but also represent the governing rule between two interaction traits in the first phase.

We consider the multiphasic growth of height and diameter as a whole according to the NGE and detect how QTLs controlled the overall growth curve from the perspective of system mapping, as well as the regulation of phase transition by QTLs (Gai et al., 2011; Sun and Wu, 2015). We design a model framework that takes n samples as the mapping population. The trait growth of each sample i was measured in a series of time points1, ⋯ , T, and y_1i = (y_1i(1), ⋯ , y_1i(T)) and y_2i = (y_2i(1), ⋯ , y_2i(T))(i = 1, ⋯ , n) represent phenotypic data of time-related traits 1 (height) and 2 (diameter), respectively. The composite growth characteristics composed of two traits approximately obey the bivariate normal distribution, where the time-dependent mean value and the symmetric covariance matrix (2T × 2T) are, respectively, expressed as

The elements on the diagonal inΣ are the variance matrix for each trait, and the nondiagonal elements are the covariance matrix between a pair of traits. The joint density function of the two-dimensional normal vectors is of the formf(y₁, y₂; θ), where θ represents the growth parameters. The joint probability density function of n sample trees constitutes the likelihood function as follows:

It is assumed that the multiphasic growth of the two traits is controlled by a set of QTLs located on the linkage map. J kinds of genotypes are on the assumed QTL, and the time-dependent mean values of the samples differ among genotypes:

The joint density function of the genotype is denoted by f_j(y₁, y₂; θ_j)(j = 1, ⋯ , J). On the assumed specific QTL, the number of samples with genotype j is n_j, satisfying ∑j=1jnj=n. The mixed likelihood function is expressed as

According to the multiphasic growth function of the NGE, we can estimate the mean vector through a set of parameters (α_jH, K_jH1, β_jH←D, α_jD, K_jD1, β_jD←H, K_jH2, p_jH, q_jH, K_jD2, p_jD, q_jD)(j = 1, ⋯ , J) instead of directly estimating 2T mean values. By comparing the genotype-dependent differences in this parameter set, we can determine whether the QTL affects multiphasic growth.

In addition, we use a highly efficient structured antedependence (SAD) (1) statistical model, which was proposed by Zhao et al. (2005), to represent the longitudinal covariance matrix Σ. This method uses a few parameters to calculate the matrix with a complex structure, which offers the advantages of simplicity and flexibility and greatly improves the computational efficiency and statistical ability of the QTL detection model. The innovation variance and the first-order pre-dependent parameter are defined as γ12, γ22  and ϕ₁, ϕ₂, respectively, to structure the residual variance of trait k (k = 1,2) at time t:

The covariance between t₁ and t₂ is expressed as

The fixed innovation variance between different time points is represented by parameter ρ, and the covariance between the two traits at different time points is expressed as

In the calculation, we maintain that the parameters of the covariance matrix under different genotypes of L₀ and L₁ are consistent, which can improve the operation efficiency of the model with little loss of precision. The maximum likelihood estimation of parameters in the mean and covariance structure of the QTL detection model is performed by applying the expectation maximization (EM) algorithm (Dempster, 1977). The fourth-order Runge–Kutta algorithm is used to calculate the parametric solution of the NGE, and the Nelder–Mead simplex algorithm (Zhao et al., 2004) is used as the optimization method to estimate the parameters of the nonlinear equation and matrix structure, making the model more efficient.

On the basis of the likelihood values L₀ and L₁, we detected the existence of QTLs that influence the multiphasic growth of traits by calculating the likelihood ratio statistic. By comparing the genotypic correlation differences within the parameter set, the existence of QTLs affecting the multiphasic growth of stem height and diameter was determined:

H₁: At least one of the above equalities does not hold true.

The null hypothesis means that growth is consistent across different genotypes, which could be described by a set of NGE parameters. The alternative hypothesis is that there are genotype differences in the NGE parameters. The form of the likelihood ratio statistic is as follows:

LR approximately obeys theχ²distribution and its degree of freedom is the difference in the number of model parameters between H₀ and H₁.

The rejection domain of the likelihood ratio statistic LR is taken as W = {LR ≥ c}, where the critical value c satisfies:

At the set test level α, if LR belongs to the rejection domain, we reject the null hypothesis H₀ and accept the alternative hypothesis H₁, which indicates that at this significance level, differences exist in the multiphasic growth of different genotypes of this marker. The p-values for hypothesis tests are converted by LRs and compared with the critical thresholds. We implemented stringent multiple testing procedures to control the false positive rate. On the one hand, we performed Bonferroni correction to adjust the critical value to ensure more rigorous results. On the other hand, a false discovery rate (FDR) correction based on the Benjamini and Hochberg method was performed on each p-value to control the proportion of false positives within a certain range.

We can further examine how QTLs control the multiphasic growth of the two traits and the characteristics in different phases. Based on the following hypotheses, we can detect whether the QTL regulates the interaction growth of stem height and diameter in the first phase:

and whether the QTL controls the growth of the two traits in the second phase by formulating the following hypotheses:

According to the biological significance of the parameter, we can also test whether QTLs control the asymptotic growths K_H1 and K_D1 of the independent growth of the two traits in the first phase, the total asymptotic growth with respect to the interaction, and the asymptotic growths K_H2 and K_D2 in the second phase. In addition, the independent growth rate of the first phases α_H and α_D, and the interaction relationships β_H←D and β_D←H of the first phase can also be tested. The initial growth values p_H and p_D and growth rates q_H and q_D of the second phase can also be tested. The phase transition time was determined by examining the differences in parameters with respect to genotypes.

We used the NGE (Equation 1) to fit the diameter and height growth data compared with monophasic, diphasic, and triphasic logistic models. In Figure 1, the fitting effects of the logistic models are good during the early stage of tree growth (approximately 7 years), but the NGE model has a better fitting effect during the later growth stage. This indicates that the whole growth trajectory of trees in different periods is changing. We compared the performances of these models by calculating the Akaike information criterion (AIC), Bayesian information criterion (BIC), Schwarz criterion, and adjusted R² (Table 1). As shown, the NGE performs better than the other growth equations, suggesting that it is crucial to consider the impact of interactive traits at different stages of growth. We also fitted the growth of 66 samples with the NGE, and the growth of both height and diameter for each hybrid showed excellent goodness of fit to the NGE (R² > 0.976; Supplementary Figures S1, S2). The residual errors of the growth data are distributed randomly over the predicted values (Supplementary Figure S3), suggesting that the NGE model is quite robust.

Obviously, the fitting results of the NGE can be divided into two parts; the first stage of which includes potential independent growth and interaction growth. Figure 1 shows the fitting results of the average growth curve. The first growth phase of stem height reached its asymptotic value at 7 years of age, and the age at which the asymptotic value of independent growth was reached was slightly later than that of the total growth of the first phase. Moreover, we found that in the first phase, the growth of stem height competed with diameter, leading to independent growth surpassing the total growth. In addition, the stem height growth of the second phase begins early, so the two phases overlap by a considerable period of time (Figure 1A). On the other hand, the stem diameter growth reached the asymptotic value of the first phase on the 10th year. The overall growth of stem diameter during this phase was higher than the independent growth curve, indicating that the growth of stem diameter was promoted by stem height. In addition, the second phase of diameter growth started later than that of the stem height and was far from asymptotic over the 24 years of the observed data (Figure 1B).

The joint process of multiphasic growth is regarded as a whole system, containing interactions between traits in the first phase and growth characteristics. By regressing the genotype-related growth trajectory of stem diameter and stem height, the differences in genotypes are expressed in the NGE parameters (α_jH, K_jH1, β_jH←D, α_jD, K_jD1, β_jD←H, K_jH2, p_jH, q_jH, K_jD2, p_jD, q_jD). The description of how a QTL affects the multiphasic growth and intrinsic interaction of traits can be characterized from genome-wide information based on a series of hypothetical tests. We performed FDR correction on the p-value results of all the SNPs and determined the critical thresholds as 10⁻⁴⁰ and 10⁻⁵⁰ after Bonferroni correction for testcross SNPs and intercross SNPs, which are strict detection levels. The Manhattan plots of the corrected p-values derived from test statistic values are shown in Figure 2. Our NGE-based mapping model identified 26 significant testcross SNPs and 39 intercross SNPs sporadically distributed over the genome, and 74% of significant SNPs were within candidate genes with known growth-related functions. For example, SNP 48,502 is within the GATA12 gene, which is a member of the GATA family of transcription factors. A set of significant SNPs are closely distributed in the same region of chromosomes, such as chromosomes 5, 8, and 9, residing within the same gene. For instance, SNPs 79,620, 79,624, and 79,625 are highly linked on chromosome 8, located within the region of the OASA1 gene, which encodes the synthesis of anthranilate synthase. Supplementary Table S3 presents basic information about all the significant SNPs, such as their chromosomal and physical locations, segregation types, allele types, p-values, and annotated candidate genes. Supplementary Table S4 displays specific genotype parameters of the significant SNPs. Supplementary Figure S4 depicts the genetic effect curves of the significant SNPs, showing the contribution to stem height and diameter growth. The heritability of stem height growth and stem diameter growth at different ages explained by each of the 65 significant SNPs is shown in Supplementary Figure S5.

Based on the 65 identified SNPs with significant effects on multiphasic growth, we estimated genotype-dependent growth parameters for two different phases and calculated the phenotypic variation explained (PVE) for 24 years by quantifying the dynamic genetic contribution of markers to growth. The PVE of these loci for stem height growth was significantly smaller than the PVE for diameter. A basic hierarchical cluster analysis was performed for the PVE of all significant SNPs over time, which was divided into two stages and represented by tags with different colors. Of the observed 24 years of growth, the first phase of stem height growth was the first 6 years, and the remaining 18 years belonged to the second phase. The first phase of diameter growth was the first 12 years, and years 12–24 represented the second phase (Figure 3). Compared with the diameter, the first phase of stem height growth was shorter, which is approximately consistent with the phase estimation results of the NGE growth curve shown in Figure 1 and which further suggests the validity of the NGE analysis.

Among the significant loci, the PVE of the same SNP for the two traits was different. For example, the genetic contribution of SNPs 17,787, 86,065, 73,524, 87,023, 28,821, and 102,765 to the growth of stem height was relatively substantial, while the genetic contribution of SNPs P28821 and 73,524 to the growth of stem diameter was the lowest among all the SNPs. In addition, SNP 77,477, with a low genetic contribution to stem height, had a PVE of 10% in diameter. The genetic contribution of some SNPs, such as 28,821 and 73,524, is mainly expressed in the growth of one trait. However, some SNPs, such as 17,787 and 102,765, make strong genetic contributions to both stem diameter and stem height, indicating that the same significant SNP affects the growth of tree height and diameter in different ways, thus revealing the pleiotropy of genes.

The growth of stem height and diameter is determined jointly by genetic effects of 26 testcross SNPs and 39 intercross SNPs as well as epistatic effects among them. According to the potential differential relationship of significant SNPs, we established a genetic network to characterize how these SNPs interact with each other to regulate the growth in diameter and height. In both the height and diameter network structures, most SNPs receive directional epistasis from other SNPs; they are activated or inhibited by other SNPs and play passive roles. Only a few significant SNPs called key QTLs were found to be in the pivotal dominant position of regulation in the genetic network, which had a strong potential to influence other QTLs with major outward links. Three SNPs, intercross SNP 48,833, intercross SNP 48,971, and testcross SNP 86,976, were the key QTLs in the stem height genetic network, and the directed links from them accounted for 27.17, 28.26, and 17.39% of all pairwise links, respectively (Figure 4A). The key QTLs in the stem diameter genetic network was testcross SNP 117,854 and intercross SNP 154,950, and links from the two directions accounted for 64.77 and 25% of the pairwise links, respectively (Figure 4B). Although these 65 significant SNPs jointly regulate the overall growth of diameter and height, they serve different roles in the structural networks for the growth of the two traits. The regulation relationship among SNPs is not singular in the network structure. For example, in the stem height network, there is reciprocal activation between SNP 86,139 and SNP 111,175. We discern a regulatory loop: SNP 48,833→SNP 86,065→SNP 48,726→SNP 48,502→SNP 48,833 in which one SNP regulates other SNPs but is influenced by another SNP. Details of the SNPs corresponding to each number are given in Supplementary Table S2.

The advantage of our NGE mapping approach lies in its capacity to analyze the genetic control of SNPs with significant effects on multiphasic growth, which can be tested by hypotheses based on related growth parameters. Here, we chose the representative testcross SNP 137,076 to interpret its genetic influence on growth in different phases of diameter and height. Both genotypes AG and AA in SNP 137,076 show that the growth of stem diameter benefits from stem height, and that stem height growth is deleteriously affected by diameter growth. This can be found from a similar form of overall growth and independent growth for stem height but greater overall growth relative to independent growth for stem diameter. The phase transition time of stem height occurs earlier than that of stem diameter. The two genotypes at this locus have similar control over the multiphasic growth of the two traits but different patterns in the growth details.

As shown in Figures 5A,B, during the juvenile phase, genotype AG had an extended period of duration of height growth compared with genotype AA (6.3 vs.5.3), and the difference was more pronounced for independent growth (8.3 vs. 6.8). Genotype AA has a larger asymptotic value (11.2 vs. 10.9) but a smaller asymptotic value for independent growth than genotype AG (16.3 vs. 16.6); thus, the latter suffers stronger inhibition (−6 vs. −5.4). The duration for the early adult phase of genotype AG was shorter (26.4 vs. 27.2), and the overlapping transition time between the two phases was also slightly shorter (5 vs. 5.1). According to the growth rate curve, the juvenile phase maximum growth rate of genotype AA was much larger than that of genotype AG (2.4 vs. 2), and there was a similar difference in the independent growth rate in this phase (3.6 vs. 2.7). The value of the early adult phase maximum growth rate was the same for both genotypes. As expected, the maximum growth rate and duration are determinants of overall growth. Different patterns in the onset, offset, duration, and maximum growth rate were found between the two genotypes for diameter growth (Figures 5C,D), suggesting that this SNP significantly affects not only the timing of growth but also the development between traits.

The influence of SNP 137,076 on the growth of the two traits is reflected in the temporal variation of the genetic effect, which also shows the characteristics of periodic change consistent with the growth curves. For stem height, the net genetic effect increases during the first 6 years and then decreases. During the juvenile phase, the overall genetic effect and the independent growth genetic effect show similar change patterns to the net effect, and then they level off at 12 years of age. The effect value of the independent growth is greater than that of the overall growth effect. The change in the interaction genetic effect shows a lag, beginning to rise at 5 years of age, and slowly declines until it stabilizes 12 years later. The genetic effect curve of the early adult phase gradually increased during years 1–12 and then decreased, which was related to the early growth of the second phase of stem height (Figure 6A). We analyzed a similar trend of genetic curves for diameter growth (Figure 6B). The pleiotropic effect between stem height and the diameter of SNP 137,076 is obvious with respect to the net genetic effect and juvenile phase genetic effect, but in the second phase, the curve shows two approximately flat trends, indicating that there is no obvious correlation between the growth of diameter and height in this phase. Thus, it is reasonable to describe the growth of the two traits without interaction in the second phase in the NGE (Figure 6C).

To verify further the statistical characteristics and mapping accuracy of the NGE, we conducted a simulation analysis based on the real growth data of tree diameter and height over 24 years under different sample size and heritability conditions. Among them, the sample size is the number of simulated samples; namely, n = 66 (the same as the real example), 100, and 200. Heritability is the proportion of genetic variance in the simulated phenotypic variance, and two levels of H²are utilized: 0.05 and 0.1. It is assumed that the growth of the two traits of trees is regulated by a set of QTLs in the genome. The simulation scale is 1,000 testcross genetic markers, and the phenotypic data are determined by one of the QTLs. In each simulation case, the proportion of simulation iterations of meaningful QTLs screened out of 100 repeated simulation experiments is the mapping accuracy (power), which is calculated and shown in Table 2. In the analysis of the multiphasic growth of tree traits, the mapping accuracy of QTL detection based on the NGE exceeded 0.49, and the accuracy of simulation increased with heritability and sample size. With a heritability of H²= 0.05, although the accuracy was only 0.49 when the sample size was 66, it significantly improved to 0.82 when the size was increased to 100. When the sample size was 200, the accuracy of the two heritability levels reached 1. The results show that the NGE achieves good accuracy in selecting QTLs for multiphasic growth and is suitable for large sample size data analysis.

On the other hand, in the absence of QTL expression, the same genetic sample sizes of n = 66 (the same as the real example), 100, and 200 and heritability levels of H²= 0.05 and 0.1 are simulated with 1,000 genetic markers and phenotypic data. At this point, the proportion of the simulation times of QTLs screened out to the number of 100 simulations is the false positive probability. The false positive rate (FPR) of QTL detection based on the NGE is small (< 0.11), indicating that the NGE has reasonably high specificity even under the conditions of large sample size and heritability (Table 2).

According to the QTL mapping accuracy and false positive probability in a series of different thresholds, we express ROC curves for different simulated sample sizes. The area under the ROC curve (AUC) was calculated to assess the accuracy of QTL mapping of the NGE (Figure 7). Under several different simulation conditions, the AUC value was less than 0.5 only when H²= 0.05 and n = 66, and other simulation results had high QTL detection significance. At the heritability level of 0.1, the AUCs of the simulated quantities of the three scales were all relatively high (>0.9659). At a certain heritability level, the accuracy of QTL mapping based on the NGE is affected by the number of genetic samples, indicating that the model is more suitable for QTL detection with large sample size.

We perform 100 simulation calculations based on the NGE parameters of the two genotypes of a significant QTL at the heritability levels of H² = 0.05 and 0.1, with sample sizes of n = 66, 100, and 200. The two genotypes are expressed as AA and Aa. The maximum likelihood estimation of NGE parameters (α_jH, K_jH1, β_jH←D, α_jD, K_jD1, β_jD←H, K_jH2, p_jH, q_jH, K_jD2, p_jD, q_jD) of simulated data is obtained by the fourth-order Runge–Kutta, EM, and Nelder–Mead algorithms. The parameter estimation results are shown in Table 3. Numerically, there are some differences between the estimated parameters and the real parameters. However, with the increase in genetic heritability and sample size, the precision of parameter estimation improves, the standard deviation of various parameters declines, and the simulation effec t becomes more stable.

Among the estimated parameters of 100 groups of simulation data obtained under the six kinds of simulation conditions, we select several groups of parameters with better effects to produce the overall growth curves and growth curves of two phases of two interaction traits. The trends of the estimated curves of different simulation scales are consistent with those of the real curves (Figures 8, 9), which suggests that the genotype-specific curves estimated from the above example are reasonably convincing. By comparing the estimated curves of the two heritability levels, the simulation effect of the estimated curve of 0.1 heritability is obviously better than that of the estimated curve of 0.05 heritability. For heritabilities of both levels 0.05 and 0.1, the green estimated curve with 200 simulated quantities is the closest to the actual curve, and the simulation effect is the best, followed by the blue curve with 100 simulated quantities, and the red curve with 66 simulated quantities. As expected, the simulation effects increase with sample size. Larger sample sizes, such as 200, can effectively minimize noise and are suitable for QTL mapping using a multiphasic growth model.