We explored the genetic structure of gravitropism and phototropism and its internal relations in a forest poplar tree. According to Zhang et al. [27], a full-sib F1 population of 345 members was derived from a cross between two dioecious P. euphratica trees in Korla, Xinjiang, China (41${}^{\circ }$10’48”–41${}^{\circ }$21’36”N, 85${}^{\circ }$14’10”–86${}^{\circ }$24’21”E). In the spring of 2014, two species of male and female flowers were selected and bred for artificial hybridization. After pollination for more than 4 months the seeds were harvested, and each seed was cultured in a sterile glass tube (400 mm in length and 40 mm in diameter). The tubes contained 350 mL 1/2 Murashige and Skoog medium (pH 6.0), and were laid out in a phytotron set at a 14-h day/10-h night cycle, 28 ${}^{\circ }$C day and 22 ${}^{\circ }$C night temperatures, with 800 $\mu$mol m${}^{-2}$ s${}^{-1}$ photosynthetically active radiation.

Among all seeds growing in culture, 345 normal-growth progeny were ultimately obtained. After planting in the glass tubes, the first observation time was individual-dependent and varied from 6 to 23 days. Main stem length and taproot length were measured repeatedly for each progeny once per week until the fastest-growing seedlings filled the tubes. A total of 12 or 16 measurements were taken over 3 months. Owing to the differences in the emerging stage of seedlings among progeny, and our efforts to reduce measurement errors by limiting the number of surveyors, the measurement schedules were slightly progeny-dependent.

Based on Zhang et al. [27], a total of 8279 single nucleotide polymorphisms (SNPs) were available, of which 6868 and 1411 were testcross and intercross markers, respectively. Testcross markers are those that segregate based on single-parent heterozygosity, whereas intercross markers segregate based on the heterozygosity of both parents.

From the perspective of the univariate model, the mechanisms of gravitropism and phototropism can be associated with QTLs through main stem height (S) and taproot length (R), respectively. As longitudinal traits, the growth of each trait can be fitted by a logistic function based on fundamental principles of biophysics and biochemistry [28]. This leads to the growth equation, expressed as:

(1)$$S(t)=\frac{a_S}{1+b_S{e}^{-r_{S}t}};$$

(2)$$R(t)=\frac{a_R}{1+b_R{e}^{-r_{R}t}};$$

where $(a_S,a_R)$ are the asymptotic values of the two component parts, main stem and taproot growth; $(b_S,b_R)$ denote the initial growth of the two component parts; and $(r_S,r_R)$ are the relative growth rates of the two component parts.

Main stem height and taproot length as longitudinal traits can be assumed to follow a normal distribution. Phenotypic values measured at a series of time points $(T_1,T_2,\mathrm{\dots },T_t)$ can be defined as vectors; e.g., $S=(S\left(T_1\right),\mathrm{\cdots },S\left(T_t\right)),R=(R\left(T_1\right),\mathrm{\cdots },R\left(T_t\right))$. Let ${\mu }_S$ and ${\mu }_R$ denote the mean vectors of these traits. The probability density functions of main stem height and taproot length for a particular progeny $i$ can be expressed as $f\left(S_i\right)$ and $f\left(R_i\right)$, respectively.

(3)$$\mathrm{f}\left({\mathrm{S}}_{\mathrm{i}}\right)=\frac{1}{\sqrt{\left|2\pi {\mathrm{\Sigma }}_{{\mathrm{S}}_{\mathrm{i}}}\right|}}\exp \left(\frac{-1}{2}{\left({\mathrm{S}}_{\mathrm{i}}-{\mu }_{{\mathrm{S}}_{\mathrm{i}}}\right)}^{\mathrm{T}}{\mathrm{\Sigma }}_{{\mathrm{S}}_{\mathrm{i}}}^{-1}\left({\mathrm{S}}_{\mathrm{i}}-{\mu }_{{\mathrm{S}}_{\mathrm{i}}}\right)\right)$$

(4)$$\mathrm{f}\left({\mathrm{R}}_{\mathrm{i}}\right)=\frac{1}{\sqrt{\left|2\pi {\mathrm{\Sigma }}_{{\mathrm{R}}_{\mathrm{i}}}\right|}}\exp \left(\frac{-1}{2}{\left({\mathrm{R}}_{\mathrm{i}}-{\mu }_{{\mathrm{R}}_{\mathrm{i}}}\right)}^{\mathrm{T}}{\mathrm{\Sigma }}_{{\mathrm{R}}_{\mathrm{i}}}^{-1}\left({\mathrm{R}}_{\mathrm{i}}-{\mu }_{{\mathrm{R}}_{\mathrm{i}}}\right)\right)$$

Where $S_i=(S_i\left(T_{i1}\right),\mathrm{\cdots },S_i\left(T_{it}\right))$ and $R_i=(R_i\left(T_{i1}\right),\mathrm{\cdots },R_i\left(T_{it}\right))$ are the longitudinal phenotypic vectors of two traits measured over t time points for progeny $i$, ${\mu }_{{\mathrm{S}}_i}=({\mu }_{{\mathrm{S}}_i}\left(T_{i1}\right),\mathrm{\cdots },{\mu }_{{\mathrm{S}}_i}\left(T_{it}\right))$ and ${\mu }_{{\mathrm{R}}_i}=({\mu }_{{\mathrm{R}}_i}\left(T_{i1}\right),\mathrm{\cdots },{\mu }_{{\mathrm{R}}_i}\left(T_{it}\right))$ are the expected means of the corresponding traits over the time points for progeny $i$. $(t\times t)$ covariance matrix modeled by SAD(1):

(5)$\stackrel{S_i}{\overbrace{{\sigma }_{T_{i1}{T}_{i1}}\mathrm{\cdots }{\sigma }_{T_{i1}{T}_{it}}}}$
${\mathrm{\Sigma }}_{S_i}=\mathrm{⋮}\mathit{\hspace{1em}\hspace{1em}}\mathrm{\ddots }\mathit{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{⋮}$
${\sigma }_{T_{it}{T}_{i1}}\mathrm{\cdots }{\sigma }_{T_{it}{T}_{it}}$

(6)$\stackrel{R_i}{\overbrace{{\sigma }_{T_{i1}{T}_{i1}}\mathrm{\cdots }{\sigma }_{T_{i1}{T}_{it}}}}$
${\mathrm{\Sigma }}_{R_i}=\mathrm{⋮}\mathit{\hspace{1em}\hspace{1em}}\mathrm{\ddots }\mathit{\hspace{1em}\hspace{1em}\hspace{1em}}\mathrm{⋮}$
${\sigma }_{T_{i1}{T}_{it}}\mathrm{\cdots }{\sigma }_{T_{it}{T}_{it}}$

Considering that a multigeneration of trees is not easy to achieve, and the genotypes of trees cannot be completely pure, we carried out QTL mapping within the F1 generation. We assumed that it is composed of n progeny, each genotyped at the genome level for molecular markers from which a genome-wide genetic linkage map can be constructed.

Because the specific measurement time for each individual was different, the likelihood of main stem height and taproot length at a marker locus of J genotypes is formulated as

(7)$$L(S)=\prod _{j=1}^J\prod _{i=1}^{n_j}{f}_j\left(S_i\mid {\theta }_j^S,{\mathrm{\Sigma }}_{S_i}\right)$$

(8)$$L(R)=\prod _{j=1}^J\prod _{i=1}^{n_j}{f}_j\left(R_i\mid {\theta }_j^R,{\mathrm{\Sigma }}_{R_i}\right)$$

where J is the number of genotype at a marker locus. $n_j$ is the number of individual of genotype j. ${\theta }_j^s=(a_{sj},b_{sj},r_{sj})$ is the parameter set of genotype j for main stem. ${\mathrm{\Sigma }}_{s_i}$ is the covariance matrix which referred from the formula Eqn. 5. ${\theta }_j^R=(a_{Rj},b_{Rj},r_{Rj})$ is the parameter set of genotype j for taproot. ${\mathrm{\Sigma }}_{R_i}$ is the covariance matrix which referred from the formula Eqn. 6.

The bivariate model is capable of locating significant QTL sites that affect both phototropism and gravitropism based on phenotypic variation in main stem height and taproot length at a given time point. Phenotypic values can be defined as vectors, i.e., ${\mu }_{\mathrm{SR}}=((S\left(T_1\right),\mathrm{\cdots }S\left(T_t\right)),(R\left(T_1\right),\mathrm{\cdots },R\left(T_t\right)))$, and ${\mathrm{\Sigma }}_{SR}$, the $(2*t\times 2*t)$ covariance matrix, is modeled by SAD(2) [29].

(9)${\mathrm{\Sigma }}_{SR}=\left(\stackrel{𝑆}{\overbrace{\begin{array}{c}\hfill \begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{i1}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{i1}}\hfill \end{array}\begin{array}{c}\hfill \mathrm{\cdots }\hfill \\ \hfill \mathrm{\ddots }\hfill \\ \hfill \mathrm{\cdots }\hfill \end{array}\begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{it}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{it}}\hfill \end{array}\hfill \\ \hfill \begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{i1}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{i1}}\hfill \end{array}\begin{array}{c}\hfill \mathrm{\cdots }\hfill \\ \hfill \mathrm{\ddots }\hfill \\ \hfill \mathrm{\cdots }\hfill \end{array}\begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{it}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{it}}\hfill \end{array}\hfill \end{array}}}\stackrel{𝑅}{\overbrace{\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{i1}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{i1}}\hfill \end{array}\begin{array}{c}\hfill \mathrm{\cdots }\hfill \\ \hfill \mathrm{\ddots }\hfill \\ \hfill \mathrm{\cdots }\hfill \end{array}\begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{it}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{it}}\hfill \end{array}\hfill \\ \hfill \begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{i1}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{i1}}\hfill \end{array}\begin{array}{c}\hfill \mathrm{\cdots }\hfill \\ \hfill \mathrm{\ddots }\hfill \\ \hfill \mathrm{\cdots }\hfill \end{array}\begin{array}{c}\hfill {\sigma }_{T_{i1}{T}_{it}}\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\sigma }_{T_{it}{T}_{it}}\hfill \end{array}\hfill \end{array}\hfill \end{array}}}\right)$

We make a few simplifications of the structure of ${\mathrm{\Sigma }}_{SR}$. First, we assume that over all the T time points, the ${\gamma }_S$ and ${\gamma }_R$, which are termed time-dependent “innovation variances” and described by a parametric function like as a polynomial of time [30], are constant. With this assumption, the residual variance for trait k (main stem or taproot) at time t is given by

(10)$$\mathrm{Var}\left(T_i\right)=\frac{1-{\mathrm{\Phi }}_k^{2T_i}}{1-{\mathrm{\Phi }}_k^2}{\gamma }_k^2,T_i=T_{i1},T_{i2},\mathrm{\dots },T_{it}$$

where ${\mathrm{\Phi }}_k$ is the antedependence parameters induced by trait k itself. For $T_{i}2\ge T_{i}1$, the correlation function of the SAD(2) model is written as

(11)$$\mathrm{Corr}(T_{i1},T_{i2})={\mathrm{\Phi }}_k^{T_{i2}-T_{i1}}\sqrt{\frac{1-{\mathrm{\Phi }}_k^{2T_{i1}}}{1-{\mathrm{\Phi }}_k^{2T_{i2}}}},T_{i2}>T_{i1}$$

According to Eqn. 10, even for constant innovation variances, the residual variance can change over time [31]. It also can be seen from Eqn. 11 that the correlation function is unsteady for the SAD model because the time intervals are not the only determinants of the correlation.

Second, we assume that the innovation correlation exists only between the two component traits at the same time point, not at different time points, and that such an innovation (i.e., $\rho$) is constant over different times. Third, we assume that the antedependence parameters of the correlated traits are equal to zero [31]. Thus, the residual correlation caused by the two traits at different times ($t_1$ for main stem and$t_2$ for taproot) is derived as

(12)
$\mathrm{Corr}({\mathrm{T}}_{\mathrm{i1}},{\mathrm{T}}_{\mathrm{i2}})={\displaystyle \frac{{\mathrm{\Phi }}_{{\mathrm{k}}_2}^{{\mathrm{T}}_{\mathrm{i2}}-{\mathrm{T}}_{\mathrm{i1}}-{\mathrm{\Phi }}_{{\mathrm{k}}_1}^{{\mathrm{T}}_{\mathrm{i1}}}{\mathrm{\Phi }}_{{\mathrm{k}}_2}^{{\mathrm{T}}_{\mathrm{i2}}}}}{1-{\mathrm{\Phi }}_{{\mathrm{k}}_1}{\mathrm{\Phi }}_{{\mathrm{k}}_2}}}{\left(\sqrt{{\displaystyle \frac{1-{\mathrm{\Phi }}_{{\mathrm{k}}_1}^{2{\mathrm{T}}_{\mathrm{i1}}}}{1-{\mathrm{\Phi }}_{{\mathrm{k}}_1}^2}}{\displaystyle \frac{1-{\mathrm{\Phi }}_{{\mathrm{k}}_2}^{2{\mathrm{t}}_{{\mathrm{T}}_{\mathrm{i}}2}}}{1-{\mathrm{\Phi }}_{{\mathrm{k}}_2}^2}}}\right)}^{-1}$
$\rho ,{\mathrm{T}}_{\mathrm{i2}}>{\mathrm{T}}_{\mathrm{i1}},{\mathrm{k}}_1\ne {\mathrm{k}}_2$

where ${\mathrm{\Phi }}_{k_1}$ and ${\mathrm{\Phi }}_{k_2}$ are the antedependence parameters induced by trait $k_1$ itself and $k_2$ itslf, respectively.

According to Eqn. 12, the across-correlation functions are not symmetrical for the SAD model. That is, the correlation between the main stem trait at time $T_{i1}$ and the taproot trait at time $T_{i2}$ is not equal to the correlation between main stem trait at time $T_{i1}$ and taproot trait at time $T_{i1}$.

Thus, by implementing these expressions in Eqn. 10,11,12 into the matrix in Eqn. 9, we can achieve a parsimonious matrix structure without losing much information.

Under the bivariate functional model, given that the specific measurement time for each individual is different, the likelihood of the two phenotypes at a marker locus of J genotypes is formulated as

(13)$$L(S,R)=\prod _{j=1}^J\prod _{i=1}^{n_j}{f}_j(S_i,R_i\mid {\theta }_j^{SR},{\mathrm{\Sigma }}_{S{R}_i})$$

where J is the number of genotype at a marker locus. $n_j$ is the number of individual of genotype j, ${\theta }_j^{SR}=(a_{sj},b_{sj},r_{sj},a_{Rj},b_{Rj},r_{Rj})$ is the parameter set of genotype j. ${\mathrm{\Sigma }}_{S{R}_i}$ is the covariance matrix which referred from the formula Eqn. 9.

Through the anatomical analysis of the plant’s main axis, a composite trait, we can determine how the significant loci affect the growth of the entire whole plant through phototropism and gravitropism. The relationship between the main axis and the two component traits is shown in Fig. 1. We partition the overall growth of the main axis, M, into two parts, expressed as:

$M=S+R$

Main stem height and taproot length can be studied allometrically. According to the allometric growth equation, the main stem equation is represented by the taproot formula. We make the summation of logistic functions fit multiple parts (main stem and taproot) of the main axis, leading to a composite growth equation, expressed as:

(14)$$M(t)=S(t)+R(t)=\alpha {\left(\frac{a_R}{1+b_R{e}^{-r_{R}t}}\right)}^{\beta }+\frac{a_R}{1+b_R{e}^{-r_{R}t}}$$

where $a_R$ is the asymptotic value of taproot length, $b_R$ denotes the initial growth of the taproot, $r_R$ is the relative growth rate of the taproot, and $\alpha$ and $\beta$ are the coefficients of root and stem relationships, respectively.

The main axis as a longitudinal trait is actually the sum of two variables, which are also longitudinal traits can be assumed to follow a normal distribution. Phenotypic values can be defined as vectors: $M=(M\left(T_1\right),\mathrm{\cdots },M\left(T_t\right))=(S\left(T_1\right)+R\left(T_1\right),\mathrm{\cdots },S\left(T_t\right)+R\left(T_t\right))=\left((S\left(T_1\right),\mathrm{\cdots },S\left(T_t\right))+(R\left(T_1\right),\mathrm{\cdots },R\left(T_t\right))\right)=S+R$. Let ${\mu }_S$, ${\mu }_R$, and ${\mu }_M$ denote the mean vectors of these traits so that we have ${\mu }_M={\mu }_S+{\mu }_R$. The probability density function of the main axis for a particular progeny i can be expressed as $f\left(M_i\right)$:

(15)$\mathrm{f}\left({\mathrm{M}}_{\mathrm{i}}\right)={\displaystyle \frac{1}{\sqrt{\left|2\pi {\mathrm{\Sigma }}_{{\mathrm{M}}_{\mathrm{i}}}\right|}}}\exp \left({\displaystyle \frac{-1}{2}}{\left({\mathrm{M}}_{\mathrm{i}}-{\mu }_{{\mathrm{M}}_{\mathrm{i}}}\right)}^{\mathrm{T}}{\displaystyle \sum _{{\mathrm{M}}_{\mathrm{i}}}^{-1}}\left({\mathrm{M}}_{\mathrm{i}}-{\mu }_{{\mathrm{M}}_{\mathrm{i}}}\right)\right)$
$={\displaystyle \frac{1}{\sqrt{\left|2\pi {\mathrm{\Sigma }}_{{\mathrm{M}}_{\mathrm{i}}}\right|}}}\exp \left({\displaystyle \frac{-1}{2}}{\left[\left({\mathrm{S}}_{\mathrm{i}}+{\mathrm{R}}_{\mathrm{i}}\right)-\left({\mu }_{{\mathrm{S}}_{\mathrm{i}}}+{\mu }_{{\mathrm{R}}_{\mathrm{i}}}\right)\right]}^{\mathrm{T}}{\displaystyle \sum _{{\mathrm{M}}_{\mathrm{i}}}^{-1}}\left[\left({\mathrm{S}}_{\mathrm{i}}+{\mathrm{R}}_{\mathrm{i}}\right)-\left({\mu }_{{\mathrm{S}}_{\mathrm{i}}}+{\mu }_{{\mathrm{R}}_{\mathrm{i}}}\right)\right]\right)$
$={\displaystyle \frac{1}{\sqrt{\left|2\pi {\mathrm{\Sigma }}_{{\mathrm{M}}_{\mathrm{i}}}\right|}}}\exp \left({\displaystyle \frac{-1}{2}}{\left[\left({\mathrm{S}}_{\mathrm{i}}-{\mu }_{{\mathrm{S}}_{\mathrm{i}}}\right)+\left({\mathrm{R}}_{\mathrm{i}}-{\mu }_{{\mathrm{R}}_{\mathrm{i}}}\right)\right]}^{\mathrm{T}}{\displaystyle \sum _{{\mathrm{M}}_{\mathrm{i}}}^{-1}}\left[\left({\mathrm{S}}_{\mathrm{i}}-{\mu }_{{\mathrm{S}}_{\mathrm{i}}}\right)+\left({\mathrm{R}}_{\mathrm{i}}-{\mu }_{{\mathrm{R}}_{\mathrm{i}}}\right)\right]\right)$

$M_i=({}_{-}M_i\left(T_{i1}\right),\mathrm{\cdots },M_i\left(T_{it}\right))$ are the longitudinal phenotypic vectors of the main axis measured over t time points, ${\mu }_{M_i}=({\mu }_{M_i}\left(T_{i1}\right),\mathrm{\cdots },{\mu }_{M_i}\left(T_{it}\right))$ are the expected means of the main axis over the time points, and ${\mathrm{\Sigma }}_{M_i}$ is the covariance matrix, which is modeled by SAD(1).

Under the composite trait model, given that the specific measurement time for each individual is different, the likelihood of the main axis at a marker locus of J genotypes is formulated as

(16)$$L(M)=\prod _{j=1}^J\prod _{i=1}^{n_j}{f}_j\left(M_i\mid {\theta }_j^M,{\mathrm{\Sigma }}_{M_i}\right)$$

where J is the number of genotype at a marker locus. $n_j$ is the number of individual of genotype j, ${\theta }_j^M=({\alpha }_j,{\beta }_j,a_{Rj},b_{Rj},r_{Rj})$ is the parameter set of genotype j. ${\mathrm{\Sigma }}_{M_i}$ is the covariance matrix, which is modeled by SAD(1).

Many algorithms have been developed for parameter estimation (Eqns. 8,16). The simplex algorithm, which is derived from the optimization principle, has been widely applied to estimate growth parameters and matrix-structuring parameters [32]. In this study, we deployed the simplex algorithm for the estimation of all parameters contained in the equations.

By plotting main stem length and taproot length data over time for the progeny, we identified a logistic equation that fit the growth curve (Fig. 2). According to the univariate functional mapping model, Eqns. 1,2 can be directly incorporated into the likelihood model to map phototropism- or gravitropism-related QTLs. Eleven phototropism-related markers and eight gravitropism-related markers were identified as QTLs within, or adjacent to, candidate genes that have been observed to function in known plant growth-related pathways (Fig. 3). The significant sites related to phototropism were distributed mainly in multiple linkage groups, such as 2, 5, 10, 15, 17, and 18. Conversely, the significant sites related to gravitropism were distributed mainly in three linkage groups: 8, 10, and 12.

The phototropism-associated QTL P. euphratica cytochrome P450 CYP72A219 (CYP72A219) is a tolerant gene that is continuously expressed at high levels in nearly all salt-tolerant crops. The gravitropism-related QTL ABC transporter B family member 8 (ABCB8) is related to ATPase activity and ATP-mediated material transport. We further characterized the genotypes and dominant genetic effects of the phototropism-related testcross QTL CYP72A219 and the gravitropism-related testcross QTL ABCB8 for each trait (Fig. 4).

QTL CYP72A219 has a clear impact on phototropic growth, gradually increasing the growth of the main stem over time. Seedling tolerance to stress increases gradually, but at an increasing rate, within the first three months of growth. QTL ABCB8, which influences geotropic growth, has increasing effects throughout the early stages of seedling germination; however, its effects begin to decline in the third month of seedling germination, the key period for seedling establishment. Gravitropism is clearly a very important process during this period. Although this method is powerful for detecting QTLs influencing either phototropism or gravitropism, the model does not take into account the biological environment of the main stem and taproot or the genetic relationship between phototropism and gravitropism.

Using the bivariate functional mapping model, eight markers were identified as QTLs (Fig. 5) significantly related to phototropism, gravitropism, or both. We chose the five most significant QTLs, all of which were annotatable, for in-depth analysis. SNP lm_ll_12166 resides in the QWRF motif, which is related to microtubule-severing activity, which in turn plays an important role in intercellular spreading. SNP lm_ll_8106 was observed to reside within the F-box gene, which directs the interaction of the protein complex with a specific target for ubiquitylation. SNP lm_ll_12477 was located within the polyubiquitin (UBQ) gene, which plays distinct roles that include signaling and the inactivation of protein kinases. SNP lm_ll_12839 was found within the ethylene-responsive transcription factor 12-like (TAF12) gene, which is essential for mediating the regulation of RNA polymerase transcription and is likely involved in the negative regulation of cytokinin sensitivity. SNP hk_hk_1862 resides within the transcription factor bHLH95-like (BHLH95) gene, which encodes a transcription factor that controls embryonic growth. Polyadenylation and cleavage factor homolog 4-like (PCFS4) is involved in plant leaf shape changes and heteromorphic leaf variations, which may be related to plant environmental adaptability. Other significant SNPs can be similarly annotated (Table 1).

We further modeled the separation curves of genotypes and the dominant genetic effects of testcross QTLs glutamate receptor 2.1-like (GLUR) and TAF12, mapped using the bivariate model, on each trait (Fig. 6, Supplementary Fig. 1). GLUR attenuates wound-induced surface potential changes and is involved in organ-organ long-distance trauma signaling. The first three months of seedling germination is a critical period for root formation. The effects of GLUR on root development initially increase in strength, then decrease. Roots can respond rapidly to their environment by transmitting internal and external environmental signals and adjusting their morphology. At the end of the first three months of germination, phototropic growth gradually increases, and GLUR begins to be actively expressed in the stem. TAF12 represents a group of plant AP2/ERF transcription factors, including the pathogenesis-related transcriptional activator PTI1-like tyrosine-protein kinase 6 (PTI6) and the ethylene-responsive transcription factor ERF. PTI6 binds to the GCC-box pathogenesis-related promoter element and activates plant defense genes. It is possible that this site is involved in the regulation of defense processes in the stem and promotes the photogenic growth of the stem. However, these effects do not explain the physiological components of main axis formation, and may therefore be limited in revealing the genetic machinery underlying how phototropism and gravitropism influence main axis growth.

We drew the dynamic trajectories of the main stem (Fig. 7A), main axis (Fig. 7B), and taproot (Fig. 7C) based on the maximum likelihood estimates for the growth parameters of the two component traits ($\alpha ,\beta$,$a_R$ ,$b_R$ ,$r_R$).

By incorporating the composite growth equation (Eqn. 14) into the likelihood estimation (Eqn. 16), the composite trait model unified the physiological components of trait formation and identified seven QTLs (Fig. 5). Two of the seven significant sites were also identified using the bivariate functional model (Table 2). We further estimated the separation curves of genotypes and the dominant genetic effects of the two testcross QTLs, PCFS4 and UBQ, on each trait (Fig. 8).

In the early stages of seedling growth, the boundary between root and stem is unclear, and root growth is critical to the successful establishment of poplar seedlings. In later stages of growth, the boundary between root and stem is obvious. PCFS4 has important effects on stem development. First, the main axis grows downward through root development. A stable root system is important for the growth of the stem. Photogenic growth causes the stem to grow upward on the main axis so that the plant can increase its exposure to sunlight and absorb more carbon. UBQ accepts ubiquitin from the E1 complex and catalyzes its covalent attachment to other proteins. UBQ can promote the bidirectional growth of the principal axis by inducing both stem and root growth.