- 1College of Science, Beijing Forestry University, Beijing, China
- 2Department of Artificial Intelligence and Data Science, Guangzhou Xinhua University, Guangzhou, China
- 3College of Biological Sciences and Technology, Center for Computational Biology, Beijing Forestry University, Beijing, China
- 4Faculty of Science, Yamagata University, Yamagata, Japan
- 5Yau Mathematical Sciences Center, Tsinghua University, Beijing, China
Introduction: The cooperative strategy of phenotypic traits during the growth of plants reflects how plants allocate photosynthesis products, which is the most favorable decision for them to optimize growth, survival, and reproduction response to changing environment. Up to now, we still know little about why plants make such decision from the perspective of biological genetic mechanisms.
Methods: In this study, we construct an analytical mapping framework to explore the genetic mechanism regulating the interaction of two complex traits. The framework describes the dynamic growth of two traits and their interaction as Differential Interaction Regulatory Equations (DIRE), then DIRE is embedded into QTL mapping model to identify the key quantitative trait loci (QTLs) that regulate this interaction and clarify the genetic effect, genetic contribution and genetic network structure of these key QTLs. Computer simulation experiment proves the reliability and practicability of our framework.
Results: In order to verify that our framework is universal and flexible, we applied it to two sets of data from Populus euphratica, namely, aboveground stem length - underground taproot length, underground root number - underground root length, which represent relationships of phenotypic traits in two spatial dimensions of plant architecture. The analytical result shows that our model is well applicable to datasets of two dimensions.
Discussion: Our model helps to better illustrate the cooperation-competition patterns between phenotypic traits, and understand the decisions that plants make in a specific environment that are most conducive to their growth from the genetic perspective.
1 Introduction
In the growth and development of plants, there are always intimate communication and connection among organs or various parts of the same organ, so that plants can maintain the optimum condition to obtain natural resources and adapt phenotypically to living environment (Niklas and Spatz, 2006; Poorter et al., 2009; Freschet et al., 2015; Weraduwage et al., 2015; Freschet et al., 2018; Tumber-Dávila et al., 2022). The communication and connection are specifically manifested as allocation patterns of biomass (Hermans et al., 2006; Poorter et al., 2012; Wang et al., 2022), which forms different allometric relationships and shapes plant morphology, such as the restrictive relation between aboveground stem and underground roots, height and diameter of stem, etc. (Niklas and Spatz, 2006; Fu et al., 2017). Besides external environmental factors, heredity is a critical driving force for individuals to show a variety of phenotypic characteristics. (Cunniff et al., 2015; Magar et al., 2021).
Quantitative trait locus (QTL) mapping can locate QTLs or genes related to complex quantitative traits on chromosome, which is always the interest of biogenetics and has been successfully applied to many breeding projects (Tang et al., 2010; Wei et al., 2012; Sonah et al., 2015; Fu et al., 2016; Xia et al., 2018). However, traditional QTL mapping mainly focuses on phenotypic data of one trait at one time point, ignoring that plant growth is a dynamic and continuous process. Functional mapping is a QTL mapping method to identify QTLs that regulate the dynamic growth change of target trait during a period of time (Wu and Lin, 2006; Sun and Wu, 2015; Fu et al., 2017).
In this paper, we establish a general and universal model framework, which takes into account the dynamic characteristics of plant growth and the interaction between two traits, screens out key QTLs regulating this pattern, and figures out the genetic structure (Figure 1). Our model is based on the Differential Interaction Regulatory Equations (DIRE) to describe the growth of two closely related phenotypic traits and the potential interaction pattern between them. The basic form of this equation comes from the Lotka-Volterra (LV) equation, which described the ecological interaction of two species (Lotka, 1920; Volterra, 1928; May, 1973). Our model identifies a number of genetic loci that determine the cooperation or competition pattern of target traits. According to the genetic effects of these key loci, the genetic network among them is further constructed. We applied the model on two sets of data, stem length - taproot length and root number - root length, which represent the vertical and horizontal relationships of plant architecture, respectively. The balance of competition or cooperation between stem length and taproot length should not only meet the requirements of plants to absorb, synthesize and transport nutrients needed for life, but also ensure plant morphology and stability. Specifically, the extension of root system under the ground makes plants absorb water, minerals elements and other substances in the soil, while stem transports carbohydrate to root and provides support for leaves to ensure effective photosynthesis. And root system shall provide enough resistance to prevent the plant from being pulled out (Modrzyński et al., 2015). The relation between root number and root length reflects the morphological characteristics of root system. Such comparison enables us to have a new understanding of the coordinated variation between traits from two dimensions of plant structure and the genetic mechanism. At the same time, it provides a reference for plant genetic breeding according to the cooperative relationship of traits by means of technical means. Finally, the simulation experiment verifies the reliability of our model.
Figure 1 Schematic diagram of our model framework. It is applied to the data of stem length and taproot length from Populus euphratica seedlings to identify key QTLs regulating the growth interaction between aboveground and underground, and clarify the architecture of genetic network. It is also used to analyze the data of root number and root length to understand the genetic mechanism.
2 Materials and methods
2.1 Plant materials
We used the published experimental data to validate the utility of our model (Zhang et al., 2017; Zhang, 2019). A full-sib F1 family of 321 members was derived from hybridization between two dioecious trees of Populus euphratica (Zhang et al., 2017). In spring 2014, cut the male and female flowering branches from the two trees and conduct artificial hybridization in water. After more than 4 months, the catkins gradually ripened, and the harvested seeds were cultured in a glass tube (40 mm in diameter and 400 mm in length) which contains 350 ml of 1/2 Murashige and Skoog medium (pH 6.0) under a sterile condition. The tube was placed in a phytotron set at 14-h-day/10-h-night cycle, 28°C day and 22°C night with 800 photosynthetically active radiation. Two phenotypes, namely stem height (mm) and taproot length (mm) were measured repeatedly for each progeny since seed germination. Measurements were undertaken 16 times: 1, 3, 5, 7, 18, 20, 22, 24 26, 28, 31, 34, 38, 47, 54, 62 days after seed germination, respectively. The full-sib population with the same dioecious parents was planted (Zhang, 2019). In 2014, the flowering branches of the male parent were cultured in water in a phytotron set, and the mature pollen was stored in the EP tube at -20°. The whole plant of the female parent was planted in the greenhouse. After the female flower matures, the male parent's pollen is used for artificial pollination. The mature seeds were obtained in the middle of June. After the seeds were planted in vitro for 4 months, the seedlings were transplanted to the substrate for cultivation, and 408 full-sib offspring were finally obtained. Through the clonal experiment and the expansion and preservation experiment of 408 individuals in the population, 156 groups that can take root on the rooting medium and grow normally into complete plants were finally obtained. Vitro inoculation experiments on 156 groups were conducted. Intercept the terminal buds (about 10 mm, 4-6 leaves) with the same growth status from each single plant and inoculate them into a cylindrical flat-bottomed glass tube (45 mm in diameter and 300 mm in length) for culture. All tubes were added with 260 ml of rooting medium, and placed under the uniform culture conditions in the tissue culture room at 16-h-day/8-h-night cycle, 22°C with 1500lx light. Root length (cm) and root number (count) were collected for each progeny once every 5 days from the 13th day to the 78th day after cultivation. Through the high-throughput sequencing, we obtained 8305 single-nucleotide polymorphisms (SNPs) distributed throughout 19 linkage groups (labeled as Q1-Q8035). All SNP markers can be divided into testcross markers (lm×ll and nn×np) and intercross hybrid markers (hk×hk) according to Mendelian genetic segregation rules, with 6886 and 1419 markers, respectively.
2.2 Differential interaction regulatory equations
The biomass allocation among phenotypic traits results in diverse allometric growth relationship. The relationship is a trade-off made by trees after they adapt to the external environment, which largely depends on their internal biological mechanisms. The multiphasic growth model shows that the growth of most organisms is composed of multiple “S”-shaped phases, which is superior to single-phase growth (Koops, 1989; Van der Klein et al., 2020; Gong et al., 2021). The seedling stage of trees is a “S”-shaped phase. According to the basic principles of biophysics and biochemistry, Logistic growth curve shows an "S" shape and consists of exponential growth, linear growth and asymptotic growth (West et al., 2001). Lotka-Volterra equation is derived from the Logistic growth curve and is originally used as dynamic model for growth and decline of amount on predators and prey in biological systems. Its feasibility and effectiveness in describing the microbial abundance under coculture and the predator-prey relationship in the ecosystem have been fully proved (Fujikawa et al., 2014; Jiang et al., 2018). We propose a DIRE model to quantify allometric growth relationship of different traits, and its basic form comes from Lotka-Volterra equation. DIRE describes the cooperation or competition between two traits in the growth process to make full use of survival resources. DIRE model is expressed as
DIRE consists of two parts:
representing the independent growth of target character, and
representing interactive growth of two traits, where and represent phenotypic values of two characters, and are the asymptotic value of independent growth, and represent the independent growth rate, scale parameters and control the independent growth rate, and is the interaction coefficient. The values of and reflect interaction types between traits. We can use strategic coordinate system to summarize all interaction types. According to the strategy coordinate system (Figure 2C), the strategies for interaction of two traits are divided into three categories: (1) Win-win cooperation: , two traits promote the growth of each other; (2) Predation-prey: , one character is conducive to the other while the latter is harmful to the former. (3) Internecine: , both traits are inhibited by the other.
Figure 2 Analysis of the growth of seedlings from a full-sib family of Populus euphratica. (A) Violin plots of stem length, taproot length, root number and root length of all samples at a series of times. (B) Growth curves of stem length, taproot length, root number and root length. Each red line is the fitting curve of average growth for corresponding trait, which consists of independent growth (broke lines) and interactive growth with another trait (dot lines). (C) Strategic coordinate system of interaction types between two traits. According to the sign of , the strategy plane is divided into four regions, representing four interaction types between traits. (D) CV curves of stem length, taproot length, root number and root length of all samples at a series of times.
2.3 Statistical model for identifying key QTLs
QTL mapping based on DIRE identifies QTLs that mediate the dynamic growth and interaction between traits (Wu et al., 2011; Sun and Wu, 2015). We design a full-sib mapping population of n seedlings. All samples are genotyped for p SNPs. Phenotypic values for trait 1 and trait 2 are obtained at a series of time points (t=1, …, T), the phenotypic value of seedling i is
where obeys bivariate normal distribution with mean vector and covariance matrix , i.e.
where the length of is 2T, and is a symmetric matrix. In other words, and represent the characteristics of the phenotype data for target population. is obtained by DIRE, corresponding to parameter set . can also be estimated by a set of parameters using first-order structured antedependence (SAD (1)) model (Jaffrézic et al., 2003; Zhao et al., 2005a; Zhao et al., 2005b).
Based on the following hypotheses, we can detect whether a QTL regulates the interaction growth:
where and . The null hypotheses means that parameters are independent of genotypes, and the alternative hypotheses indicates that there are genotype differences in the parameter sets.
Hypotheses test is realized by maximum likelihood estimation (MLE). The likelihood function of n samples for the null hypotheses is
where is a bivariate probability density function. And the likelihood function for the alternative hypotheses is expressed as
where J is the number of genotypes of target QTL, is the number of samples with genotype j , which satisfies . The likelihood values and are calculated, respectively, and the test statistic log-likelihood ratio (LR) is calculated as
LR follows the law of chi-square distribution according to its construction principle, and its degree of freedom is the difference between the number of parameters of and . p value is compared with the critical threshold after FDR (False Discovery Rate) correction to reduce false positives.
In the process of hypotheses test solution, parameter estimation is realized by Expectation Maximization (EM) algorithm (Moon, 1996; Do and Batzoglou, 2008). According to EM algorithm, we give an initialization parameter, and estimate the likelihood function value, then iterate and optimize it so that the likelihood function value approximates the local optimal value to obtain new parameters. When the increased value of the likelihood function is less than the target threshold, the iteration ends. We can obtain parameter sets and of each locus within certain precision, in which the solution of DIRE depends on the fourth-order Runge-Kutta algorithm (Blum, 1962).
2.4 Construction of genetic network
The significant QTLs obtained from hypotheses test play key roles in regulating the growth of phenotypic traits. By constructing the genetic network of these QTLs, we can better understand the genetic mechanism of phenotypic variation. We introduce the replication factor equation (Wu and Jiang, 2021), which regards the observed value of the target variable as the sum of its own strategy and interactive strategies with its counterparts. The two parts are referred to as the independent part and the interactive part, respectively, and they are expressed in the form of a nonlinear differential equation system. If these significant QTLs are regarded as a system, the overall genetic effect of each QTL in this system can be decomposed into independent part and the interactive part affected by other QTLs.
、 denote the overall genetic effects of QTL k (k=1, 2, …, K) for trait 1 and trait 2, respectively. We only retain QTLs that have a greater impact on QTL k through variable selection in the LASSO regression model:
where represents the overall genetic effect of QTL k, is the regression coefficient of QTL , is a constant and represents residual. LASSO (Wang and Leng, 2008) selects the most important set of QTLs for the target QTL, and the connections among them are determined by means of a nonlinear differential equation:
where the overall genetic effect is decomposed into two parts: describes the independent genetic effect of QTL k; represents interaction produced by other QTLs on QTL k, and are parameter sets, respectively. The fourth-order Runge-Kutta algorithm is required to solve the differential equation to determine the extent of interaction among QTLs.
On the other hand, the number of nodes in the network may be more than that of time points. To ensure the accuracy, Legendre orthogonal polynomial (LOP) is firstly applied to fit the genetic effect curve of this QTL and obtain more time point information by interpolation on the curve.
2.5 Coefficient of variation
Coefficient of Variation (CV) is a normalized measure to evaluate the dispersion of probability distribution. It can be used to compare the dispersion of data among different phenotypic traits on the basis of eliminating influence of dimension. We use it to quantitatively describe the change of dispersion of each phenotypic trait over time. The calculation formula is as follows:
where represents standard deviation and represents mean.
2.6 Phenotypic variation explained
We define the proportion of phenotypic variation caused by genetic factors (or heritability explained by chosen QTLs) as phenotypic variation explained (PVE), and it is calculated by means of variance analysis. First of all, take the target trait phenotype as the dependent variable and genotypes as the independent variables to conduct generalized linear regression, and then perform variance analysis on the regression result to obtain the PVE value (Xu et al., 2016; Gong et al., 2021).
3 Results
3.1 Phenotypic variation analysis
Phenotypic traits are expected to be variable among individuals within a species, which is one kind of manifestation of intraspecific variability (Cianciaruso et al., 2009; Albert et al., 2010) and may be caused by genetic factors or phenotypic plasticity induced by environment (Miner et al., 2005). Here, specific performance of this variation is discrepancy in quantitative characters among progeny seedlings, which may not be clear at the seed germination stage, but morphological variation and structural differences within the population become more and more obvious over time. Figure 2A illustrates the distribution of stem length - taproot length and root number – root length from F1 population of Populus euphratica at a series of time points. Phenotypic values of them both display bigger dispersion with time. The width of the violin gradually narrows, while the vertical height constantly increases, indicating that the differences among progenies increase over time. Taking taproot length as an example for detailed analysis, the height range of the violin at 3 days is 0-20mm, that is, taproot length of all samples does not exceed 20mm at this time. The widest width of the violin shows the concentrated distribution of taproot length in F1 population, and the taproot length of most seedlings is between 1 and 5mm. By the 18th day, the height of the violin has a significant increase, meaning that some samples grow very fast during this period, while some grow slowly. At the moment, the range of the violin is about 1-10mm, i.e., taproot length of most samples is relatively close, about 1-10mm. By the 62nd day, the distribution range of the taproot length has expanded to 0-150mm, and the shape of the violin is very close to a line, indicating that the differences among samples further enlarge, some of the taproots are 150mm long, while others are only about 10mm long. The growth of underground roots also shows similar characteristics, the difference among samples increases with time (Figure 2A). This phenomenon is shown on the violin diagram that the violins are both wide and short in the early stage but narrow and high in the late stage. One potential cause of phenotypic variation may be that those samples have different genotypes, suggesting that there may be QTLs regulating woody growth.
The CV curve of stem length varies widely, with the maximum value reaching 2 mm at the 3rd day (Figure 2D). And it decreases rapidly from the 3rd day to the 18th day, then rises slowly after the 18th day. Above variation trend is caused by the characteristic of sample growth. At the 3rd day, the stem length of most samples is 0, which results in data dispersion. As time goes by, the stems of some samples begin to grow, and the data dispersion begins to decrease. However, there are still some samples with stem length of 0 before the 18th day, so the CV values before the 18th day are always large. After 18 days, basically all samples grow stems, and the CV value is the lowest in the whole growth process. Then, affected by individual differences, some of the samples grow well, while the others grow slowly, with slight phenotypic changes. Therefore, the degree of data dispersion gradually increases, and the CV shows an upward trend. Compared with the stems with slow growth in the early stage, the development of root is earlier, that is, during the growth of seedlings, young individuals tend to allocate more biomass to root growth rather than stem growth. During the whole culture process, the CV value of taproot growth remains at about 1, and decreases slightly after 38 days. The CV curve of root number rises slowly before the 33rd day and remains basically unchanged after the 33rd day. The CV curve of root length decreases slowly before the 33rd day, and remains basically unchanged after the 33rd day, which is completely opposite to that of root number.
3.2 Curve fitting and QTL detection of DIRE
Based on the least square method, we used DIRE to fit the growth trajectory of target traits. Unlike some classical growth equations, such as Logistic equation, Gompertz equation and Richards equation, which describe the target character separately, DIRE considers the growth of two traits as a whole, gives the potential internal interaction pattern between them, that is, how trait 1 promotes or inhibits the growth of trait 2, and conversely, how trait 2 affects the growth of trait 1. DIRE can fit the growth data of stem length, taproot length, root number and root length for most progenies well (Supplementary Figure 1, Supplementary Figure 2, Supplementary Figure 3 and Supplementary Figure 4). The residual distribution diagram shows that the residual value randomly distributed, indicating the robustness of the fitting results (Supplementary Figure 5). Table 1 shows the estimated parameters and evaluation information of the average growth value of the two groups of data by Logistic equation, Gompertz equation and Richards equation and DIRE equation. The evaluation criteria include the residual sum of squares (RSS), the coefficient of determination (), Akaike information criterion (AIC) and Bayesian Information Criterion (BIC). The evaluation result shows that our model has a netter fitting effect than most traditional equations, except that Gompertz equation fits the growth of root number – root length slightly better than DIRE, which indicates that internal interaction indeed affects the growth of the two traits (Table 1).
Our model divides the overall growth of the target trait into two parts, namely, the independent growth part and the interaction part. According to the symbols of , which are determined by the symbols of interaction parameters and , we can judge the interaction strategies between two traits (Table 1; Figure 2C). The interaction strategy of stem and taproot belongs to predation-prey strategy. Stem length growth is inhibited by taproot length growth; conversely, taproot growth benefits from stem length (Figure 2B). The growth pattern of root number and root length follows a win-win cooperation strategy, they promote the growth of each other. However, the degree of mutual benefit is not completely equal. It is obvious that the interactive curve of root length is significantly higher than 0 horizontal line, while that of root number is close to 0 horizontal line, meaning that root length growth gains more benefits from root number growth (Figure 2B).
DIRE quantitatively describes the growth of two characters as a whole interacting with each other. Different growth patterns show different growth curve trajectories, corresponding to different parameter sets . Significant QTLs that regulate the growth of target traits can be identified from the whole genome through a series of hypothesis tests on parameter sets. We obtained 94 significant QTLs regulating the growth of stem length - taproot length and 93 significant QTLs regulating the growth of root number - root length at the threshold of (Figure 3A). These QTLs are sporadically distributed in each linkage group. QTLs regulating the growth of stem length - taproot length mainly locate in Linkage Group 8, 11, 12 and 17, accounting for 10.64%, 22.34%, 10.64% and 12.77% of all key QTLs, respectively (Figure 3B).
Figure 3 Manhattan plot (A) of p-values over 19 linkage groups of Populus euphratica, where the threshold (red horizontal line) is determined by FDR correction. The proportion of significant QTL regulating stem length - taproot length (B) and root number - root length (C) in 19 linkage groups.
A number of QTLs are annotated to biological functions closely related to the growth of Populus euphratica, or are homologous with some genes that play key roles in the growth of other trees. For example, the gene annotated by Q3009 (nn_np_12066, Linkage Group 5) encodes pentapeptide repeat containing protein At4g14850, which is a homologous protein of LOI1 (Zhang et al., 2017). LOI1 is a component of mitochondria and a regulator of isoprenoid biosynthesis, affecting the synthesis of ATPQ, RNA and other biological macromolecules (Kobayashi et al., 2007). LOI1 interacts with MEF14 and participates in mitochondrial editing and cytochrome c, while the cytochrome c, as an important component of the electron transport chain, has the oxidation-reduction ability, participates in the oxidative phosphorylation KEGG pathway, and is active in mitochondria, organelle envelope and cytoplasm. Q2345(nn_np_11278, Linkage Group 4) is homologous with LOC7486368, the gene encoding sugar transport protein (STP) 10 in Poplus tricocarpa. STP is a plant-specific transport protein, which is responsible for absorbing glucose from the apoplast into plant cells. It is an important member of monosaccharide transporters and plays a crucial role in the division and morphogenesis of organs such as seeds. They are the components of organ development in co-plastic isolated tissues, such as seeds, pollen and fruit. (Cheng et al., 2015; Rottmann et al., 2018; Bavnhøj et al., 2021). The key QTLs that control the growth of root number - root length are mainly distributed in Linkage Groups 1 and 4, in which Linkage Groups 1 contains the largest number of QTLs, reaching 39, accounting for 41.94% of all the key QTLs (Figure 3C). These QTLs are also annotated to rich biological functions. In particular, some QTLs are closely related to the growth of plant roots. The gene annotated by Q910 (nn_np_11836) in Linkage Group 1 encodes a leucine-rich repetitive receptor-like protein kinase family protein, which is mainly located in the membrane system. This gene is involved in ATP synthesis, cytokinin regulation and endocytosis transport of cell membrane. In Arabidopsis research, it has been proved that this gene controls root growth by mediating cytokinin. The increase of At2g33170 will promote the growth of root length in Arabidopsis (Hove et al., 2011). Supplementary Table S1 and Supplementary Table S2 show the basic information of these key QTLs, including the linkage group, genetic distance, annotated gene function information, etc.
3.3 Genetic architecture analysis
Genetic factors contribute to phenotypic variation. We calculated heritability i.e. phenotypic variation explained (PVE) of 93 and 94 significant QTLs by quantifying the dynamic genetic contribution of markers to growth (Figure 4). There is dramatic variation in the temporal pattern of heritability (PVE) for stem length, and these patterns can be roughly clustered into three categories (Figure 4A): ① PVE increases and then decreases (pink); ② PVE increases monotonously over time (yellow); ③ PVE decreases consistently over time (green). And the temporal patterns of heritability for taproot length can be sorted into similar three groups (Figure 4B). Among the significant loci, the variation trends of PVE at the same QTL for two traits are different and even completely opposite. Taking Q5033 as an example, for taproot length, PVE of it increases monotonously over time, while PVE of it for the growth of stem length decreases with time, displaying a completely opposite trend. It means that the genetic contribution to the growth of taproot length increases with time, and the genetic contribution to the growth of stem length decreases gradually. At the same time, the PVE variation trend of a pair of characters on the same QTL may be similar, for example, PVE of Q4017 for stem length and taproot length (Figures 4A, B) both increases with time. The temporal patterns of heritability for other traits (Figures 4C, D) can also be sorted into three clusters.
Figure 4 Heat maps and hierarchical clustering trees of heritability are explained by 93 significant QTLs for stem length (A), taproot length (B), and 94 significant QTLs for root number (C), root length (D) of Populus euphratica, respectively. Heat maps of heritability are explained by several selected QTLs. Significant QTLs are clustered into 3 categories (pink, green and yellow).
The growth of phenotypic traits is largely influenced by key QTLs. Although these identified QTLs jointly regulate the overall growth of target trait, there are great differences in their roles and importance in the genetic network. These QTLs do not exist independently, they are subject to epistatic effects from other QTLs besides direct genetic effects produced by themselves. We renumbered these QTLs and constructed their genetic effect network to clarify their interactions and how these QTLs directly or indirectly affect the growth of target traits. These genetic networks, including stem length network, taproot length network, root number network and root length network, share similarities in structure (Figure 5). Several dominant core QTLs regulate most other QTLs, while most QTLs are in a relatively secondary position, they receive epistatic influence from core QTLs in the form of stimulation or inhibition. For example, in stem length network (Figures 5A, B), only 28 core QTLs play major regulatory roles, accounting for 29.8% of all nodes. In the root number network (Figures 5C, D), only 6 QTLs play regulatory roles, accounting for only 6.45% of all QTLs, in which the roles of S42 (Q1538, hk_hk_1218) and S82 (Q7296, nn_np_9326) are particularly critical. The number of links generated by the two QTLs accounts for 79.6% of the total network links, and there are also links among them, which means that the leading QTLs not only produce direct genetic effects, but also indirectly affect the expression of other QTLs, thus directly and indirectly affect the growth of root. In the genetic network of other traits, the situation is similar. On the other hand, genetic networks contain a variety of mutual regulatory relationships between QTLs, including unidirectional and bidirectional relationships through promotion or inhibition, which are determined by the characteristics of genetic effects produced by the QTL itself and the QTL interacting with it. In the taproot network (Figure 5A), there are mutually promotion QTL pairs S16↔S93. In the genetic network of root number (Figure 5C), there is S42↔S20, in which S42 stimulate S20, while S20 inhibits the expression of S42. By counting the number of link type between QTLs in each genetic network, Figure 5B and Figure 5D shows that the number of activation links is more than the number of inhibition links, and the former is 1-2 times that of the latter.
Figure 5 Genetic effect networks (A), bar chart for activating links and inhibitory links and bar chart for outgoing nodes and incoming nodes in networks (B) for stem length and taproot length. Genetic effect networks (C), bar chart for activating links and inhibitory links and bar chart for outgoing nodes and incoming nodes in networks (D) for root number and root length. Each node represents a QTL. The edges connecting nodes represent the interaction among QTLs, the dotted links represent the inhibitory effect, the solid links represent facilitation. The red nodes represent QTLs that stimulate or inhibit other QTLs.
3.4 Simulation
In order to further verify the accuracy and validity of our model, we simulated and analyzed the growth data of Populus euphratica under different sample size (n) and heritability (). The sample size n is set to 100 and 300 respectively, and the heritability , that is, the proportion of genetic variation in the simulated phenotypic variation, is set to 0.05 and 0.1 respectively. According to the analysis results on real data, the growth of stem length and taproot length is regulated by Q681 of linkage group 1, and the growth of root number and root length is regulated by Q5518 of linkage group 11. In the simulation experiment, we conducted 100 effective simulation calculations based on the parameters of the two genotypes of above two significant QTLs under different sample sizes and heritability. Since Q681 and Q5518 are markers of test cross and have two genotypes, we record these two genotypes as ll and lm respectively. Our model is applied to analyze phenotypic data and genotype data obtained from simulation experiment. The Maximum Likelihood Estimation algorithm is implemented to obtain the likelihood function values of the null hypothesis and alternative hypothesis in hypothesis testing and parameter set combined with fourth-order Runge-Kutta, Expectation Maximization (EM), Least Square Method and local search optimization algorithm BFGS. We can find that there are some differences between the parameters obtained from the simulation experiment and the real parameters (Tables 2, 3), but the differences come smaller with the increase of heritability or sample size. Here, the parameter values obtained from the simulation experiment are the average values of 100 groups of parameters.
Table 2 Comparison between the real parameter set and the parameter set of simulation experiment results for stem length (mm) and taproot length (mm).
Table 3 Comparison between the real parameter set and the parameter set of simulation experiment results for root number (count) and root length (cm).
The parameter set of the simulation result is presented in the form of growth curve, which helps us more intuitively observe the similarity and difference between the simulation experiment result and the real growth data (Figure 6). According to the average parameter values of estimated parameters of 100 groups of simulation experiments under four simulation conditions respectively (n=100, 300; =0.05, 0.1), we depicted their corresponding growth curves, including the overall growth curve, independent growth curve and interactive growth curve. By comparing the estimation curves of different heritability levels and different simulated sample sizes with real growth curves, we find that the simulation effect of the estimation curve with heritability of 0.1 is better than that of 0.05, and the simulation effect of the estimation curve with simulation quantity of 300 is obviously better than that of 100. In Figure 6, the simulation curve with larger heritability or sample size is closer to the real curve, especially the curve with heritability of 0.1 and sample size of 300, which is the closest to the real curve. This is consistent with our expected results, that is, the simulation effect increases with the increase of sample size and heritability.
Figure 6 Growth curves of stem length (A), taproot length (B), root number (C) and root length (D) of two genotypes with a heritability of 0.05 and 0.1. The overall growth (solid line) for each trait is decomposed into independent part (broken line) and interactive part (dotted line). The sample sizes are 100 and 300. Black line represents the real curve, red represents the estimation curve with the sample size of 100 and a heritability of 0.05, blue represents the estimation curve with the sample size of 100 and a heritability of 0.1, the green represents the estimation curve with the sample size of 300 and a heritability of 0.05, the orange represents the estimation curve with the sample size of 300 and a heritability of 0.1.
4 Discussion
In the process of tree growth and development, different organs or different tissues of the same organ do not exist independently. We cannot ignore the interaction between them. DIRE not only quantitatively describes cooperation and confrontation between two traits under a certain environment, but also explains the genetic regulation mechanism behind this collaboration. It is calculated on two groups of data with different characteristics: stem length - taproot length, root number - root length. Simulation experiments verify that all groups of data have achieved good results, indicating the validity and universality of our model. The characters that affect and restrict each other in the growth process of tree include but are not limited to above traits. With appropriate modifications and adjustments, DIRE can be extended to describe the relationship between more characters such as height and diameter of stem, stem and leaves, then tap the genetic driving force behind them.
The competition and cooperation between different characters exist in the whole process of tree growth. In this paper, we mainly focus on the seedling stage, which is only a very short period for trees whose life span in years, and its time units are days. But this does not mean that our model is only suitable for describing the growth of the seedling stage. We can extend it to describe the dynamic changes of characters in the process of tree growth for several years or even more than ten years, or consider different growth stages of trees.
The data materials in this paper are cultivated under the condition of controlling environmental variables, but in reality, the environment in which trees grow is uncontrollable, and the genetic effect in the whole growth and development of trees is also affected by environmental stimulus. Therefore, our model can be applied to the response of the genetic structure of tree growth to environmental changes, including in different light, temperature, soil and so on (Srikanth and Schmid, 2011; Xie et al., 2012; Boyce et al., 2020). In particular, by comparing the genetic regulation mechanism of tree growth under some extreme conditions, such as drought, salt stress, etc., it can help to explore the survival mechanism of plants under extreme conditions from a genetic perspective (Wang et al., 2021).
In addition, the model can be extended from two-dimensional to multi-dimensional model, so as to mine the key quantitative trait loci that regulate the synergetic growth of multiple complex traits. It is more reasonable to consider the integrity of plant growth, but this expansion will largely increase the complexity of the model and computational difficulty to a certain extent.
We have carried out simulation experiments based on real data, and the simulation results show that our model has good statistical properties, which provides an effective analytical framework for describing the dynamic interaction patterns between different traits during plant growth and development.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/Lukaiyan/2datasets, github.
Author contributions
KL performed data analysis and wrote the manuscript. X-YZ and RW conceived of the idea, designed the model and contributed to the revision of the manuscript. XW participated in the design of the data analysis. HG and DY participated in data analysis and gene annotation. MY provided data of Populus euphratica, interpretation of data and guidance for data analysis. QF provides guidance for data analysis and simulation. All authors contributed to the article and approved the submitted version.
Funding
This research was supported by the National Key Research and Development Program of China (No. 2021YFB2301601). The work is partially supported by the National Natural Science Foundation of China (Grant No. 32070601) and the Japan Society for the Promotion of Science (Grant No. 19K03613).
Acknowledgments
We thank the members of the Center for Computational Biology, Beijing Forestry University, for their contributions to this project.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2023.1149879/full#supplementary-material
References
Albert, C. H., Thuiller, W., Yoccoz, N. G., Soudant, A., Boucher, F., Saccone, P., et al. (2010). Intraspecific functional variability: Extent, structure and sources of variation. J. Ecol. 98 (3), 604–613. doi: 10.1111/j.1365-2745.2010.01651.x
Bavnhøj, L., Paulsen, P. A., Flores-Canales, J. C., Schiøtt, B., Panyella, B. P. (2021). Molecular mechanism of high affinity sugar transport in plants unveiled by structures of glucose/h+ symporter STP10. Nat. Plants 7, 1409–1419. doi: 10.1101/2020.11.05.369397
Blum, E. K. (1962). A modification of the runge-kutta fourth-order method. Mathematics Comput. 16 (78), 176–187. doi: 10.1090/S0025-5718-1962-0145661-4
Boyce, W. T., Sokolowski, M. B., Robinson, G. E. (2020). Genes and environments, development and time. Proc. Natl. Acad. Sci. 117 (38), 23235–23241. doi: 10.1073/pnas.2016710117
Cheng, J., Wang, Z., Yao, F., Gao, L., Ma, S., Sui, X., et al. (2015). Down-regulating CSHT1, a cucumber pollen-specific hexose transporter, inhibits pollen germination, tube growth, and seed development. Plant Physiol. 168 (2), 635–647. doi: 10.1104/pp.15.00290
Cianciaruso, M. V., Batalha, M. A., Gaston, K. J., Petchey, O. L. (2009). Including intraspecific variability in functional diversity. Ecology 90 (1), 81–89. doi: 10.1890/07-1864.1
Cunniff, J., Purdy, S. J., Barraclough, T. J., Castle, M., Maddison, A. L., Jones, L. E., et al. (2015). High yielding biomass genotypes of willow (Salix spp.) show differences in below ground biomass allocation. Biomass Bioenergy 80, 114–127. doi: 10.1016/j.biombioe.2015.04.020
Do, C., Batzoglou, S. (2008). What is the expectation maximization algorithm? Nat. Biotechnol. 26, 897–899. doi: 10.1038/nbt1406
Freschet, G. T., Swart, E. M., Cornelissen, J. H. (2015). Integrated plant phenotypic responses to contrasting above- and below-ground resources: Key roles of specific leaf area and root mass fraction. New Phytol. 206 (4), 1247–1260. doi: 10.1111/nph.13352
Freschet, G. T., Violle, C., Bourget, M. Y., Scherer-Lorenzen, M., Fort, F. (2018). Allocation, morphology, physiology, architecture: The multiple facets of plant above- and below-ground responses to resource stress. New Phytol. 219 (4), 1338–1352. doi: 10.1111/nph.15225
Fu, Z., Li, W., Xing, X., Xu, M., Liu, X., Li, H., et al. (2016). Genetic analysis of arsenic accumulation in maize using QTL mapping. Sci. Rep. 6 (1), 21292. doi: 10.1038/srep21292
Fu, L., Sun, L., Hao, H., Jiang, L., Zhu, S., Ye, M., et al. (2017). How trees allocate carbon for optimal growth: Insight from a game-theoretic model. Briefings Bioinf. 19 (4), 593–602. doi: 10.1093/bib/bbx003
Fujikawa, H., Munakata, K., Sakha, M. Z. (2014). Development of a competition model for microbial growth in mixed culture. Biocontrol Sci. 19 (2), 61–71. doi: 10.4265/bio.19.61
Gong, H., Zhang, X.-Y., Zhu, S., Jiang, L., Zhu, X., Fang Q and Wu, R. (2021). Genetic architecture of multiphasic growth covariation as revealed by a nonlinear mixed mapping framework. Front. Plant Sci. 12. doi: 10.3389/fpls.2021.711219
Hermans, C., Hammond, J. P., White, P. J., Verbruggen, N. (2006). How do plants respond to nutrient shortage by biomass allocation? Trends Plant Sci. 11 (12), 610–617. doi: 10.1016/j.tplants.2006.10.007
Hove, C. A. T., Bochdanovits, Z., Jansweijer, V. M., Koning, F. G., Berke, L., Sanchez-Perez, G. F., et al. (2011). Probing the roles of LRR RLK genes in arabidopsis thaliana roots using a custom T-DNA insertion set. Plant Mol. Biol. 76 (1-2), 69–83. doi: 10.1007/s11103-011-9769-x
Jaffrézic, F., Thompson, R., Hill, W. G. (2003). Structured antedependence models for genetic analysis of repeated measures on multiple quantitative traits. Genetical Res. 82 (1), 55–65. doi: 10.1017/s0016672303006281
Jiang, L., He, X., Jin, Y., Ye, M., Sang, M., Chen, N., et al. (2018). A mapping framework of competition–cooperation QTLs that drive community dynamics. Nat. Commun. 9 (1), 3010. doi: 10.1038/s41467-018-05416-w
Kobayashi, K., Suzuki, M., Tang, J., Nagata, N., Ohyama, K., Seki, H., et al. (2007). Lovastatin insensitive 1, a novel pentatricopeptide repeat protein, is a potential regulatory factor of isoprenoid biosynthesis in arabidopsis. Plant Cell Physiol. 48 (2), 322–331. doi: 10.1093/pcp/pcm005
Lotka, A. J. (1920). Analytical note on certain rhythmic relations in organic systems. Proc. Natl. Acad. Sci. 6 (7), 410–415. doi: 10.1073/pnas.6.7.410
Magar, B. T., Acharya, S., Gyawali, B., Timilsena, K., Upadhayaya, J., Shrestha, J. (2021). Genetic variability and trait association in maize (Zea mays l.) varieties for growth and yield traits. Heliyon 7 (9), 2405-8440. doi: 10.1016/j.heliyon.2021.e07939
May, R. M. (1973). Stability and complexity in model ecosystems (Princeton, NJ: Princeton University Press).
Miner, B. G., Sultan, S. E., Morgan, S. G., Padilla, D. K., Relyea, R. A. (2005). Ecological consequences of phenotypic plasticity. Trends Ecol. Evol. 20 (12), 685–692. doi: 10.1016/j.tree.2005.08.002
Modrzyński, J., Chmura, D. J., Tjoelker, M. G. (2015). Seedling growth and biomass allocation in relation to leaf habit and shade tolerance among 10 temperate tree species. Tree Physiol. 35 (8), 879–893. doi: 10.1093/treephys/tpv053
Moon, T. K. (1996). The expectation-maximization algorithm. IEEE Signal Process. Mag. 13 (6), 47–60. doi: 10.1109/79.543975
Niklas, K. J., Spatz, H. (2006). Allometric theory and the mechanical stability of large trees: Proof and conjecture. Am. J. Bot. 93 (6), 824–828. doi: 10.3732/ajb.93.6.824
Poorter, H., Niinemets, Ü, Poorter, L., Wright, I. J., Villar, R. (2012). Causes and consequences of variation in leaf mass per area (LMA): A meta-analysis. New Phytol. 182 (3), 565–588. doi: 10.1111/j.1469-8137.2009.02830.x
Poorter, H., Niklas, K. J., Reich, P. B., Oleksyn, J., Poot, P., Mommer, L. (2011). Biomass allocation to leaves, stems and roots: Meta-analyses of interspecific variation and environmental control. New Phytol. 193 (1), 30–50. doi: 10.1111/j.1469-8137.2011.03952.x
Rottmann, T., Fritz, C., Sauer, N., Stadler, R. (2018). Glucose uptake via STP transporters inhibits in vitro pollen tube growth in a hexokinase1-dependent manner in Arabidopsis thaliana. Plant Cell 30 (9), 2057–2081. doi: 10.1105/tpc.18.00356
Sonah, H., O'Donoughue, L., Cober, E., Rajcan, I., Belzile, F. (2015). Identification of loci governing eight agronomic traits using a GBS-GWAS approach and validation by QTL mapping in soya bean. Plant Biotechnol. J. 13 (2), 211–221. doi: 10.1111/pbi.12249
Srikanth, A., Schmid, M. (2011). Regulation of flowering time: All roads lead to Rome. Cell. Mol. Life Sci. 68 (12), 2013–2037. doi: 10.1007/s00018-011-0673-y
Sun, L., Wu, R. (2015). Mapping complex traits as a dynamic system. Phys. Life Rev. 13, 155–185. doi: 10.1016/j.plrev.2015.02.007
Tang, J., Yan, J., Ma, X., Teng, W., Wu, W., Dai, J., et al. (2010). Dissection of the genetic basis of heterosis in an elite maize hybrid by QTL mapping in an immortalized F2 population. Theor. Appl. Genet. 120 (2), 333–340. doi: 10.1007/s00122-009-1213-0
Tumber-Dávila, S. J., Schenk, H. J., Du, E., Jackson, R. B. (2022). Plant sizes and shapes above and belowground and their interactions with climate. New Phytol. 235 (3), 1032–1056. doi: 10.1111/nph.18031
Van der Klein, S. A. S., Kwakkel, R. P., Ducro, B. J., Zuidhof, M. J. (2020). Multiphasic nonlinear mixed growth models for laying hens. Poult. Sci. 99, 5615–5624. doi: 10.1016/j.psj.2020.08.054
Volterra, V. (1928). Variations and fluctuations of the number of individuals in animal species living together. ICES J. Mar. Sci. 3 (1), 3–51. doi: 10.1093/icesjms/3.1.3
Wang, H., Dong, T., Ye, M., Yang, D., Jiang, L., Wu, R. (2021). Modeling genome-wide by environment interactions through comnigenic interactome networks. Cell Rep. 35 (6), 109114. doi: 10.2139/ssrn.3765612
Wang, T., Huang, L., Zhang, X., Wang, M., Tan, D. (2022). Root morphology and biomass allocation of 50 annual ephemeral species in relation to two soil condition. Plants 11 (19), 2495. doi: 10.3390/plants11192495
Wang, H., Leng, C. (2008). A note on adaptive group lasso. Comput. Stat Data Anal. 52 (12), 5277–5286. doi: 10.1016/j.csda.2008.05.006
Wei, D., Cui, K., Ye, G., Pan, J., Xiang, J., Huang, J., et al. (2012). QTL mapping for nitrogen-use efficiency and nitrogen-deficiency tolerance traits in rice. Plant Soil 359 (1-2), 281–295. doi: 10.1007/s11104-012-1142-6
Weraduwage, S. M., Chen, J., Anozie, F. C., Morales, A., Weise, S. E., Sharkey, T. D. (2015). The relationship between leaf area growth and biomass accumulation in arabidopsis thaliana. Front. Plant Sci. 6. doi: 10.3389/fpls.2015.00167
West, G. B., Brown, J. H., Enquist, B. J. (2001). A general model for ontogenetic growth. Nature 413, 628–631. doi: 10.1038/35098076
Wu, R., Cao, J., Huang, Z., Wang, Z., Gai, J., Vallejos, E. (2011). Systems mapping: How to improve the genetic mapping of complex traits through design principles of biological systems. BMC Syst. Biol. 5 (1), 84. doi: 10.1186/1752-0509-5-84
Wu, R., Jiang, L. (2021). Recovering dynamic networks in big static datasets. Phys. Rep. 912, 1–57. doi: 10.1016/j.physrep.2021.01.003
Wu, R., Lin, M. (2006). Functional mapping - how to map and study the genetic architecture of dynamic complex traits. Nat. Rev. Genet. 7 (3), 229–237. doi: 10.1038/nrg1804
Xia, W., Xiao, Z., Cao, P., Zhang, Y., Du, K., Wang, N. (2018). Construction of a high-density genetic map and its application for leaf shape QTL mapping in poplar. Planta 248 (5), 1173–1185. doi: 10.1007/s00425-018-2958-y
Xie, X., Li, S., Zhang, R., Zhao, J., Chen, Y., Zhao, Q., et al. (2012). The bHLH transcription factor MDBHLH3 promotes anthocyanin accumulation and fruit colouration in response to low temperature in apples. Plant Cell Environ. 35 (11), 1884–1897. doi: 10.1111/j.1365-3040.2012.02523.x
Xu, M., Jiang, L., Zhu, S., Zhou, C., Ye, M., Mao, K., et al. (2016). A computational framework for mapping the timing of vegetative phase change. New Phytol. 211, 750–760. doi: 10.1111/nph.13907
Zhang, M. M. (2019). QTL mapping and epistatic analysis of the response of populus euphratica root growth dynamics to salt stress (Beijing: Beijing Forestry University). dissertations.
Zhang, M. M., Bo, W., Xu, F., Li, H., Ye, M., Jiang, L., et al. (2017). The genetic architecture of shoot-root covariation during seedling emergence of a desert tree, Populus euphratica. Plant J. 90 (5), 918–928. doi: 10.1111/tpj.13518
Zhao, W., Chen, Y. Q., Casella, G., Cheverud, J. M., Wu, R. (2005a). A non-stationary model for functional mapping of complex traits. Bioinformatics 21 (10), 2469–2477. doi: 10.1093/bioinformatics/bti382
Keywords: Populus euphratica, differential equation, cooperation-competition, QTL mapping, genetic regulatory network
Citation: Lu K, Wang X, Gong H, Yang D, Ye M, Fang Q, Zhang X-Y and Wu R (2023) The genetic architecture of trait covariation in Populus euphratica, a desert tree. Front. Plant Sci. 14:1149879. doi: 10.3389/fpls.2023.1149879
Received: 23 January 2023; Accepted: 20 March 2023;
Published: 05 April 2023.
Edited by:
Lifeng Xu, Zhejiang University of Technology, ChinaReviewed by:
Milton Brian Traw, Nanjing University, ChinaZitong Li, Commonwealth Scientific and Industrial Research Organisation (CSIRO), Australia
Copyright © 2023 Lu, Wang, Gong, Yang, Ye, Fang, Zhang and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xiao-Yu Zhang, eHl6aGFuZ0BiamZ1LmVkdS5jbg==; Rongling Wu, cnd1QGJqZnUuZWR1LmNu