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.

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^®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¹.

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²). 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 (Zhang T. et al., 2015), 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 (Li et al., 2009) 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.

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.

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 pre-adjusted 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.

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.

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.

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).

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 physical map. The Dt subgenome showed good coverage of the physical map, representing 96.81% of the genome assembly length, while the At subgenome showed a lower coverage of 88.33% (Supplementary Table S2; Figure 6).

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.

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.

Chr09 contained one cluster and one hotspot (Table 6). The Chr09-cluster-1, identified at 46∼62 cM, carried four QTLs for FU and BW while Chr09-hotspot-1, identified at 47∼62 cM, carried three QTLs of FU and overlapped the cluster Chr09-cluster-1.

Chr10 contained one cluster (Chr10-cluster-1) pertaining to FL, BW, and LP QTLs (Table 6). This cluster was located at 32∼45 cM, containing four QTLs.

Chr14 contained one cluster and three hotspots (Table 6). The Chr14-cluster-1 was located at 6∼26 cM and carried eight QTLs for FL, FS, MIC, and FE. The two hotspots, Chr14-hotspot-1 and Chr14-hotspot-2, identified fell into the range of the Chr14-cluster-1, and were FL and FS related hotspots, respectively. Furthermore, the three hotspots Chr14-hotspot-1, Chr14-hotspot-2, and Chr14-hotspot-3 for FE, FS, and FL overlapped at 1∼17, 6∼22, and 14∼26 cM, respectively, indicating a possible fiber quality hotspot cluster in the region.

Chr19 contained one cluster, Chr19-cluster-1, which ranged from 20∼28 cM, and carried four QTLs for FL, FU, and FS (Table 6).

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).

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.

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., 2005, 2006; Zhang et al., 2005, 2009; Wu 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., 2005, 2006), and the ratio of segregation distortion in the RIL population was higher than that in the F₂ 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 maps of upland cotton (16.6∼36.7%; Tan et al., 2014; Cao et al., 2015; Liu et al., 2015). Most distorted loci skewed to HS46 alleles, which might be the result of chromosomal elements markers diverged via female drive in the meiosis (Fishman and Willis, 2005). Furthermore, most segregation distortion loci occurred in the clusters, like previously reported for interspecific and intraspecific populations (Guo et al., 2007; Zhang et al., 2009; Yu et al., 2011b; Tan et al., 2014; Liu et al., 2015). The loci clustering on the same chromosome or within the same SDR skewed to the same allele suggested that genetic hitchhiking effects existed in upland cotton (Yu et al., 2011b).

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 four-way 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³ 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.

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.

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).