A diversity panel consisting of 355 representative upland cotton accessions (Supplementary Table S1) obtained from cotton germplasm collections in our laboratory and from the low-temperature germplasm genebank of the Cotton Research Institute, Chinese Academy of Agricultural Sciences (CRI-CAAS), was examined. The population consisted of 331 cultivars or accessions developed in China and 24 introduced from abroad. These accessions were divided into the following five groups according to their ecological areas: (1) YR group (162 accessions from the Yellow River region in China), (2) YZR group (51 accessions from the Yangtze River region in China), (3) NW group (98 accessions from the northwest inland region in China), (4) LN group (20 accessions from the Liaoning province in China) and (5) foreign group (20 accessions from the United States of America and 4 accessions from central Asia countries).

All of the 355 upland cotton accessions were planted at Anyang, Henan, China (36⋅ 08′ N, 114⋅ 48′ E) in 2014 and 2015 (designated AY-14 and AY-15, respectively) and at Shihezi (SHZ), Xinjiang, China (44⋅ 31′ N, 86⋅ 01′ E) in 2014 and 2015 (designated SHZ-14 and SHZ-15, respectively). The field experiments followed a randomized complete block design with three replications. At Anyang, Henan, each accession was grown in a single-row plot (5.00 m long and 0.80 m row wide) with 18–23 plants, while at Shihezi (SHZ), Xinjiang, each accession was grown in a plot with 30–40 plants in two rows, with 0.10 m between the plants in each row and 0.45 m between the rows. The trial management followed standard breeding field practices. The fully opened random 20 cotton middle bolls for each accession were manually harvested annually in September. The lint percentage (LP, %) of each harvested sample was evaluated. The correlation coefficients for the lint percentage between environments were calculated and the analysis of variance (ANOVA) was conducted using R software. The broad-sense heritability of lint percentage was calculated using the R software package “lem4.”

Genomic DNA from all of the accessions was extracted from young leaf tissue using a modified cetyl trimethyl ammonium bromide (CTAB) method (Paterson et al., 1993). SNP genotyping of the association panel was performed using a SLAF-seq approach (Sun et al., 2013). Two restriction enzymes (Rsa I and Hae III, New England Biolabs, NEB, USA) were used for library preparation. Paired-end sequencing (125 bp at each end) was performed on an Illumina HiSeq 2500 system (Illumina, Inc., San Diego, CA, USA) according to the manufacturer's recommendations. The GATK and Samtools packages were used for SNP calling. BWA software was used to map the raw paired-end reads onto the reference genome (Gossypium hirsutum v 1.0) (Li et al., 2015). SNPs with a missing rate < 10% and a minor allele frequency (MAF) ≤ 0.05 were excluded. A total of 81,675 SNP markers were identified and used for the subsequent analysis.

The best linear unbiased predictions (BLUPs) for lint percentage across the four environments were estimated using R software. The BLUP values for the four environments and the phenotypic values for each environment's lint percentage were used for GWAS. A mixed linear model (MLM) was used to calculate the associations in all of the analyses by incorporating both principal components (PCs) and kinship data (Lipka et al., 2012). The suggestive and significant P thresholds were 6.12E–06 (P = 0.5/n, where n = the number of markers used; −log₁₀ (0.5/81,675) = 5.21) and 6.12E–07 (P = 0.05/n, where n = the number of markers used; −log₁₀ (0.05/81,675) = 6.21), respectively, for the entire population (Holm, 1979; Mao et al., 2015). Manhattan plots were performed using the R software package “CMplot.”

The phenotypic value of each haplotype was estimated by averaging the phenotypic values over the accessions with each type of SNP locus associated with a specific target trait. The favorable haplotypes (FHs) were subsequently identified according to the breeding objectives of target trait. Box plots for the relative phenotypic values were performed using R software.

Linkage disequilibrium (LD) blocks containing SNP loci associated with target traits were generated using the R software package “LDheatmap” (Shin et al., 2013). LD blocks harboring significantly associated SNPs were then defined as the candidate genome regions. The genes distributed in these regions were collected. The gene annotation closest to the Arabidopsis homolog was obtained from the Arabidopsis information resource (http://www.arabidopsis.org/). Transcriptome sequencing data (SRA: PRJNA248163) from 10 upland cotton (TM-1) tissues (including root, stem, leaf, torus, seed, cotyledon, fiber-5 DPA (days post-anthesis), fiber-10 DPA, fiber-20 DPA and fiber-25 DPA) were available on the EBI website (http://www.ebi.ac.uk/). Expression analysis of the high-throughput mRNA sequencing was performed using TopHat and Cufflinks software (Trapnell et al., 2012). Heat maps of the putative candidate gene expression patterns were created using the R package “pheatmap.”

Significant variation was observed for lint percentage among the 355 upland cotton accessions, ranging from 26.58 to 48.13%, with an average of 41.03%. In the AY-14, AY-15, SHZ-14, and SHZ-15 environments, the natural population exhibited average lint percentage values (±SD) of 41.05 ± 2.87, 39.70 ± 3.18, 41.78 ± 2.78, and 41.58 ± 2.64%, respectively (Figure 1). The ANOVA showed that the genotype (G), environment (E) and genotype × environment (G × E) interactions have significant effects on the lint percentage (P < 0.001, Supplementary Table S2). The phenotypic trait also showed significant differences (P < 0.01) among the four environments (Supplementary Table S2), suggesting that the lint percentage was significantly influenced by the environment. However, significant positive correlations (P < 0.001) between the genotypes across the four environments were observed, with correlation coefficients ranging from 0.7462 to 0.8620 (Supplementary Table S3). Additionally, the broad-sense heritability of lint percentage was as high as 69.72% (Supplementary Table S2). The results indicated that the stability of the lint percentage was high, even though significant G × E was observed.

In our previous studies, 355 upland cotton accessions were grouped into two major subpopulations by principal component analysis (PCA) and clustering analysis, and LD was estimated as r² (the squared Pearson correlation coefficient) between all pairs of SNP markers; the approximated LD decay distance was 100 kb (Su et al., 2016). In this study, GWASs were conducted for lint percentage using the BLUPs across four environments and an individual environment in an MLM. The threshold of −log₁₀ (P) > 5.21 was also derived from the quantile-quantile (QQ) plots because most of the upward deviation from the linear line occurred around −log₁₀ (P) = 5.21, which presumably indicates true positives (Figure 2B). Twelve SNPs were significantly associated with lint percentage; these SNP loci were positioned on chromosomes A_t3 (3), A_t4 (4), A_t10 (1), A_t11 (2), and D_t3 (2) (Figure 2A). Most importantly, five of the twelve SNPs significantly associated with the trait were detected with the lowest P-value (−log₁₀(P) > 6.21); these five loci were positioned on chromosomes A_t3 and A_t4 of upland cotton and explained 9.42–12.88% of the total phenotypic variance (Figure 2A, Table 1). Interestingly, two SNP loci (rsA_t4:15572813 and rsA_t4:15573052) on chromosome A_t4 showed significant marker-trait associations with the highest −log₁₀(P) values, which were detected in BLUP and each of the four environments simultaneously (Figure 2A, Supplementary Figure S1, and Table 1). Additionally, three SNP loci (rsA_t3:43538238, rsA_t3:43631774, and rsA_t3:43631819) on chromosome A_t3 were associated with the target trait in BLUP and three of the four environments at the same time (Figure 2A, Supplementary Figure S1, and Table 1). These results indicated that there were two major-effect QTLs for lint percentage on chromosomes A_t3 and A_t4.

To identify favorable QTL alleles for lint percentage, the five significantly associated SNP loci (−log₁₀(P) > 6.21) in two major-effect QTLs were selected, which exhibited the minimum P-value and could explain the maximum phenotypic variation. Three SNP loci from chromosome A_t3 (rsA_t3:43538238, rsA_t3:43631774, and rsA_t3:43631819) were associated with lint percentage via GWAS; these SNP alleles were G/T, C/T, and A/G, respectively. Three types of haplotypes (GG-TT-GG, TT-CC-AA, and GT-CT-AG) were found because of a close linkage relationship among the three SNP loci. The GG-TT-GG haplotype included 308 lines and was considered the favorable haplotype (FH) because the average lint percentage (41.32%) of the haplotype was significantly higher than the mean value (35.81%) of the other corresponding haplotype (TT-CC-AA, the unfavorable haplotype, UFH) and included 26 lines. The remaining 21 lines were heterozygous haplotypes (HH, GT-CT-AG), with a mean lint percentage of 38.77% (Figure 3A). Furthermore, the FH accounted for a very large proportion of the accessions with higher lint percentage than that with lower lint percentage. For example, there was no FH in the accessions with a low lint percentage (< 30.50%), and there was no UFH in the lines with a high lint percentage (>45.50%), (Figure 3D). The other two SNP loci (rsA_t4:15572813 and rsA_t4:15573052) associated with the lint percentage were closely linked and were located on chromosome A_t4, separated by 238 bp. Similarly, the lint percentage of accessions with the FH (AA-AA, 41.36%), was higher on average than that of accessions with the UFH (GG-GG, 34.79%; Figure 3B). The FH accounted for a large proportion of the accessions with a higher lint percentage, and the UFH was not found among the upland cotton lines with lint percentage values >39.50% (Figure 3E). These results further confirmed that there were two major QTLs controlling lint percentage in the adjacent regions of the five associated SNP loci.

To further understand the cumulative effect of FHs on lint percentage, the 355 upland cotton lines were grouped into three groups according to the number of FHs. The cotton lint percentage among the three groups obviously increased with an increasing number of FHs (Figure 3C). For example, accessions with two FHs showed a 41.53% mean lint percentage, which was higher than accessions with one FH (37.33%) and those without a FH (33.09%). These results suggested that the genetic control of lint percentage exhibits a largely additive effect in upland cotton.

Two upland cotton (TM-1) reference genome sequences are currently available in a public database (http://www.cottongen.org). To identify the corresponding chromosomal locations of the SNP loci associated with the target trait between two upland cotton reference genomes, 250 bp of the sequence on either side of each associated SNP were extracted from one upland cotton reference genome (Li et al., 2015), and 501 bp of the sequences was aligned with the other upland cotton reference genome (Zhang et al., 2015). There was a clear correspondence between the chromosomal positions of the SNP loci associated with lint percentage between these two upland cotton reference genomes. For example, positions 15572813 and 15573052 on chromosome A_t4 corresponded to 52098956 and 52098717 on chromosome A08, respectively, while positions 43538238, 43631774 and 43631819 on chromosome A_t3 corresponded to 74801531, 74713290 and 74713245 on chromosome A02, respectively (Table 2).

To look for putative candidate genes in the neighboring regions of the SNP loci associated with lint percentage, we further determined LD blocks harboring five significant SNPs (−log₁₀(P) > 6.21). Two LD blocks were found on chromosomes A02 (73.25–74.75 Mb) and A08 (51.95–52.12 Mb) (Figures 4A,B). Based on the results of the genes annotated in the CottonGen database (http://www.cottongen.org), a total of 45 and 2 putative candidate genes were contained in two LD blocks on chromosomes A02 and A08, respectively. Furthermore, the expression levels of 47 putative candidate genes were analyzed using RNA-seq data from ten upland cotton (TM-1) tissues (Zhang et al., 2015). Six genes (Gh_A02G1229, Gh_A02G1262, Gh_A02G1259, Gh_A02G1263, Gh_A02G1268, and Gh_A02G1228) presented higher expression levels in the fiber-5 DPA and fiber-10 DPA than in the other eight tissues. Specifically, only one gene (Gh_A02G1268) was predominantly expressed in fiber-5 DPA, and its expression levels in fiber-5 DPA were 300-fold higher than in the seeds in TM-1 (Figure 4C), suggesting its potential role in increasing the lint percentage and improving fiber quality in upland cotton. Gh_A02G1268, a member of the myo-inositol-1-phosphate synthase (MIPS) gene family, was located 95.81 kb away from the peak SNP (A02:74713245). Its closest Arabidopsis homolog (At2g22240) is known to regulate embryogenesis and seed development (Mitsuhashi et al., 2008; Luo et al., 2011).

To further confirm these SNP loci associated with the lint percentage, our new GWAS was compared with previous linkage and association studies. Regarding previous QTL mapping studies, 281 SSR markers (Supplementary Tables S4, S5) containing QTLs for lint percentage were selected from 23 reports of QTL mapping (Supplementary Table S6), and 183 primer sequences corresponding to these markers (Supplementary Table S4) were obtained from the CottonGen Database (http://www.cottongen.org). The physical locations of these SSR primer sequences were mapped to the reference genome sequence (Zhang et al., 2015) via electronic PCR (e-PCR), among which 6, 4, 7, 8, and 18 SSR primers were mapped to chromosomes A02, A06 A08, A10, and D03, respectively (Figure 5). Three SNP loci from chromosome A02 (A_t3) (rsA_t3:43538238, rsA_t3:43631774, and rsA_t3:43631819) were mapped to an adjacent region of HAU2014; two SNP loci (rsA_t4:15572813 and rsA_t4:15573052) from chromosome A08 (A_t4) were positioned between MUSS167 and MUCS531; rsA_t11:9054901 near CIR166 was located on chromosome A10 (A_t11); and two SNP loci (rsD_t3:18792729 and rsD_t3:18809554) from chromosome D03 (D_t3) were mapped in the vicinity of many QTL-SSR markers, such as BNL2443, NAU3479 and NAU5444.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.