Genetic analysis of Upland cotton dynamic heterosis for boll number per plant at multiple developmental stages

Yield is an important breeding target. As important yield components, boll number per plant (BNP) shows dynamic character and strong heterosis in Upland cotton. However, the genetic basis underlying the dynamic heterosis is poorly understood. In this study, we conducted dynamic quantitative trait loci (QTL) analysis for BNP and heterosis at multiple developmental stages and environments using two recombinant inbred lines (RILs) and two corresponding backcross populations. By the single-locus analysis, 23 QTLs were identified at final maturity, while 99 QTLs were identified across other three developmental stages. A total of 48 conditional QTLs for BNP were identified for the adjacent stages. QTLs detected at later stage mainly existed in the partial dominance to dominance range and QTLs identified at early stage mostly showed effects with the dominance to overdominance range during plant development. By two-locus analysis, we observe that epistasis played an important role not only in the variation of the performance of the RIL population but also in the expression of heterosis in backcross population. Taken together, the present study reveals that the genetic basis of heterosis is dynamic and complicated, and it is involved in dynamic dominance effect, epistasis and QTL by environmental interactions.


Relationship between heterozygosity and dynamic performance. Most of correlation coefficients
are not significant between heterozygosity of whole-genome and the dynamic performance of BC(V) and MPH data in terms of BNP trait at four developmental stages in two hybrids ( Table 3). The result suggested that overall genome heterozygosity alone possessed little effect on dynamic trait performance. Heterosis might derived from just a small of genome heterozygosity in the BC(V) data at not only early stage but also maturity of cotton plant. The low correlation coefficients may be the result of low density loci for the whole genome and only half of heterozygosity of whole-genome existed in backcross population. These results support the previous studies 4, 33,34 . Single-locus QTL. Genetic maps for the two hybrid populations were previously constructed by the polymorphic loci ( Figure S1). For the XZ hybrid, 623 loci were mapped to 32 linkage groups and the genetic map spanned 3,889.9 cM. For the XZV hybrid, 308 loci were mapped to 39 linkage groups and the genetic map spanned 3,048.4 cM 17 .
Scientific RepoRts | 6:35515 | DOI: 10.1038/srep35515 QTLs detected using composite interval mapping in XZ and XZV hybrids are shown in Table S2 at the single-locus level. Table 4 listed numbers for different effect of QTLs identified by composite interval mapping at four developmental stages at three environments in backcross population. A total of 38 and 22 QTLs were respectively detected for BNP at four developmental stages of XZ and XZV hybrids.
In the XZ hybrid, a total of 38 QTLs were detected in RIL, BC, and MPH data sets (Table S2). Twenty-one QTLs were identified in more than two developmental stages or environments or populations. Totally, 22 and 14 QTLs were respectively detected at stage t1, t2, t3, and t4 in three data sets. In backcross population, 19 QTLs with an additive effect, nine with partial dominant effect, and five with overdominant effect were observed. In the XZV hybrid, a total of 22 QTLs were detected in RILV, BCV, and MPH data sets (Table S2). Eleven QTLs were identified in more than two developmental stages or environments or populations. Interestingly, we observed one QTL qBNP-Chr23-3 for BNP which were detected at all four time points. Totally, 12, 14, 13, and 9 QTLs, were respectively detected at stage t1, t2, t3, and t4 in three data sets. In backcross population, eleven QTLs with an additive effect, eight with partial dominant effect, and twelve with over dominant effect were observed. Conditional QTL refers to the QTL expression at the specific stage from t-1 to t at the locus 31,32 . QTL identified at t stage conditioned on t-1 stage revealed net effects of QTL expression at the stage from t-1 to t stage. A total of 34 and 14 conditional QTLs for BNP were identified for the adjacent stages in XZ and XZV hybrids respectively (Tables S3 and S4). In the XZ hybrid, nineteen, eight, and seven conditional QTLs were detected at Δ t1-2, Δ t2-3, and Δ t3-4 stages respectively; seven, five, and two conditional QTLs were respectively detected at Δ t1-2, Δ t2-3, and Δ t3-4 stages in the XZV hybrid. For MPH data set, eight, four, and three conditional QTLs were respectively found at Δ t1-2, Δ t2-3, and Δ t3-4 stages in two hybrids. These results indicated that QTL controlling BNP and heterosis might expressed more actively at the early stage. No conditional QTL and heterotic QTL (hQTL) was simultaneously identified during the entire stage of growth. Most of the conditional QTL detected at special stage showed that the QTL and hQTL expressed selectively at certain stage during reproductive growth.

QTLs and QE interactions resolved by two-locus analyses. A total of 110 and 70 M-QTLs and QEs
were respectively detected by inclusive composite interval mapping (ICIM) at four developmental stages of XZ and XZV hybrids ( Table 5, S5, S6 and S7). In the XZ hybrid, a total of 57, 36 and 17 M-QTLs and QEs were detected in the RILs, BC hybrids, and MPH data at all of developmental stages, respectively. On average, the M-QTL explained 2.52%, 1.84%, and 1.08% of the phenotype variation (PV), and the QEs explained 0.51%, 0.73%, and 1.32% of the phenotype variation in the RILs, BC hybrids, and MPH data respectively. In the XZV hybrid, a total of 26, 21 and 23 M-QTLs and QEs were respectively detected in the RILVs, BCV hybrids, and MPH data at all of developmental stages. On average, the M-QTL explained 2.86%, 1.43%, and 1.87% of the phenotype variation, and the QEs explained 0.75%, 1.29%, and 0.80% of the phenotype variation in the RILVs, BCV hybrids, and MPH data, respectively.
Totally, at two-locus levels, 215 and 393 E-QTLs and QEs were detected by ICIM in three data sets of XZ and XZV hybrids at four developmental stages respectively (Table 5, S8, S9 and S10). In the XZ hybrid, a total of 145, 54, and 16 E-QTLs and QEs were detected in the RILs, BC hybrids, and MPH data, respectively. On average, the E-QTL explained 3.02%, 2.88%, and 1.78% of the phenotype variation, and the QEs explained 0.60%, 0.91%, and 1.56% of the phenotype variation in the RILs, BC hybrids, and MPH data respectively. In the XZV hybrid, a total of 288, 40 and 65 E-QTLs and QEs were detected in the RILVs, BCV hybrids, and MPH data, respectively. On average, the E-QTL explained 3.44%, 1.61%, and 2.40% of the phenotype variation, and the QEs explained 0.64%, 1.70%, and 1.19% of the phenotype variation in the RILVs, BCV hybrids, and MPH data respectively.

Discussion
Yield is an important breeding target in Upland cotton. Biologically, plant height equals to all of internodes' lengths above the ground, showing the rate of vegetative growth in Upland cotton 31    performance were observed for four developmental stages in three environments (Table S11). Early fruit number has indication of BNP at maturation. Meanwhile, BNP trait as yield components shows the best contribution to the yield (Table S11). Dynamic traits varied considerably among individuals in two hybrids. The incremental value of RIL and BC populations at early developmental stages were more than that at the final stage. The incremental values of the parents, RIL and BC populations displayed a dynamical developmental process during plant growth. The MPH (%) of backcross population gradually reduces with the plant growth, and heterosis performance showed time dependent and dynamic character. The BNP dynamic performance was conducive for studying the dynamics of heterosis in Upland cotton. The genetic control underlying the dynamic heterosis is poorly understood in crop. The hQTL mapping provided information of cumulative effects at special stage in previous studies 5,35 . Based on final value of a quantitative trait, the genetic effects of QTLs at different developmental stages were overlooked 31,32,36 . In present research, a total of 14 and 9 QTLs for BNP were identified at maturity in XZ and XZV hybrids, while 60 and 39 QTLs were identified across other three development stages other than maturity respectively. It is obvious that lot of hQTLs identified at early stage were not identified at maturity stage (Table 5). Similar results were observed for the developmental behavior of tiller number in rice 37 and plant height in Upland cotton 31 . Moreover, our results showed that the QTL and hQTL possessed temporal character 30,31,38,39 . At the single-locus level, the most number of QTLs was detected at t2 stage, and the t2 stage was the period of rapid reproductive stages. This phenomenon was supported by the result that the most incremental values and heterosis of BC population were observed at t1-t2 stages. In addition, the most number of conditional QTL and hQTL detected at t1-2 stage also verified this phenomenon. It is possible that more genes controlling heterosis of BNP trait largely expressed at the period of rapid reproductive stage in two hybrids. These results mean that the t1-t2 stage might be a more reasonable time for fine mapping of the hQTL.
Several QTLs at t1 and t4 stages in the MPH data set overlapped with QTLs in the RIL(V) and BC(V) populations (Table S12). At t2 and t3 stages, few QTLs were identified in the MPH data set that were also observed in the RIL(V) and BC(V) populations (Table S12). Previous studies have suggested that hQTLs of yield-related traits were independent and were different between QTLs and hQTL 22,40,41 . While, a recent study showed that many QTLs detected for grain yield in the MPH data set overlapped with QTL in the IF 2 population and hQTLs were not independent 42 . In present study, we observe that hQTLs were not independent at the early and maturity stages, however, hQTLs were independent at rapid growth stage of t2 and t3. It suggested that phenotypes and heterotic traits might be controlled by different loci at different developmental stages in Upland cotton. In addition, we observed that the QTL qBNP-Chr23-3 which showed high LOD value and phenotype variation at early stage can be repeatedly identified at four stages in multiple populations and environments in the XZV hybrid. Like the results from dynamic QTL for plant height at different developmental stages of rapeseed, we might predict this QTL at early developmental stages with no need for measurement of BNP trait at maturity stage 43 . These results added additional information and supplied important reference on QTL mapping for BNP and heterosis performance at different developmental stages.
Permanent genetic population such as RIL and permanent BC populations is required in order to detect convincing QTLs in more than one population or environment 4 . Using RIL and BC populations could mutually increase the power of detecting QTLs 34 . The RIL population appears to be able to detect more QTLs than BC population in the XZ hybrid (Table S2). Similar results were available in previous studies in rice 34 and oilseed rape 10 . In the XZV hybrid, the opposite happened. In present study, the QTLs detected with the MPH data were relatively fewer than that of the QTLs identified with the RIL(V)s and BC(V) data. The reason may be that QTLs with an intermediate mode of inheritance and the QTL with lacking dominance effect could not be identified in the MPH data 10 . In previous study, three QTLs for BNP were identified in F 2:3 and F 2:4 populations derived from Upland cotton cross 'GX1135' × 'GX100-2' . Present RIL population of F 9 generation consisting 177 RILs along with corresponding backcross population were developed from F 2: 3 population by single seed descent 17 . These three QTLs for BNP identified previously were once again detected simultaneously using RIL and backcross   population. A large number of same QTLs between XZ and XZV hybrids were not detected. One possible reason is that the density of genetic map of the XZV hybrid is low, compared with the XZ hybrid. A total of 21 and 11 stable QTLs were identified in more than two developmental stages or environments or populations in XZ and XZV hybrids respectively. Specially, qBNP-Chr23-3 could be detected across four developmental stages in the XZV hybrid, and the QTLs in different developmental stages showed same direction of genetic effect. These stable QTLs might be important for marker-assisted selection (MAS) for developing varieties with high BNP and yield in Upland cotton. Pyramiding the stable targeted loci is an efficient and feasible strategy for improving BNP and yield in Upland cotton. The QTL identified simultaneously with RIL(V)s, BC(V) and MPH data allowed an assessment of the degree of dominance in two hybrids at the single-locus level 10 . In the XZV hybrid, a total of one (25.0%) QTLs for a PD effect and three (75.0%) QTLs for an OD effect were identified at t1 stage. At final t4 stage, a total of two (100.0%) QTLs for a PD effect and zero (0.00%) QTLs for an OD effect were identified. It is obvious that overdominance played a relatively more important role in controlling heterosis than partial dominance at early stage, while the result was opposite at final maturation stage. The study of maize heterosis showed that QTL for trait with low heterosis mainly existed in the additive to dominance range and QTL for trait with high heterosis had effects with the dominance to overdominance range 9 . Grain yield with the highest level of heterosis has the largest number of QTL exhibiting overdominance in all traits measured in rapeseed heterotic study 10 . In our study, BNP trait show the high level of heterosis in t1 and t2 stage, and the largest number of QTL exhibiting overdominance in all traits measured in rapeseed heterotic study. These results showed that dynamic dominance and overdominance played a role in controlling the expression of heterosis during growth in Upland cotton.
In XZ and XZV hybrids, a host of epistatic interactions and QEs were observed with the three data sets at four developmental stages. On average, the variation explained by E-QTL for most of stages was much greater than that by M-QTL in three data sets (Table 5). These results suggest that universal epistasis played an important role not only in the variation of the performance of the RIL(V) population but also in the expression of heterosis in BC(V) population 10,34 . Epistasis as the genetic basis of heterosis was supported by a series of previous studies 17,33,40,44,45 . In addition, QEs of M-QTL and QQEs of E-QTL in BC(V) population were much greater than that in RIL(V) population, showing that the backcross population was more susceptible to environment than RIL population. Recently, heterosis study of maize revealed that hybrids had a significant but moderate association between sensitivity to the environment and mean genotype value, whereas inbred lines did not show association 46 . That is why hybrid has superior performance than the parental lines in bad environment. Moreover, the result that QEs (QQEs) of M-QTL and E-QTL at early stages were much greater than that at final stage was found. It indicated that genotype environment interaction was dynamic at various developmental stages and heterosis performance was more susceptible to environment at early stage. Therefore, genotype by environment interaction was important for the stability of hQTL, and it should be taken into account in plant breeding programs 47 .
In present study, we find that the genetic effect of BNP heterosis is dynamic and complicated. It is involved in the different numbers of loci, dominance effect, epistatic interactions and QTL by environmental interactions at different developmental stages.
Totally, four populations were used based on experimental design. The first population was an RIL population. It was developed from the F 1 of 'Xinza 1' 31,49 . The second population was an RILV population from the XZV hybrid. One hundred and eighty RILs were developed through 10 consecutive selfing generations 17,50 . The third population was a backcross developed from RIL population of the XZ hybrid. One hundred and seventy-seven BC hybrids, each hybrid was obtained from a cross where one RIL was used as the female parent and the common parent, 'GX1135' , was used as male parent, respectively 51 . The fourth population was another backcross population (XZV). One hundred and eighty BC hybrids were developed from crosses between RILs from the RILV population used as the female parent and the common parent, GX1135, was used as male parent respectively 17,50 .
Two commercial hybrids of Upland cotton (G. hirsutum) were used as control respectively 17 . The F 1 hybrid 'Ruiza 816' was used as control (CK 1 ) at Handan (E1) and Cangzhou (E2), Hebei Province. The F 1 hybrid "Ezamian 10" was used as control (CK 2 ) at location Xiangyang (E3), Hubei Province. Additionally, two special plots, each consisting of two rows of the XZ hybrid, and its parents 'GX1135' and 'GX100-2' respectively, were used as controls in the experiments of population 1 and population 3. Similar controls were set for the experiments of population 2 and population 4; each plot consisting of the XZV hybrid F 1 , and its parents 'GX1135' and 'VGX100-2' .
Field arrangement and trait evaluation. The four populations and controls were planted at three locations (E1: Handan, E2: Cangzhou, Hebei Province; E3: Xiangyang, Hubei Province) in 2012. In experiment of population 1 and population 2, two-row plots with each plot were used. However in experiment of population 3 and population 4, six-row plots with each plot consisting of two rows of BC hybrid [RIL(V) × GX1135] in the middle, and two rows in both sides for the corresponding parents: one side two rows for the corresponding RIL(V)′ , and the other two rows for 'GX1135' . Each line in RIL(V)′ population was used as the female parent in BC(V) population and was the same one in RIL(V) population. For ease of description, we will refer to the RIL(V) Scientific RepoRts | 6:35515 | DOI: 10.1038/srep35515 s in BC(V) population as RIL(V)′ population, respectively. So in both population 3 and 4, the female were marked as RIL(V)′ , respectively.
The field planting followed a randomized complete block design with two replications at each location. Field management followed the local conventional standard field practices. Boll samples of eight typical plants were measured for boll number per plant (BNP) at four time points, respectively (t1: August 1th, t2: August 15th, t3: August 30th, and t4: September 15th). The data for different stages of t1, t2, t3 and t4 were used for QTL mapping. The incremental values in interval time Δ t1-2, Δ t2-3 and Δ t3-4 were used to map conditional QTL.
DNA isolation, genotype analysis and linkage map construction. Young leaves were collected from labeled parents and two RIL individuals. Extraction of individual genomic DNA and population genotype analysis were carried out following the methods of Liang et al. 48 . A total of 48,836 pairs of SSR primer were used to screen polymorphic loci between three parents. Totally, 653 polymorphic loci for the XZ hybrid and 382 polymorphic loci for the XZV hybrid were acquired and used to conduct genotype analysis of population. and MPH data. A stringent LOD threshold of 3.0 was used to declare suggestive QTL 17,48 . The graphic representation of the linkage group and QTL marked were created by Map Chart 2.2 54 . QTL nomenclature used in rice was employed 55 . At the single-locus level, the genetic effects in BC were defined as follows: a = (P1P1 − P2P2)/2; MPH = d = (BC − (P1P1 + P2P2)/2) and BC = (a + d) (P1 is the recurrent parent). QTL detected only in RIL or BC and not for MPH were referred to as additive. QTLs with d/a ≤ 1 were considered as being complete or partial dominant loci. QTLs with d/a > 1 or only detectable for MPH data were considered as over-dominant loci 10 . Two-locus analysis that tests the main-effect QTL (M-QTL), and digenic epistatic QTL (E-QTL) and their environmental interactions (QTL × environment, QE), was conducted using three data sets by the software ICIMapping 4.0 (http://www.isbreeding.net/). LOD thresholds were set at 2.5 and 5.0 for declaring the presence of M-QTLs and E-QTLs mapping, respectively, same as Shang et al. 17 . Basic statistical analysis was conducted by the software SPSS version 19.0 (SPSS, Chicago, USA).