Dynamic analysis of QTLs on plant height with single segment substitution lines in rice

Dynamic regulation of QTLs remains mysterious. Single segment substitution lines (SSSLs) and conditional QTL mapping and functional QTL mappings are ideal materials and methods to explore dynamics of QTLs for complex traits. This paper analyzed the dynamics of QTLs on plant height with SSSLs in rice. Five SSSLs were verified with plant height QTLs first. All five QTLs had significant positive effects at one or more developmental stages except QTL1. They interacted each other, with negative effects before 49 d after transplanting and positive effects since then. The five QTLs selectively expressed in specific periods, mainly in the periods from 35 to 42 d and from 49 to 56 d after transplanting. Expressions of epistasis were dispersedly in various periods, negative effects appearing mainly before 35 d. The five QTLs brought the inflexion point ahead of schedule, accelerated growth and degradation, and changed the peak plant height, while their interactions had the opposite effects. The information will be helpful to understand the genetic mechanism for developmental traits.

www.nature.com/scientificreports/ owing to the incorporation of biological principles 22 . Studies using these strategies have provided a wealth of information about QTLs controlling the developmental behavior of traits, such as QTL accumulated effects at a point in time, QTL net effects in a period of time, and functional QTLs changing the process of development.
Over the past 30 years, great efforts have been made to dissect the genetic basis of plant height [8][9][10]16,[23][24][25] . Studies have indicated that plant height in rice is generally controlled by both qualitative and quantitative genes, more than 1000 QTLs have been mapped on rice chromosomes, and their action mechanisms manifest mainly accumulation and interaction of QTL effects, i.e., additive effect and dominant or epistatic effect with the features of spatio-temporal expression (http:// www. grame ne. org/ qtl). However, most of the research materials applied in previous studies were conventional mapping populations such as F 2 self-pollinating populations (F 2 ), back-crossing populations (BC), double haploid lines (DHLs), and recombinant inbred lines (RILs), in which inconsistent genetic backgrounds among individuals or lines could disturb the analysis results 26 . Moreover, systematic studies on developmental traits have rarely been reported using various approaches simultaneously, especially those lacking to tail after the expression pattern for an identified QTL. In this paper, we applied five single segment substitution lines (SSSLs) as experimental materials to first identify whether putative QTLs were carried out on plant height for each SSSL, to then estimate conditional effects y(t|t − 1) and four functional parameters in Wang-Lan-Ding's model and to finally carry out unconditional, conditional and functional QTL mapping. Here, the conditional effects were estimated based on phenotypic values of plant height measured at nine time points of development. Wang-Lan-Ding's model is about the second-order ordinary differential equation of insect developmental rate with respect to temperature 27 (Wang et al. 1982). The model, similar to the well-known logistic model, was verified to be able to more accurately describe the entire developmental process of insects 28 and can also be applied to fit the growth curve of developmental traits in crops 29 . We aimed to provide the first quantification of genetic patterns affecting plant height, including QTL effects and interactions, accumulated effects and net effects, as well as how functionally to regulate the growth curve of plant height, and a comprehensive understanding of dynamic genetic mechanisms for developmental traits such as plant height.

Materials and methods
Materials and mapping population. Similar experimental materials as a previous study 30 , Huajingxian 74 (HJX74) and its five single segment substitution lines (SSSLs) ( Table 1), were applied in this trial. HJX74 is an elite indica variety with many excellent properties cultured by our laboratory in South China. Each SSSL possessed only a single substituted segment from a donor under the HJX74 genetic background and was distributed in related molecular marker regions on corresponding chromosomes with given lengths (Fig. 1). Double segment substitution lines (DSSLs) were polymerized based on the F 2 populations from the crossing combinations of two SSSLs. A half diallel mating scheme was conducted using HJX74 and their SSSLs, DSSLs as crossing parents to generate the mapping population that included HJX74, 5 SSSLs, 7 DSSLs and 26 crossing combinations. Some crossing combinations were lacking since the seeds of F 1 were scarce.
Field experiments and measurement of plant height. The same phenotypic trial as the previous study 30 was applied in this study. The trial site was located at the teaching and experimental station of South China Agricultural University in Guangzhou, China (23° 79ʹ N, 113° 159ʹ E). In the early season (duration from March to July) of 2018, 39 genotypic materials were grown in a completely randomized block design with three replications. Germinated seeds were sown in a seedling bed, and then seedlings were transplanted to a rice field 20 days later, with one plant per hill and a density of 16.7 cm × 16.7 cm. Each plot consisted of four rows with ten plants per row. Local standard practices were used for the management of the trial. Plant height per hill on 10 central plants, the distance between the base and the top of the main stem of the plant, was measured in each plot from seven days after transplanting onwards, and data every 7 days once were continuously recorded for nine weeks (denoted by t 1 to t 9 ). The averages of plant height in each plot for the nine stages were used as input data for the subsequent analysis.
Statistical analysis and estimation of QTL effects. The model y ij = µ + G i + B j + e ij was adopted to analyze the variances on the phenotypic performances (y) of plant height measured for each plot at various stages, where µ, G, B and e were the estimations of the population mean, genotype, block and residual error effect, respectively. i and j represent serial numbers of genotypes and blocks, respectively. The additive effect (a) or dominant effect (d) of QTLs was calculated by the estimations of (S i − HJX74) , and the epistatic effect (e) between QTL pairs was calculated by (D ij + HJX74 − S i − S j ), where S, D and HJX74 indicate the phenotypic means of single-, dual-segment substitution materials and HJX74, respectively, and i, j represents two homozygotes or heterozygotes of SSSLs.
Conditional variable y t|t−1 was estimated by the formula of y t|t−1 = y t − b t/t−1 (y t−1 − y t−1 ) , presenting phenotypic value at time t conditional on the phenotypic value at given time t − 1, where y t−1 , y t−1 and y t , ⇀ y t were the phenotypic values and the means at times t − 1 and t, respectively. b t/t−1 is the regression coefficient for phenotypic values at time t versus those at time t − 1. QTL mapping was imposed on the conditional variable y t|t−1 to generate conditional QTLs. Statistical analysis and estimation of QTL effects were carried out with aov ( ) and lm ( ) functions in R language (https:// www.r-proje ct. org/).
The study was in compliance with relevant institutional, national, and international guidelines and legislation.

Results
Unconditional QTL mapping on plant height. Plant height approximately approached the "S" type of growth curve (Fig. 2). The figure drawn from all 39 genotypic materials indicated that after slow growth, plant height reached its peak and then started to decrease slightly. Separate analysis of variances at various stages revealed that a significant difference in plant height existed among genotypes (Supplementary Table 1), supporting the existence of QTLs for plant height in the mapping population. The contrast tests found that each of the five SSSLs harbored plant height QTLs ( Table 2). All carried with additive and/or dominant effects detected at one or more stages (see QTL(t) in Table 2). QTL on SSSL S 5 (denoted by QTL 5, similarly hereinafter) detected with significant dominant effects just at one of stages perhaps was unreliable. The other QTLs www.nature.com/scientificreports/ repeatedly appeared guaranteed their truth. Only QTL 1 exhibited negative effects, and the others increased plant height. During the early period (from t 0 to t 3 ), few QTLs were detected, while more QTLs were present in the middle-late period. The variations in QTL effects with time implied the dynamics of expression for these QTLs.
To understand the interaction mechanism among these QTLs, we first partially aggregated two SSSLs to generate a dual segment substitution line (DSSL) and then carried out a half diallel mating scheme to achieve the various genotypes required to estimate epistasis. Four epistatic components, additive-additive (aa), additivedominance (ad), dominance-additive (da), and dominance-dominance (dd), were estimated according to configurations of genotypes for seven QTL pairs (Table 3). All seven pairs of S i /S j holding significant interaction effects further confirmed the prevalence of epistasis (see QTL(t) in Table 3). Two epistatic components, d 1 d 2 (denoted dd of S 1 /S 2, similarly hereinafter) and a 1 a 5 were detected only at one of stages, and reliability was subject to further verification. The other epistatic components were significant at least at two stages, which indicated the validity of these interactions. Negative epistatic effects were major, while positive epistases mainly appeared at the periods of t 8 and t 9 . The causes of negative (positive) epistases perhaps were due to positive (negative)  www.nature.com/scientificreports/ additive and dominance effects of QTLs, which will be discussed later. All epistatic effects dynamically changed with developmental stages (see QTL(t) in Table 3).

Conditional QTL mapping on plant height.
To acquire the net effect of a given QTL on plant height during a certain period, we carried out conditional QTL mapping. The conditional effects y (t|t−1) , the net effects of phenotypic values at time t given the phenotypic values at time t − 1 , were first estimated (Supplementary  Table 2). Then, conditional QTL effects, the net effects of QTLs from time t − 1 to time t , were calculated based on the conditional effects y (t|t−1) (see QTL(t|t − 1) in Table 2). Conditional QTLs revealed the quantities and stages of QTL expression. QTL 1 had two expressions, one was in the stage from t 5 to t 6 , exhibiting a − 6.57** additive effect, and the other was from t 7 to t 8 with a − 3.75* additive effect and a − 4.96** dominant effect. Similarly, QTL 4 expressed a 3.39* dominant effect from t 3 to t 4 , a 5.60** additive effect and a 3.82* dominant effect from t 5 to t 6 , and a QTL 2 5.03** dominant effect from t 8 to t 9 . Although QTL 2 , QTL 3 and QTL 5 exhibited significant accumulated effects at certain stages, the concentrated expression stages of these QTLs were not detected due to insufficient large expression quantities. There were no significant expressions to be detected in the early period (from t 0 to t 3 ), and QTL expressions occurred mainly in the middle period (from t 3 to t 6 ) and the late period (from t 6 to t 9 ). Table 3. Epistatic effects between QTLs estimated on plant height at various developmental stages. SSSL is the abbreviation of the single segment substitution line. S i represented ith SSSL. aa, ad, da and dd were the abbreviations of epistatic components of additive-additive, additive-dominance, dominance-additive and dominance-dominance, respectively, estimated by the mean of (D ij + HJX74 − S i − S j ) (where D ij , S i , S j indicated dual segment and its two single segment materials, respectively, which might be homozygotes or heterozygotes). QTL(t) and QTL(t|t − 1) represent unconditional and conditional QTLs at stage t, respectively. t i indicates various developmental stages 7 days apart. Sign "-" indicates descending plant height due to alleles from donors. Superscripts "* and **" indicate significance at the 5% and 1% levels, respectively. www.nature.com/scientificreports/ QTL interactions also exhibited different dynamic models (see QTL(t|t − 1) in Table 3). In the early, middle and late periods, there were six, seven and thirteen significant epistatic expressions, respectively. In the three periods, t 2 − t 3 , t 4 − t 5 and t 6 − t 7 , QTLs hardly expressed. There were 14 significant positive epistatic effects and 12 negative epistatic effects. Mostly, negative expressions appeared in the early period, while positive expressions appeared in the late period. Some epistatic components had significant accumulated effects at certain stages, but their expression periods were not detected due to dispersed expression being insignificant. Conversely, some epistatic components had significant net effects in certain stages but did not detect significant accumulated effects due to reverse expressions. Many epistatic expressions were feeble and failed to be detected, while some large expressions became invisible because of large errors.
Functional QTL mapping on plant height. Functional QTL mapping is an appropriate method that passes a mathematical equation to describe a biological developmental process with the genetic mapping framework 18 . We first applied the Wang-Lan-Ding model 27 to fit curves of plant height and to estimate four functional parameters--the optimum time (t 0 ) , the growth rate (r) , the maximum value (K) and the degradation rate (c)(Supplementary Table 3). Then, based on these estimations, we carried out QTL mapping ( Table 4). The five SSSLs were found to harbor QTLs with additive and/or dominance to regulate the four parameters. Any one of the SSSLs was associated with at least two functional parameters, for which pleiotropy or close linkage of genes were responsible. S 1 and S 5 were involved in all four parameters, and S 2 , S 3 and S 4 regulated t 0 or K and c . For the parameter t 0 , QTLs shortened the time of the inflexion point on a curve. For r and c , QTLs improved not only the growth rate but also the degradation rate. The impact of QTLs on parameter K differed, enabling the potential of plant height to increase or decrease.
All pairs of S i /S j had significant interaction effects, involving one or more parameters by various epistatic components. Interactions between QTL 2 and QTL 3 , QTL 4 , and QTL 5 influenced one, two and three parameters, respectively. Epiststatic interactions between QTL 1 and the other QTLs were associated with all four parameters. Epistases always regulated t 0 and K positively, while r and c negatively ( Table 4). The relationship in which positive (negative) epistasis was always derived from the interaction of negative (positive) QTLs was confirmed once again.

Discussions
Dynamic mapping for developmental traits. Most traits of agricultural importance are under the control of an interacting network of genes, which grow and develop through the dynamics of gene expression during the whole growth period 13,14 . Studies via QTL mapping on different kinds of data can provide a wealth of information about the dynamics of QTLs regulating the developmental behavior of quantitative traits 11 . Unconditional QTL mapping of the directly measured phenotypes at various developmental stages can detect the accumulated effects of QTLs before a certain static time point but cannot estimate the expression quantity due to the correlations of phenotypes between two time points 13,14 . Conditional QTL mapping on the indirect estimations of conditional phenotypes can provide net expression of QTLs in a time interval and the stages of QTL expression 16,17 . Functional QTL mapping of the parameters defined in a mathematical function that describes trait variation with biological significance can reveal the QTLs regulating the shape and trajectory of developmental curves [18][19][20][21][22] . In this paper, we carried out a systematic analysis of the dynamics of QTLs regulating the developmental behavior of plant height in rice. We detected that the five SSSLs carried significant additive and/ or dominant effects of QTLs on plant height at multiple developmental stages. These QTLs were credible due to their repeatability. Except for QTL 1 , which exhibited negative effects in the middle-late periods, the other QTLs showed enhanced plant height (see QTL(t) in Table 2). These QTLs interact with each other to form a genetic network to regulate plant height. The seven pairs of SSSL combinations tested all had one or more significant epistatic components to mostly reduce plant height (see QTL(t) in Table 3). QTLs were characterized by temporal expression, selectively appearing significant effects in specific stages of development. The five QTLs turned on mainly in the middle-late periods (see QTL(t|t − 1) in Table 2), whereas the seven QTL interactions dispersed in various periods (see QTL(t|t − 1) in Table 3). Some expressions were too small to be statistically detected. Plant height varied approximately, followed by the logistic curve of the Wang-Lan-Ding model (Fig. 1), which was determined by the parameters of t 0 , r , K and c (Wang et al. 1982). These parameters changed the trajectory of the growth curve of plant height, including the inflection point, the growth rate, the peak value and the degradation rate. Our research indicated that the four functional parameters were regulated by the QTLs and the QTL interactions on the five SSSLs, each of which regulated at least two parameters (Table 4).

Dynamic patterns of QTL expressions.
One of the major goals in developmental genetics is to explore gene expression 13,14 . Conditional QTL mapping makes this possible, which can estimate the net expression of QTLs in a certain time interval 16,17 . In theory, the unconditional QTL effect at time point t is the accumulation effect of QTLs from initial time to time t , which can be divided into several conditional QTL components, i.e., QTL t = QTL 1 + QTL 2|1 + QTL 3|2 + · · · + QTL t|t−1 . Conditional QTL effects were independent of each other and thus were additive. According to the formula, it is possible to generate follow a few cases at stage t , being significant for both QTL t and QTL t|t−1 , either QTL t or QTL t|t−1 , and neither QTL t nor QTL t|t−1 . The relationship between unconditional QTLs and conditional QTLs was discussed in our previous paper 30 and was well validated by the results estimated in this paper. The correlation coefficient between QTL effects at the final stage t 9 and the sum of all conditional QTL effects before t 9 reached 0.9379** in the previous paper 30 and 0.7208 ** in this paper (data not shown). Only a series of conditional QTLs truly reflected the expression periods and quantities of a QTL throughout the whole developmental stage. Conditional QTL mapping has widely been applied to reveal dynamic gene patterns for developmental traits 8 www.nature.com/scientificreports/ for the genetic control of growth trajectories, permanent QTLs, early QTLs, late QTLs and inverse QTLs 21 . This knowledge derived from the accumulated effects of QTLs, QTLs being permanent, early, late and inverse when one genotype was better than the other in entire growth process, at early stages, at late stages and one genotype showed inverse effects with the other since a particular stage, respectively. However, the accumulated effects of QTLs could not reflect the expression stages and quantities of QTLs. This paper indicated that QTLs for plant height expressed all additive, dominant and epistatic effects according to the temporal expression pattern (see QTL(t|t − 1) in Tables 2 and 3). QTLs and their interactions expressed significant effects only in one or more periods and sometimes even had no significant expression periods while remaining silent all the time. Permanent expression of QTLs was rare. QTL 1 and QTL 2 were expressed mainly in the late period, QTL 4 was expressed in the middle period, and QTL 3 and QTL 5 were not significantly expressed (see QTL(t|t − 1) in Table 2). Similarly, QTL 1 /QTL 2 , QTL 1 /QTL 3 , QTL 1 /QTL 5 and QTL 2 /QTL 3 were expressed mainly since period t 5 , QTL 1 /QTL 4 Table 4. QTL effects for the four functional parameters on plant height. SSSL is the abbreviation of the single segment substitution line. S i represents the homozygote or heterozygote ith SSSL. The additive effect (a) or dominant effect (d) of QTLs was estimated by the mean of (S i − HJX74) , where HJX74 was the abbreviation of Huajingxian 74. aa, ad, da and dd represented the additive-additive, additive-dominance, dominanceadditive and dominance-dominance epistatic components, respectively, which were estimated by the mean of (D ij + HJX74 − S i − S j ) (where D ij , S i , S j indicated dual segment and its two single segment materials, respectively, which might be homozygotes or heterozygotes). t 0 , r, K and c were the optimum time, the growth rate, the maximum value and the degradation rate, respectively. Sign "-" indicates descending parameters due to alleles from donors. Superscripts "* and **" indicate significance at the 5% and 1% levels, respectively. www.nature.com/scientificreports/ and QTL 2 /QTL 4 were dispersed in various periods, and QTL 2 /QTL 5 had inverse effects between the early period and the late period (see QTL(t|t − 1) in Table 3). In fact, QTLs and their interactions expressed net effects in various stages, just some of which reached the levels of significance statistically. Small expressions of QTLs were considered as no expressing or experimental error.

QTLs regulated developmental trajectories of temporal traits. Developmental theory considers
that if different genotypes at a given QTL correspond to different developmental trajectories, the QTL must affect the differentiation of this trait 21 . Therefore, by estimating the functional parameters that define the trait curve of each QTL genotype and testing the differences in these parameters among genotypes, one can determine whether a QTL affects the formation and expression of a trait during development. In the Wang-Lan-Ding model, there were four parameters to regulate the growth curves of developmental traits, which might change the inflexion point ( t 0 ), the upper limit ( K ), the rise speed ( r ) and the descent speed ( c ) of curves 27 . In this paper, genotypes of five SSSLs differed from that of HJX74 at a given QTL (Fig. 3), implying that a putative QTL existed on each of the SSSLs. Both unconditional QTL mapping and conditional QTL mapping confirmed the existence of QTLs (Table 2). How did these QTLs affect the development of plant height? Functional QTL mapping based on the estimations of the four parameters indicated that QTL 1 and QTL 5 regulated all four parameters by additive and/or dominant effects, and the other three QTLs influenced two of them, t 0 or K and c , respectively. These QTLs brought the inflexion point ahead of schedule and accelerated the growth and degradation of plant height. QTL 1 and QTL 5 made the maximum plant height shorter, while QTL 3 and QTL 4 had higher plant heights (Table 4). Similarly, the interactions among these QTLs also influenced the four parameters, which always positively regulated t 0 and K , while r and c negatively regulated by various epistatic components ( Table 4).

Impact of epistasis on plant height.
In multiple gene systems, gene interactions are inevitable except for gene additives, which include allelic interactions (dominance) and nonallelic interactions (epistasis) 32,33 . For statistical purposes, genotypic effects can be divided into single locus effects (additive or dominance) and interaction effects among segregating loci (epistasis) 34 . Thus, epistasis is an important genetic component and a plausible feature of the genetic architecture of quantitative traits. Mapping epistatic interactions is challenging experimentally, statistically and computationally, which requires large sample sizes, severe penalties and a large number of tests 35 . QTL mapping studies using primary mapping populations such as F 2 , BC, DHLs and RILs cannot clearly support the existence of specific interactions among QTLs because of the limitation of inconsistent genetic backgrounds among individuals or lines 26 . SSSLs or NILs (near isogenic lines) have huge advantages for QTL identification in general and characterization of epistasis in particular 26,36,37 . On the one hand, target QTLs can be detected by testing the difference between any one of the SSSLs and the receptor parent; on the other hand, pyramiding of SSSLs enables estimation of epistatic effects [38][39][40][41][42][43][44] . In this paper, we first detected five SSSLs that carry significant additive and/or dominant effects of QTLs on plant height (Table 2), and then four components of epistases were estimated via analysis of pyramiding effects derived from two SSSLs ( Table 3). The information is reliable due to the repeated emergence of putative QTLs and their interactions at multiple stages of development. All seven pairs of SSSLs were tested with two or more significant effects of four epistatic components, further confirming the prevalence of epistasis (Table 3). Epistasis may be brought about by modification of gene function due to alterations in the signal-transducing pathway. Epistatic genes are more deleterious in combination than separately, which are often accompanied by inverse epistatic interactions as homeostatic (that is, canalizing) mechanisms 35 . This role of epistasis was first observed by Eshed and Zamir 37 when they found lessthan-additive interactions between QTLs in tomato. We also confirmed an inverse epistastic role of yield traits in rice, negative epistasis being derived mainly from interactions between positive QTLs, while positive epistasis from negative QTLs 30,43 . In this paper, most QTLs were detected with positive additive and/or dominance; thus, the estimations of epistatic components were mainly negative. Positive epistasis occurred only after the stage of t 6 since negative QTLs appeared (Tables 2 and 3). The property of epistasis was stipulated by the calculated formula e ij = g ij − a i − a j , where e, g, a indicated epistatic, genotypic and two of single locus effects, respectively, and i, j represented two loci. Because the value of e ij is inversely proportional to the sum of single locus effects, it is more likely to gain negative (positive) epistasis when there are positive (negative) single locus effects. Additionally, one QTL always interacted with multiple other QTLs, forming a genetic network. In this paper, five QTLs were detected to interact with each other with one or more significant epistatic components (Table 3). Of the seven combinations of SSSL sets, QTL 1 and QTL 2 interacted with the other four QTLs, while QTL 3 , QTL 4 , and QTL 5 interacted with at least the other two QTLs. Five QTLs had various interaction magnitudes, displaying different epistatic intensities. QTL 2 and QTL 4 seemed to have larger average epistatic effects and greater interoperability than the other QTLs (Supplementary Table 4). Of the four epistatic components, the average estimation of dd was seemingly larger than those of the others (Supplementary Table 4). Knowledge of epistatic interactions will improve our understanding of genetic networks and mechanisms that underlie genetic homeostasis and improve predictions of responses to artificial pyramiding breeding for quantitative traits in agricultural crop species. In the future, we must assess the effects of pairwise and higher-order epistatic interactions between polymorphic DNA variants on molecular interaction networks and, in turn, evaluate their effects on organismal phenotypes to understand the mechanistic basis of epistasis 35 . Only then will we be able to go beyond describing the phenomenon of epistasis to predicting and testing its consequences for genetic systems.