Genome-Wide SNP Linkage Mapping and QTL Analysis for Fiber Quality and Yield Traits in the Upland Cotton Recombinant Inbred Lines Population

It is of significance to discover genes related to fiber quality and yield traits and tightly linked markers for marker-assisted selection (MAS) in cotton breeding. In this study, 188 F8 recombinant inbred lines (RILs), derived from a intraspecific cross between HS46 and MARCABUCAG8US-1-88 were genotyped by the cotton 63K single nucleotide polymorphism (SNP) assay. Field trials were conducted in Sanya, Hainan Province, during the 2014–2015 cropping seasons under standard conditions. Results revealed significant differences (P < 0.05) among RILs, environments and replications for fiber quality and yield traits. Broad-sense heritabilities of all traits including fiber length, fiber uniformity, micronaire, fiber elongation, fiber strength, boll weight, and lint percentage ranged from 0.26 to 0.66. A 1784.28 cM (centimorgans) linkage map, harboring 2618 polymorphic SNP markers, was constructed, which had 0.68 cM per marker density. Seventy-one quantitative trait locus (QTLs) for fiber quality and yield traits were detected on 21 chromosomes, explaining 4.70∼32.28% phenotypic variance, in which 16 were identified as stable QTLs across two environments. Meanwhile, 12 certain regions were investigated to be involved in the control of one (hotspot) or more (cluster) traits, mainly focused on Chr05, Chr09, Chr10, Chr14, Chr19, and Chr20. Nineteen pairs of epistatic QTLs (e-QTLs) were identified, of which two pairs involved in two additive QTLs. These additive QTLs, e-QTLs, and QTL clusters were tightly linked to SNP markers, which may serve as target regions for map-based cloning, gene discovery, and MAS in cotton breeding.


INTRODUCTION
The Gossypium genus is the most important source of natural textile fiber. It consists of 50 species approximately, including four cultivated species, G. arboreum, G. herbaceum, G. hirsutum, and G. barbadense. G. hirsutum, accounting for 95% of overall cotton production, is characterized by high yield, moderate fiber quality, and wide adaptability (Cai et al., 2014). In the past few years, area under cotton cultivation has declined worldwide, mainly due to high production costs and strong market competition with other crops (Mei et al., 2013). Thus, developing new cotton cultivars with superior fiber quality and high yields are immense to meet the demand of textile industry demand and to maintain the profitability of cotton for its growers.
Fiber quality and yield traits are complex traits, controlled by a multitude of quantitative trait locus (QTLs; Said et al., 2013). There is a complicated genetic correlation between fiber quality and yield due to different population types and parental lines (Qin et al., 2008;Yu et al., 2013;Zhang et al., 2014;Cao et al., 2015;Wang et al., 2015a). Therefore, improving yield and fiber quality, simultaneously, is a long-term challenge for cotton breeders. Traditional breeding procedures are increasingly difficult because of long duration and low selective efficiency (Shen et al., 2005). Marker-assisted selection (MAS) could be one of the more efficient approaches for breeding elite upland cotton cultivars. To date, majority of cotton genetic maps have been developed based on interspecific populations (Reinisch et al., 1994;Lacape et al., 2003;Nguyen et al., 2004;Rong et al., 2004;Guo et al., 2007;Yu et al., 2011bYu et al., , 2013, which had little usage in upland cotton breeding programs (Ulloa et al., 2002). The intraspecific genetic maps of upland cotton have been constructed to detect QTLs for fiber quality and yield traits (Zhang et al., 2003(Zhang et al., , 2009Shen et al., 2006;Wang et al., 2006Wang et al., , 2015bQin et al., 2008;Wu et al., 2008;Sun et al., 2011;Liang et al., 2013;Tan et al., 2014;Islam et al., 2015;Zhang Z. et al., 2015). However, due to low levels of intraspecific DNA marker polymorphisms in upland cotton, most intraspecific maps had a relatively low density and could not satisfy for MAS and mapbased cloning. Therefore, it is necessary to develop new type of markers which can enable mapping upland cotton populations to obtain high polymorphism.
In comparison with other molecular markers, single nucleotide polymorphisms (SNPs) marker which provides the most abundant form of genetic variations, and is characterized by lower mutation rates, higher numbers, and higher accuracy (Ball et al., 2010;Yu et al., 2011a). Previous research has showed that SNPs plays role in phenotypic changes and act as functional marker for traits in MAS, when located in a gene or promoter region (Beales et al., 2005;Konishi et al., 2006). It led to the discovery of superior high-density SNP gene-chip technology was then developed as a superior method for linkage mapping and QTL detection. Now, it is being used extensively to detect QTL in bi-parental populations of many crop species (McNally et al., 2009;Gao et al., 2015). But only a few studies had reported its use in intraspecific populations of upland cotton, either with a limited number (Byers et al., 2012;Yu et al., 2012;Gore et al., 2014;Zhu et al., 2014) or with lower density (Hulse-Kemp et al., 2015;Wang Y. et al., 2015). So, there is a huge knowledge gap to be filled by a comprehensive study.
Herein, a 63 K Illumina SNP assay was used to screen 188 recombinant inbred lines (RILs) derived from the cross of HS46/MARCABUCAG8US-1-88, and a final map with 2618 loci and 0.68 cM high density map was constructed. The aims of this study were to identifying stable QTLs for fiber quality and yield traits and their tightly linked SNP markers for MAS in upland cotton breeding.

Plant Materials
A RIL population of 188 individual lines was developed following a modified single-hill procedure (bulked progeny row; Wu et al., 2008) by crossing two upland cotton cultivars, HS46 and MARCABUCAG8US-1-88. The former is a commercial cultivar with good fiber qualities and higher yield, and the later, a male parent was a germplasm with good resistance. The RILs and their parents were kindly provided by USDA-ARS, Starkville, MS, USA in 1999 (Liu et al., 2012).
In September 2014, the seeds of 188 RILs and two parents were sown in Yacheng and Baogang (two different environments were hereinafter referred to as Yc and Bg, respectively) in Sanya, Hainan Province, China. A completely randomized block design with two replications was applied in each location. Plot size was one row with 7.0 m long and 0.8 m wide. Standard cultivation, weed and insect control practices were followed throughout the growing season.

Phenotypic Measurement and Analysis
A total of 20 normally open bolls were hand-harvested from each line. Approximately 20 g of fiber from each sample was measured by HVI 1000 (Uster R Hvispectrum, Spinlab, USA) under controlled environmental conditions (20 • C and 65% RH) in the Cotton Quality Supervision, Inspection and Testing Center, Ministry of Agriculture, Anyang, Henan province, China. The fiber quality traits include fiber length (FL, mm), fiber length uniformity (FU, %), micronaire (MIC), fiber elongation (FE, %), and fiber strength (FS, cN.tex −1 ). Yield traits consist of boll weight (BW, g) and lint percent (LP, %).
The basic statistics for the phenotypic data of the RILs, the significance of differences for each trait between the two parents, and the correlation among different traits were calculated by SPSS20.0. The variance components for fiber quality and yield traits were estimated by QTModel 1 .

SNP Maker Analysis and Genotyping
Genomic DNA was extracted from young leaves of the 188 RILs and two parents using modified CTAB method (Paterson et al., 1993).
The 188 RILs and their parents were genotyped with cotton 63K SNP array (Hulse-Kemp et al., 2015) from Emei Tongde Technology Development, Co. Ltd (EMTD; Beijing, China 2 ). The array, consisted of 63,058 SNPs, were derived from published literatures (Van Deynze et al., 2009;Byers et al., 2012;Lacape et al., 2012;Rai et al., 2013). Candidate SNPs suitable for further analysis were identified as follows: (1) SNPs were filtered by excluding those with monomorphic markers or with poor quality data; (2) SNPs which the parental genotypes were inconsistent with progeny genotypic ratios or parental genotypes data had missing information were removed from the dataset; (3) SNPs of 188 RILs with missing values more than 40% were removed. Candidate SNPs, obtained from the array, were further aligned to the tetraploid upland cotton (TM-1) reference genome , using BWA software (Li and Durbin, 2009). Only SNPs with less than two mismatches based on a high quality sequence and a mapping Q-Value > 20 were used. The retained SNPs were sent to samtools  and then screened for polymorphism between the mapping parents. Polymorphic SNPs were classified based on Illumina GenTrain score and call frequencies across samples (Hulse-Kemp et al., 2015). Minor allele frequencies of polymorphic markers were only used to genotype 188 RILs.

Map Construction
Linkage maps were constructed by JoinMap 4.0 Version Software (Van Ooijen, 2006), using a regression approach with the log of odds (LODs) score of 3∼10 and the jump threshold of 0.5. Converting recombination frequencies into map distances were calculated using Kosambi's mapping function (Kosambi, 1944).
The chi-square analysis was performed to test segregating markers which deviated from 1:1 expected segregation ratio. A segregation distorted region (SDR) was defined as region with at least three adjacent loci showing significant segregation distortion (P < 0.05; Yu et al., 2011b).
Based on the results of SNPs aligned to the G. hirsutum reference genome by BWA (Li and Durbin, 2009), software CIRCOS 0.69 was used to compare the collinearity of SNPs based on their genetic positions and physical positions.

QTL Analysis
Fiber quality and yield traits related QTLs detection were performed by WinQTLCart2.5 software (Wang et al., 2001), using composite interval mapping (CIM) approach. The LOD threshold of significant QTL was calculated by 1,000 permutation tests with a significance level of P < 0.05, a mapping step of 1.0 cM, and five control markers. LOD score values between 2.5 and permutation test LOD threshold were used to declare suggestive QTL. E-QTLs were detected by IciMapping ver. 4.0 software (Li et al., 2007) using multi-environment trials (METs) function and the inclusive composite interval mapping (ICIM) method. The e-QTLs identification was done with preadjusted IciMapping parameters: Scan = 5 cM, LOD = 5.0, and PIN = 0.0001. A graphical representation of the linkage groups and QTLs was created by Map Chart 2.2 (Voorrips, 2002).
QTLs were named as following: q + trait abbreviation+ chromosome number+ QTL number.

Phenotypic Evaluation of RIL Populations
Descriptive statistics for the fiber quality and yield traits of the RIL population, as well as their parents across two environments were presented in Table 1. FL, FU, BW, FS, and MIC of HS46 showed significantly higher than those of MARCABUCAG8US-1-88 in one or both locations; however, there were no significant differences were found for FE and LP between parental lines. In the RIL population, all analyzed traits presented continuous variation and transgressive segregation, and accorded with normal distributions.
Results showed that the FL had the highest broad-sense heritability among all traits, indicating that it could be mainly controlled by genotype, while remaining six traits had lower broad-sense heritability with values of 0.26 for FU, 0.27 for MIC, 0.37 for FE, 0.28 for FS, 0.31 for BW, and 0.38 for LP, suggesting the environmental effects were important for the performance of these traits ( Table 2).  The correlation analysis for fiber and yield traits based on RILs data over two environments ( Table 3). Results revealed significant negative correlation of LP with FS and FU, and significant positive correlation with MIC and FE. Similarly, BW was significant positively correlated with MIC and FS. Among the five fiber traits, all trait pairs presented significantly correlation except for FU -FE, FU -MIC, and FS -FE.

Map Construction
Among 63,058 SNPs used for screening RIL population, 3120 SNP markers (4.9%) were polymorphic between the two parents. Of those, 2618 were mapped on 26 chromosomes of upland cotton. The total length of this map was 1784.28 cM with average marker density of 0.68 cM ( Table 4; Supplementary  Table S1; and Figures 1-5). There were 101 SNPs on each chromosome averagely, with 1198 SNPs on At subgenome and 1420 SNPs on Dt subgenome, respectively. Uneven distribution of SNP markers on cotton chromosomes was observed. Chr14 had the highest number of SNPs (274 loci), while Chr12 had the lowest SNPs (15 loci). The average chromosome length was 68.63 cM, and the longest chromosome was Chr18 with 119.97 cM and the shortest one was Chr23 with 30.93 cM. The total lengths of At and Dt subgenomes were 888.61 and 902.67 cM, respectively. There were more loci on Dt subgenome than At subgenome. Twenty-two gaps (marker interval > 10 cM) were found on this genetic map, in which, 10 on At subgenome and 12 on Dt subgenome were observed.

Segregation Distortion
Among the 2618 mapped SNPs, 13.29% (348) showed segregation distortion and most loci (71.26%) showed a higher allelic frequency from the female parent (Table 4; Figures 1-5). These SNPs were unevenly distributed on the 26 cotton chromosomes and ranged from 0 to 35 loci on each chromosome ( Table 4). As previous report (Yu et al., 2011b), more distorted loci were located on the Dt subgenome than on the At subgenome (180 versus 168). Segregation distortion was non-random across the linkage map. Three chromosomesChr12, Chr17, and Chr25, showed serious segregation distortions of 33.33, 32.20, and 41.18%, respectively. Furthermore, a total of 35 SDRs were found on 20 chromosomes with 18 SDRs on the At subgenome and 17 SDRs on the Dt subgenome (Table 4; Figures 1-5). Interestingly, the distorted loci in some of the SDRs (SDR03-2, SDR06-1, SDR20-1, and SDR25-2) skewed toward the same allele and showed similar degree of segregation (Figures 1-5).

Collinearity Analysis
All the mapped 2618 SNPs were aligned to the G. hirsutum reference genome to validate the genetic map. Alignments indicated that the genetic map constructed in the present study had good collinearity with the physical map (Figure 6), suggesting the high quality of the RIL map. However, several deviation on Chr04, Chr05, and Chr10 in the At subgenome and Chr16, Chr17, Chr18, Chr21, Chr23, and Chr25 in the Dt subgenome were detected between the genetic map and the   Figure 6).

QTL Analysis of Fiber and Yield Traits
A total of 71 QTLs for fiber quality and yield traits with 4.70∼32.28% of the total explained phenotypic variance (PV) were identified by CIM analysis (Supplementary Table S3; Figures 1-5). In these QTLs, 35 were positive additive indicating HS46 contributed alleles leading to an increase in relevant traits, while 36 were negative additive which meant MARCABUCAG8US-1-88 contributed alleles to increase fiber quality and yield traits. These QTLs were detected on 21 chromosomes, except Chr02, Chr06, Chr08, Chr13, and Chr26. Sixteen QTLs explaining 5.35∼32.28% of the PV were detected in both locations ( Table 5). Among them, three were three QTLs each for FL, FU, MIC, FE, and LP, and one for BW.

Fiber Strength
Eight QTLs, explaining 5.14∼9.96% of the total PV, were detected on Chr05, Chr14, Chr19, and Chr20 (Supplementary Table S3; Figures 1-5). Alleles for increasing FS on Chr05 and Chr19 were contributed by HS46, and positive alleles on Chr14 and Chr20 came from MARCABUCAG8US-1-88. All of these eight QTLs were detected in single environment, indicating that the environmental effects were important for the performance of FS.

(Supplementary
FIGURE 4 | Genetic maps of Chr10/Chr20, Chr11/Chr21, and Chr12/Chr26 homoeologous chromosomes and QTL detection for fiber quality and yield traits in RIL population. All legends are same as described for Figure 1.

QTL Clusters and Hotspots
Quantitative trait locus were not randomly distributed across chromosomes and chromosomal regions. Some QTLs were identified as "cluster" and "hotspot, " where clusters and hotspots were defined to contain multiple QTLs within 20 cM regions, approximately, for different and same traits, respectively (Guo et al., 2007;Rong et al., 2007;Said et al., 2013).
In the current study, there were two QTL clusters on Chr05 which contained three and five QTLs, respectively ( Table 6). The Chr05-cluster-1, which possessed three QTLs, was found at 11∼13 cM for FL, MIC, and LP and Chr05-cluster-2 with five QTLs was located at 40∼55 cM for FL, FU, and FS. The FL hotspot, Chr05-hotspot-1, carrying three QTLs, was located at 40∼53 cM. It should be noted that the position of the hotspot coincided with the second identified cluster.
Chr20 had one cluster Chr20-cluster-1, which was identified at 41∼60 cM and carried three QTLs related to FE and FS ( Table 6).

Identification of the Epistatic and QTL × Environment Interactions Loci
Nineteen e-QTLs were identified by the MET analysis of the multi-environment phenotypic values (Supplementary Table S4; Figure 7). Fifteen of them had both an epistatic main effect and minor epistasis × environment interaction effect while remaining four had the former effect only. The explained PV by each QTL ranged from 3.89 to 6.22%, while that by each QTL × environment interaction only ranged from 0.0 to 0.67%. Among them, the highest number of e-QTLs (9) was detected for FL, three for MIC, two for FE, two for FS, one for BW, and two for LP. No e-QTL was detected for FU. For FL, there were two loci i38186Gh-i11502Gh and i15345Gh-i18849Gh overlapped with two additive QTLs, qFL-Chr10-1 and qFL-Chr14-2, indicating both additive and epistatic value of these two loci had played important role in FL. One pair of interacting marker intervals, i08786Gh-i00558Gh and i09371Gh-i09643Gh, was detected for FL on the same chromosome of Chr19, whereas other interacted loci were located on different chromosomes. In addition, some marker intervals had interactions with other multiple marker intervals to control same or different traits. The marker interval, i22642Gh-i41613Gh on Chr21, had interaction with two marker intervals, i35903Gh-i20966Gh on Chr03 and i44474Gh-i03341Gh on Chr17, for FL. The marker interval i16566Gh-i08941Gh on Chr19 interacted with two marker intervals, i14920Gh-i51624Gb on Chr17 and i08832Gh-i09452Gh on Chr19, for MIC. The marker interval i10502Gh-i36496Gh on Chr04 had interactions with two marker intervals, i08933Gh-i28797Gh on Chr19 and i05035Gh-i22015Gh on Chr14, to control two traits, FL and FS, respectively. The marker interval i32883Gh-i13851Gh on Chr18 had interactions with two marker intervals including i38186Gh-i11502Gh on Chr10 and i02298Gh-i42430Gh on Chr01 pertaining to FL and FE, respectively.

SNP Discovery and Map Construction
With the development of theoretical and applied genetic breeding research, high-density genetic maps are becoming more and more important. SNPs were proved to be the most abundant form of genetic variation, providing a rich source of DNA markers (Agarwal et al., 2008). Recently, SNP arrays for many crops had been developed and utilized in MAS breeding (Fang et al., 2013;Gao et al., 2015;Wang X. et al., 2015). In present work, 2618 polymorphic SNP markers from a 63K SNP assay (Hulse-Kemp et al., 2015) were used to construct a relative high-density genetic map for 188 RILs derived from the combination of HS46/MARCABUCAG8US-1-88.
The total length of the linkage map was 1784.28 cM, far longer than previously map lengths reported with the same cross (Shappley et al., 1998;Wu et al., 2008). The interval of the map was 0.68 cM/marker, increasing the density of the previous linkage map to the thickest linkage map in upland cotton and representing a considerable advance over previously map researches based on RFLP, SSR, AFLP, ITISJ, SRAPetc (Shappley et al., 1998;Ulloa et al., 2002;Shen et al., 2005Shen et al., , 2006Zhang et al., 2005Zhang et al., , 2009Wu et al., 2008;Tan et al., 2014;Liu et al., 2015). The collinearity analysis showed that the constructed map had good collinearity with the G. hirsutum reference genome, indicating the high quality and accuracy of the map. One main reason for the lower coverage of At subgenome than Dt subgenome was that two linkage groups, Chr12 and Chr10, only represent 10.91 and 70.34% of the corresponding chromosomes, respectively. Consistent with previous studies (Tan et al., 2014;Cao et al., 2015;Hulse-Kemp et al., 2015), more markers were found on Dt subgenome (54.24%) than At subgenome (45.76%) in the present map, which was attributed to the lower level of polymorphism in At subgenome of upland cotton. The Dt subgenome was longer than At subgenome in this map, duo to Dt subgenome with more loci experienced a higher frequency of recombination (Guo et al., 2007). Although the average density was high, there were still some gaps in several chromosomes. Moreover, like previous maps (Shen et al., 2006;Tan et al., 2014;Cao et al., 2015), although the markers of our map distributed evenly on the entire genome, there were still some chromosomes anchored more markers than others, which might be attributed to the non-random distribution of markers and the lack of marker polymorphism between mapping parents on some chromosomes because of relative narrow range of genetic diversity in upland cotton.
Segregation distortion, regarded as the important source of plant evolution, widely exists in plant populations. The species, population types, crosses, and marker types of plants will lead to significant variance in segregation distortion with different origin, genetic effects, and degree (Xu et al., 1997). The reports from intraspecific population of upland cotton indicated that the ratio of segregation distortion enhanced with the increase of divergence level of the parents (Shappley et al., 1998;Ulloa et al., 2002;Shen et al., 2005Shen et al., , 2006, and the ratio of segregation distortion in the RIL population was higher than that in the F 2 population, which might mainly result from genetic drift (Shen et al., 2006). In present study, 13.29% of the total loci showed segregation distortion (P < 0.05), which was similar to previous

Significance and Potential Application of QTL Mapping
The genotypic value for HS46 was greater (P < 0.05) than that for MARCABUCAG8US-1-88 with respect to FL, FU, BW, FS, and MIC. Similar to the same cross used by Wu et al. (2008), no significant difference between the two parents was detected for FE and LP. Even though the two parents were phenotypically similar regarding these two traits, due to genetic dissimilarities between the two parents, significant differences in this RIL population existed. In the present research, FS had a low heritability and easily affected by environment. Similar results could be found in previous literatures about a fourway cross population in upland cotton (Qin et al., 2008). The obvious difference indicated FS was highly influenced by the experimental environment and difficult in genetic improvement. Wang et al. (2015a) suggested FS was a moderate heritability trait and that all the fiber quality and yield component traits presented significant environmental effects. Moreover, in the research based on same population, FS did not show a higher heritability in all the measured traits but a significant VGE/VP instead (Wu et al., 2008). Therefore, the heritability of FS was not stable and affected by different population and environment. Similarly, BW was highly influenced by environments. Furthermore, BW was significantly and positively correlated with MIC and FS. Consequently, it is feasible to improve BW by selecting these correlated traits, having more accurate repeatability across environments in breeding. Stable QTLs such as qFL-Chr14-3, qFL-Chr15-1, qFU-Chr09-1, qFE-Chr14-1 qFE-Chr20-1, qBW-Chr10-1, qLP-Chr10-1, and qLP-Chr12-1 could be used in breeding. It was more likely that these QTLs could be used to identify candidate genes for these related traits because of the availability of high-density SNP markers. In such case, they will be the potential candidates for fine mapping and ultimate candidate gene discovery.
Although there were many differences in parental lines, mapping populations, and markers type, our results were comparable with earlier identified QTLs. Searching for QTLs of upland cotton in a CottonQTLdb database 3 developed by Said et al. (2015) and statistical data on cotton QTL previously, there were seven QTLs for FL, two QTLs for FU, two QTLs for MIC, four QTLs for FE, two QTLs for FS, and five QTLs for LP in present work, sharing same genetic position (spacing distance < 5 cM) and physical position with earlier reports (Wu et al., 2008;Tan et al., 2014;Liu et al., 2015;Zhang Z. et al., 2015). In addition, five QTLs for FL, six QTLs for FU, eight QTLs for MIC, four QTLs for FE, seven QTLs for FS, six QTLs for BW, and seven QTLs for LP had been mapped on the same chromosomes but not on the same position as reported in the previous researches. These inconsistencies might be due to different genetic backgrounds and DNA markers (Shappley et al., 1998;Wu et al., 2008;Said et al., 2015). Remaining nine additive QTLs (qFE-Chr17-1, qBW-Chr09-1, qBW-Chr09-2, qBW-Chr23-1, qBW-Chr24-1, qBW-Chr25-1, qLP-Chr20-1, qLP-Chr21-1, and qLP-Chr21-2) might be novel loci, due to unavailability of any reports for these traits on those chromosomes. As the number of identified QTLs for fiber quality and yield traits increased, the genetic control of fiber quality and yield will be better understood.

Molecular Mechanism of Trait Correlation and Linkage Drag
Co-localization of QTLs on chromosomes, referred to as "QTL cluster/hotspot, " has previously been reported in cotton (Shappley et al., 1998;Qin et al., 2008;Said et al., 2013) and many other species (Gonzalez et al., 2015;Li et al., 2015). In the present study, 12 certain genomic regions, especially, Chr05, Chr09, Chr10, Chr14, Chr19, and Chr20, were investigated for their involvement in controlling one (hotspot) or more (clusters) fiber quality or yield traits, the similar result was also reported in the publications (Qin et al., 2008;Sun et al., 2011;Said et al., 2013;Yu et al., 2013). The existence of QTL clusters explained why so many traits were highly interrelated.
Based on the comprehensive analysis of clusters and hotspots in this study, breeding programs targeting fiber quality or yield traits can focus on hotspot clustering regions and select around the region. Notably, almost all the hotspots overlapped QTL clusters. The presence of QTL clusters and hotspots proved that genes related to certain traits were more heavily concentrated in certain regions of genome than others (Said et al., 2013). The discovery of cluster and hotspot may be useful in MAS breeding program and may help breeders to select the traits of interests and find novel QTLs once the markers have anchored these regions.

Interaction between Loci within and Across Chromosomes
It was successful to identify e-QTLs with both additive and epistatic effects, e-QTL pairs and epistasis × environment interactions in the present work, which were often neglected in some complex trait studies. Generally, if the proportion of PV explained by the identified additive QTL is close to broad-sense heritability, epistasis is less important (Li et al., 2007). However, for FL, the total PV explained by additive QTLs were much lower than the broad-sense heritability (66%, Table 2), indicating that there were epistatic interactions in these loci. Finally, nine pairs of e-QTLs for FL were detected in our work, of which there were two loci, i38186Gh-i11502Gh on Chr10 and i15345Gh-i18849Gh on Chr14, locating on the same position with two additive QTLs (qFL-Chr10-1 and qFL-Chr14-2), suggesting that both addictive and epistatic effects had played important roles in genetic control of FL. In fact, interactions among loci or QTL × environmental factors contributed a substantial effect to complex trait phenotypic variation (Carlborg and Haley, 2004). Several novel QTLs and specific trait relationships between loci, within and across chromosomes, could be considered as the interactions between loci (Grosse-Brinkhaus et al., 2010).

CONCLUSION
A high-density linkage map was constructed in the upland cotton RIL population using the 63K cotton SNP array. Nine novel QTLs, seven pleiotropic QTL clusters, five hotspots, and 19 e-QTLs for fiber quality and yield traits were identified with tightly linked SNP markers. These QTLs could serve as target regions for map-based gene cloning and MAS in cotton breeding.

AUTHOR CONTRIBUTIONS
CoL and SZ designed and conducted the experiments, analyzed data and wrote the manuscript. CoL and YD performed phenotyping and data analysis, CoL, TZ, and LL participated in field trials. ChL, EY, LM, QH, and MD prepared and reviewed the manuscript. SZ and JC designed and supervised the experiments. All authors have read and approved the final manuscript.

ACKNOWLEDGMENTS
We are grateful to Dr. XU (Zhejiang University, China), Li Chen (Zhejiang University, China), and Zhen Zhang (Chinese Academy of Agricultural Sciences, China) for technical assistance.