Upland cotton (Gossypium hirsutum L.) is one of the most important economic crops in the world. One of the main objectives for the cotton-producing countries is to improve fiber production. Vegetative growth is the basis for life circle of the plants,which provides the essential nutrition for initiating cotton flowers,bolls,and fiber. Revealing genetic basis and estimating effects for quantitative trait loci (QTLs) related to plant morphological traits^[1] by linkage mapping can provide useful information for conducting marker-assisted selection (MAS).

Association mapping is an efficient method for detecting polymorphic genotypes to phenotypic variations,and the genome-wide association study (GWAS) is a special analysis method for association mapping. It can scan the molecular markers of entire genome for particular traits of interest^{[2, 3]}. Mapping genetic loci can be performed by using RAPD (random amplification polymorphism DNA),SCAR (sequence characterized amplified region),SSR (simple sequence repeat),and SNP (single nucleotide polymorphism) markers,etc. In this study,we utilized PCR-based molecular markers SSRs,which came mostly from the cotton marker database (CMD).

Linkage mapping has been widely used in cotton among the previous studies. But most of them concentrated on the fiber traits and yield traits,such as micronaire value,lint yield,and boll mass^{[4, 5, 6, 7, 8]}. These traits were direct factors related to cotton industry,so the breeders wanted to understand their genetic basis and find ways to improve them. Few works had been reported on plant morphological traits^{[9, 10, 11]}. Indirect factors including plant height,stem diameter,and stem node,also influenced the growth of cotton and hence the yield. Therefore,we aimed to dissect the genetic basis,find the QTSs (quantitative trait SSRs) associated to the traits and provide more information for cotton breeding.

In this study,we performed association mapping for 872 marker loci of 39 cultivars and their 178 F₁ hybrids in three environments. A newly developed software QTXNetwork (http://ibi.zju.edu.cn/software/QTXNetwork/) based on mixed linear model approaches^[12] was applied in this study. A full genetic model with genetic effects of additive,dominance,epistasis and their environment interactions was built. The association mapping results revealed the importance of genetic architecture on plant morphological traits in cotton.

1 Materials and methods 1.1 Parent selection and experimental design

In this study,the mapping population was provided by the Institute of Cotton Research,Chinese Academy of Agricultural Sciences and Tarimu University. These cultivars were representatives of cotton-growing regions in China. Thirty-nine cotton cultivars were sampled and their 178 F₁ hybrids were analyzed for four morphological traits at three environments (e₁: Anyang of Henan Province in 2012; e₂: Anyang of Henan Province in 2013; e₃: Alar of Xinjiang Autonomous Region in 2012). Anyang of Henan Province located in the Yellow River valley cotton-growing region,and Alar of Xinjiang Autonomous Region located in the northwestward cotton-growing region. For field experiments,a completely randomized block design with three replications was applied. Each block was settled with two rows,and each row was kept in a plot of 5.00 m in length and 0.70 m in width. The distance between plants was 0.25 m in width in Anyang,and 0.10 m in width in Alar. Plant height (PLH),stem diameter (STD),stem first (STF) and stem node (STN) were evaluated to dissect the genetic network of these morphological traits.

1.2 DNA extraction and SSR markers screening

Two cotyledons from two plants were collected randomly for each line in the plots in March 2012. DNA from 39 parent samples was extracted according to the method of PATERSON et al.^[13]. The sequences of SSR primers,which came mostly from the cotton marker database (CMD,http://www.cottonmarker.org/cgi-bin/panel.cgi),were downloaded to obtain specific information about the markers. This database is an integrated web-enabled database,providing a centralized access to all publicly available cotton microsatellites and single nucleotide polymorphisms,which provides more efficient applications of molecular markers and facilitates basic and applied research in molecular breeding and genetic mapping in cotton. A standard screening panel contained the polymorphic information of upland cotton cultivars and other tetraploid species. Among the 39 lines,there were 5 052 SSR primer pairs in total to be examined for screening polymorphisms. Those markers included DPL,GH,HAU,JESPR,NAU,BNL,CIR,CGR,DC,MUSS,MGHES,C,CER and CM series. The marker nomenclature consisted of a letter that indicated the origin of the marker,followed by the primer number. SSR analysis was carried out according to the procedure proposed by YAO et al.^[14]. To sum up,872 marker loci were obtained to conduct association analysis.

1.3 Statistical analysis

To analyze different effect values,a full genetic model was used for the phenotypic value of the k-th individual in the h-th environment (y_hk),which can be presented by the following mixed linear model,

$\begin{gathered} {y_{hk}} = \mu + \sum\limits_i {{a_i}x{A_{ik}} + } \sum\limits_i {{d_i}x{D_{ik}} + } \sum\limits_{i{\text{ < }}j} {a{a_{ij}}xA{A_{ijk}} + } \hfill \\ \sum\limits_{i{\text{ < }}j} {a{d_{ij}}xA{D_{ijk}} + } \sum\limits_{i{\text{ < }}j} {d{a_{ij}}xD{A_{ijk}} + } \hfill \\ \sum\limits_{i{\text{ < }}j} {d{d_{ij}}xD{D_{ijk}} + } {e_h} + \sum\limits_i {a{e_{ih}}uA{E_{ihk}} + } \hfill \\ \sum\limits_i {d{e_{ih}}uD{E_{ihk}} + } \sum\limits_{i{\text{ < }}j} {aa{e_{ijh}}xAA{E_{ijhk}} + } \hfill \\ \sum\limits_{i{\text{ < }}j} {ad{e_{ijh}}uAD{E_{ijhk}} + } \sum\limits_{i{\text{ < }}j} {da{e_{ijh}}uDA{E_{ijhk}} + } \hfill \\ \sum\limits_{i{\text{ < }}j} {ad{e_{ijh}}uDD{E_{ijhk}} + } {\varepsilon _{hk}} \hfill \\ \end{gathered} $

(Equation1)

where μ is the population mean; a_i is the additive effect of the i-th locus with coefficient x_Aik; d_i is the dominance effect of the i-th locus with coefficient x_Dik; aa_ij,ad_ij,da_ij and dd_ij are the digenic epistasis effects with coefficients x_AAijk,x_ADijk,x_DAijk and x_DDijk; e_h is the random effect of the h-th environment; ae_ih is the additive × environment interaction effect of the i-th locus in the h-th environment with coefficient u_AEihk; de_ih is the dominance × environment interaction effect of the i-th locus in the h-th environment with coefficient u_DEihk; aae_ijh,ade_ijh,dae_ijh and dde_ijh are the digenic epistasis × environment interaction effects with coefficients u_AAEijhk,u_ADEijhk,u_DAEijhk and u_DDEijhk; and ε_hk is the residual effect of the k-th individual in the h-th environment population.

The mixed linear model can be expressed in matrix notation as follows,

$\begin{gathered} y = 1\mu + {X_A}{e_A} + {X_D}{e_D} + {X_{AA}}{e_{AA}} + {X_{AD}}{e_{AD}} + \hfill \\ {X_{DA}}{e_{DA}} + {X_{DD}}{e_{DD}} + {U_E}{e_E} + {U_{AE}}{e_{AE}} + \hfill \\ {U_{DE}}{e_{DE}} + {U_{AAE}}{e_{AAE}} + {U_{ADE}}{e_{ADE}} + {U_{DAE}}{e_{DAE}} + \hfill \\ {U_{DDE}}{e_{DDE}} + {e_\varepsilon } \hfill \\ = Xb + \sum\limits_{v = 1}^7 {{U_v}{e_v} + {e_\varepsilon }~{\text{MVN}}} \left( {Xb,\sum\limits_{\nu = 1}^7 {\sigma _v^2{U_v}U_v^T + I\sigma _\varepsilon ^2} } \right) \hfill \\ \end{gathered}$

(Equation2)

where y is an n×1 column vector of phenotypic values and n is the sample size of observations; μ is a vector of population mean; X is the known incidence matrix related to the fixed effects; b is a vector of all fixed genetic effects in the experiment; U_v is the known coefficient matrix for the v-th vector of random environment-specific genetic effects e_v ~ MVN(0,Iσ_v²); and e_ε ~ MVN(0,Iσ_ε²) is an column vector of residual effects.

The total heritability can be estimated by the following equation,

$\begin{gathered} h_T^2 = h_G^2 + h_{GE}^2 \hfill \\ = \left( {\sum\limits_i {h_a^2} + \sum\limits_j {h_d^2} + \sum\limits_{i{\text{ < }}j} {h_{aa}^2} + \sum\limits_{i{\text{ < }}j} {h_{ad}^2} + \sum\limits_{i{\text{ < }}j} {h_{da}^2} + \sum\limits_{i{\text{ < }}j} {h_{dd}^2} } \right) + \hfill \\ \left( {\sum\limits_i {h_{ae}^2} + \sum\limits_j {h_{de}^2} + \sum\limits_{i{\text{ < }}j} {h_{aae}^2} + \sum\limits_{i{\text{ < }}j} {h_{ade}^2} + \sum\limits_{i{\text{ < }}j} {h_{dae}^2} + \sum\limits_{i{\text{ < }}j} {h_{dde}^2} } \right) \hfill \\ = \left( {h_A^2 + h_D^2 + h_{AA}^2 + h_{AD}^2 + h_{DA}^2 + h_{DD}^2} \right) + \hfill \\ \left( {h_{AE}^2 + h_{DE}^2 + h_{AAE}^2 + h_{ADE}^2 + h_{DAE}^2 + h_{DDE}^2} \right) = \hfill \\ \left( {h_A^2 + h_D^2 + h_I^2} \right) + \left( {h_{AE}^2 + h_{DE}^2 + h_{IE}^2} \right) \hfill \\ \end{gathered}$

(Equation3)

where h_T² is the total heritability; h_G² is the heritability for genetic effects; h_GE² is the heritability for genetic-environment interaction effects; h_A² is the heritability for additive effects; h_D² is the heritability for dominance effects; h_I² is the heritability for epistasis effects including AA,AD,DA,DD; h_AE² is the heritability for environment-specific additive interaction effects; h_DE² is the heritability for environment-specific dominance effects; h_IE² is the heritability for environment-specific epistasis effects including AAE,ADE,DAE,DDE.

1.4 Association analysis

The marker-trait association analyses were carried out by the above mixed linear model. All the detected marker-trait association loci were fitted by a full genetic model to predict the additive effects,dominance effects,epistasis effects,and their environment interaction effects. The mixed linear model framework with Henderson method Ⅲ^[15] was used to construct the F-statistic test to implement the association analysis. Permutation test was set for a total of 2 000 times to calculate the critical F-value to control the experiment-wise type Ⅰ error (α_EW < 0.05). The QTS effects were predicted by using the MCMC (Markov Chain Monte Carlo) algorithm with 20 000 Gibbs sample iterations^[16]. Based on the information of genetic effects of detected-QTSs for the four plant morphological traits,the breeding values ($\hat u + \hat G + \hat GE$) for the superior lines and superior hybrids of the mapping population were predicted^[17].

2 Results 2.1 Estimated heritability and genetic effects

A total of 25 experiment-wise significant QTSs (P_EW-value < 0.05) were detected for the four plant morphological traits (5-7 QTSs for each trait) in cotton (Fig. 1,Table S1 available at http://www.journals.zju.edu.cn/agr/EN/article/showSupportInfo.do?id=10359). Estimated heritability was quite high (63.08%-78.28%),and mainly contributed by genetic main effects (h_G²= 48.09% for PLH,36.09% for STD,40.18% for STF,and 56.72% for STN) (Table 1). It was suggested that plant morphological traits were quite stable across environments. The additive effects were the primary genetic component for STD and STF,while dominance effects were dominated for PLH and STN.

2.2 Plant height trait

Plant height (PLH) is one of the most common traits among the plant morphological traits,which plays an important role in formation of plant morphology type as a key trait. Cotton plant height is also considered to be related with yield and fiber quality traits. In this study,it was discovered that dominance and its environment interaction effects were the main components of total heritability (h_D²≈28.97%,h_DE²≈20.97%). No epitasis or their environment interaction was detected.

A total of seven QTSs were detected for PLH (Fig. 1,Table S1 available at http://www.journals.zju.edu.cn/agr/EN/article/showSupportInfo.do?id=10359). There were seven individual QTSs highly significant (P_EW-value < 1×10^-5) (Table 2),among which three QTSs (GH111-5,GH638-3 and JESPR153-6) were found with additive effects. And the other four QTSs were detected with dominance effects. The GH638-3 was detected with environment-specific dominance effects in e₁,while NAU1269-3 was discovered with environment-specific additive effects in e₁. No pair epistasis interaction effect was identified. Selecting genotypes qq of GH111-5 and QQ of GH638-3,JESPR153-6 for different environments could improve PLH. Due to large environment-specific dominance effect,heterozygote Qq of GH638-3 in e₁ could increase PLH. Besides,selecting genotype qq of NAU1269-3 in e₁ could also result in higher PLH.

2.3 Stem traits

For STD,a total of five QTSs were detected (Fig. 1,Table S1 available at http://www.journals.zju.edu.cn/agr/EN/article/showSupportInfo.do?id=10359). Four highly significant QTSs (P_EW-value < 1×10^-5) were listed in Table 2. All the four highly significant QTSs were detected for additive effects. Environment-specific effects were identified for CGR5534-3 by de₁ and DPL282-1 by ae₁ and ae₃. Increasing thick diameter of stem could be expected by selecting genotypes QQ of BNL3569-6 and qq of MGHES44-3. Selecting genotypes qq of DPL282-1 in e₁ and Qq of CGR5534-3 in e₁ could also improve this trait due to environment-specific additive/dominance effects.

For STF,there were six QTSs (CGR5704-1,DPL268-1,DPL0354-1,HAU0377-5,HAU2119-1 and JESPR153-5) highly significant with positive additive effects (Fig. 1,Table 2). For this trait,additive effects contributed nearly half heritability. In addition,CGR5704-1,HAU2119-1 and JESPR153-5 had environment-specific additive effects in e₁ and e₃,and the same was DPL0354-1 in all three environments. Besides,DPL268-1 was detected with positive effects. Among these six QTSs,all the effects were negative in e₁ but positive in e₃; therefore STF could be reduced via selecting genotypes QQ of CGR5704-1,DPL0354-1,HAU2119-1 and JESPR153-5 in e₁ or qq of CGR5704-1 in e₃.

The number of STN was determined by genetic and environmental factors,and had an essential impact on the final plant height. A total of eight QTSs were significantly dissected for STN (Fig. 1,Table S1 available at http://www.journals.zju.edu.cn/agr/EN/article/showSupportInfo.do?id=10359). Most genetic effects were contributed by additive and dominance effects,ranging from 22.18% to 26.47% (Table 1). For the highly significant QTSs (Table 2),CGR6270-2 had negative dominance effect,while other six QTSs (CGR6150-1,DPL0504-4,GH189-1,HAU1523-1,NAU862-3 and NAU1221-5) were additive effects. NAU862-3 was discovered with environment-specific additive effect in e₁. The only one epistasis was between GH189-1 and NAU1221-5,and environment-specific dominance-additive effect in e₁ was discovered (Table S1 available at http://www.journals.zju.edu.cn/agr/EN/article/showSupportInfo.do?id=10359),which/ had the highest heritability in STN. Selecting genotypes qq of DPL0504-4,NAU862-3 or QQ of CGR6150-1,GH189-1,HAU1523-1 and NAU1221-5 could increase STN. Genotype QQ of NAU862-3 in e₁ could also be chosen to further increase STN.

2.4 Predicted genetic effects for different genotypes

Based on the estimated genetic effects of each QTS for three environments,the maximum and minimum genotypic effects of the superior lines and superior hybrids were predicted (Table 3). Moreover,the genotypic effects of homozygotes (QQ,qq) and heterozygote (Qq) were also predicted for the four morphological traits in the three environments. For the traits of PLH,STD,STF and STN,the predicted genotypic effects were positive for major-allele homozygote QQ,but negative for minor-allele homozygote qq. It was suggested that the selection had the potential to increase plant height,stem diameter,first stem,and stem node. The predicted genotypic effects of heterozygote Qq were positive for PLH and STN in all environments,but various for STF and STD across the environments,suggesting that the PLH and STN could be improved through hybrid vigor.

These four morphological traits are moderate,and need to be balanced through breeding and cultivation techniques. Except for the superior lines (±) of STF,most of the absolute breeding values of superior lines (±) in the mapping population were larger than those of homozygote genotypes QQ and qq. It was indicated that the breeders could increase or decrease three morphological traits (PLH,STN and STD),through selecting superior lines with a mixture of homozygote QQ and qq genotypes of detected QTSs. And the breeding value of the superior hybrid could also be reinforced for each trait by selecting a hybrid with a mixture of heterozygote (Qq) and homozygotes (QQ,qq).

3 Discussion

Association mapping is a powerful tool to dissect quantitative trait loci underlying complex traits,and has been widely used in different studies. But due to the lack of suitable statistical methodology and heavy computational demand,epistasis and environment interaction effects were often missed in current GWAS^{[18, 19, 20]}. Ignoring epistasis and environment interaction might lead to biased estimation of effects on loci and lower detecting power^[21]. Furthermore,the heritability of target traits would be missed if reduced models were used^{[22, 23]}. To solve the computational burden problems,LOU et al.^[24] developed a new approach,which combined GPU (graphics processing unit) with GMDR (generalized multifactor dimensionality reduction). GPU could accelerate the operating speed in order to deal with millions of markers,and GMDR was able to detect high dimension interactions.

In this study,the software QTXNetwork was utilized to perform association mapping^[12]. QTXNetwork is a software for GPU computation of linkage and association analyses on full genetic model,including epistasis and environment interaction of complex traits based on omics genotypes. Compared with the reduced models used in the previous research,a full mixed linear model (MLM) could dissect total heritability into different parts: additive,dominance,epistasis and their environment interactions. Although diverse softwares such as PLINK^[25],TASSEL^[26],R^[27] were used,and different models including MLM,GLM were applied in previous studies,only the narrow or broad-sense heritability could be calculated^{[28, 29, 30, 31]}.

Several highly significant QTSs detected by this analysis for morphological traits were also found by other mapping analyses in cotton,which were DPL0504^[32],CGR5534^[33],NAU1269^[34],HAU2119^[35],JESPR153^[36]. This could validate the efficiency of the mixed linear model approach for association mapping to identify relevant loci for complex traits. Furthermore,we identified different alleles of one SSR that had effects on STF (JESPR153-5) and PLH (JESPR153-6). This provided evidence that some QTSs could have pleiotropic or physiological associations for accelerating the process of MAS programs in cotton breeding. We had also calculated the R-values (correlation coefficient) of the predicted genotypic effects in Table 1. For these four traits,the R-values ranged from 0.79 to 0.93,which were close to the total heritability. It was suggested that the genotypic variation could be predicted very well by the full genetic model.