The 355 upland cotton germplasm resources collected by laboratories worldwide represent a natural population. Previous studies focused on early maturity (Li L. et al., 2021), fiber quality (Su et al., 2016b), fiber yield (Su et al., 2016a; Feng et al., 2022), and plant architecture component traits based on abundant phenotypic variations in this population (Su et al., 2018). These upland cotton varieties are from different countries and represent accessions resulting from more than 100 years of global upland cotton breeding. Seeds of the GWAS population used for phenotyping SV-related traits were collected from three environments, including Huanggang in Hubei Province (30°57′ N, 114°92′ E) in 2021 (E1: Huanggang-2021) and Sanya in Hainan Province (18°36′ N, 109°17′ E) in two consecutive years (2021 and 2022) (E2: Sanya-2021 and E3: Sanya-2022). The field experiments in Sanya and Huanggang were conducted following a randomized complete block design with two and three replications, respectively.

The phenotyping of SV-related traits was carried out by the sandponic method based on previously described methods (Si et al., 2022). Cotton seeds collected from the field were ginned, and cotton fuzz was removed by concentrated sulfuric acid. Then, all seeds were sun-dried for 2 days to break dormancy uniformly. A total of 150 plump seeds with uniform size and full grain were selected, disinfected with 15% sodium hypochlorite for 10 min, and then washed clean with distilled water. Then, each line was evenly planted in a plastic sand box containing 800 g of dry quartz sand with a size of 13 cm × 19 cm × 12 cm. Subsequently, the seeds were covered with 250 g of dry quartz sand, and 200 mL of distilled water was added. The number of germinated seeds was counted each day until the seventh day. All experiments were conducted in a phytotron with 16 h of light (25°C) and 8 h of darkness (18°C). Three biological replicates were included for each accession, and 50 seeds were used for each replicate. Moreover, three SV-related traits (GR, GP, and GI) were selected for measurement. The full name, abbreviation, and measurement method of each trait are listed in Table 1 as described by Yuan et al. (2019). The statistical analysis of the maximum value, minimum value, average value, etc., was performed using R software (version: 4.2.2).

The resequencing data (PRJNA389777) of the 355 upland cotton germplasms used in this study were reported in a previous study (Li L. et al., 2021). The Illumina HiSeq4000 platform was used for paired-end read sequencing, with an average sequencing depth of more than 10×. Based on previously released data, the new variation map of the natural population was employed in the ‘HaplotypeCaller’ module of GATK (version: 4.2.6.1) (Mckenna et al., 2010). Briefly, the variation detection process was as follows: (1) The quality of paired-end reads from 355 accessions was assessed using FastQC (version: 0.11.9) (Andrews, 2010); (2) Sequencing quality control was carried out with fastp software (version: 0.23.2) to obtain high-quality reads with the following parameters: ‘-w 16 -c -l 80 -5 -3 -W 4 -M 20 -f 10 -F 13 -t 3 -T 3 -q 20 -u 40’ (Chen et al., 2018); (3) All high-quality reads were mapped to the ‘TM-1’ (version: HAU_v1.1) reference genome using BWA (version: 0.7.17-r1188) (Li, 2013; Wang et al., 2019); (4) Then, Picard software (https://github.com/broadinstitute/picard) was used to sort the BAM file and mark duplicate reads; (5) The ‘HaplotypeCaller’ module of GATK (version: 4.2.6.1) was used to identify variant sites and perform SNP filtering with the following conditions: ‘QUAL <30, DP <1,340, DP >10,050, QD <2.0, MQ <35, FS >70, SOR >3, MQRankSum <-12.5, and ReadPosRankSum <-4.0’; (6) The SNP clusters with at least three SNPs detected within a 10-base window were removed; (7) SNPs within five base pairs of an InDel were filtered out by BCFtools software (version: 0.1.19-44428cd) (Danecek et al., 2021); and (8) SNPs with a minor allele frequency (MAF) <5% and missing rate <20% were discarded by VCFtools (version: 0.1.16) (Danecek et al., 2011).

Genome-wide association analysis was performed by combining 2,262,367 high-quality SNPs with the phenotype data of 355 upland cotton accessions collected in three environments for SV-related traits using linear mixed models in GEMMA (version: 0.98.3) and executed by vcf2gwas software (version: 0.8.7) (Zhou and Stephens, 2012; Vogt et al., 2022). P <1 × 10^-6 was used as the threshold to detect significant SNP loci. Additionally, the PVE by each marker was calculated as previously reported (Feng et al., 2022). The nucleotide diversity (π)a value was calculated using VCFtools based on the release years (before the 1950s, 1950s–1970s, 1980s–1990s, and 2000s–2020s) and geographical distribution (early maturity region: NSER, Yellow River region: YRR, Yangtze River region: YZRR, and Northwest Inland region: NIR) of the 355 accessions. The packages ‘CMplot’ (https://github.com/YinLiLin/CMplot), ‘LDheatmap’ (Shin et al., 2006), and ‘ggplot2’ (Wickham, 2011) in R software were used to generate Manhattan plots and for linkage disequilibrium (LD) block analysis and haplotype analysis.

Based on the ‘TM-1’ reference genome (HAU_v1.1) (Wang et al., 2019), the genes in the interval located 200 kb upstream and downstream of the significant SNPs were defined as candidate genes. The protein sequences of the candidate genes were obtained from Cottongene (https://www.cottongen.org/). Then, local BLAST software was used to compare the protein sequence of the candidate gene with the Arabidopsis protein database (https://www.arabidopsis.org) to obtain the homologous sequence, and the criterion was set to less than E^-60 (Johnson et al., 2008). The expression patterns of SV candidate genes in upland cotton were determined by RNA-seq and quantitative reverse-transcription PCR (qRT-PCR) analysis. RNA isolation method was performed as described by Feng et al. (2022). GhUBQ7 was used as an internal control. Quantitative analysis method was performed using a Roche real-time qPCR system (Light Cycler 480 II) and SYBR with three biological repeats. The public RNA-seq data (PRJNA248163) including SRR1695160, SRR1695161, and SRR1695162 were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/bioproject/). The Illumina Hiseq2000 platform was used to perform RNA sequencing on ‘TM1’ seeds soaked in water for 0, 5, and 10 h, and the paired-end clean reads length was more than 100 bp. The gene expression values were normalized by the average expression levels (log2) based on transcripts per million values. The clustered heat map was drawn by the R package ‘pheatmap’ (Kolde, 2012).

Resequencing of the natural population libraries by the Illumina HiSeq 4000 platform with 150 bp paired-end reads, as described in previous reports (Li L. et al., 2021), yielded approximately 65,013 million reads in total for the 355 cotton genotypes. Approximately 88.3% of the total bases were successfully mapped to the cotton reference genome, and the statistical sequencing depth corresponded to 11.7-fold in the 355 upland cotton accessions. A total of 2,262,367 SNPs distributed across the cotton genome with a MAF >0.05, and missing rate of resequencing data of less than 20% was used for the GWAS of the 355 cotton germplasm accessions, of which the At and Dt subgenomes contained 1,404,637 and 857,730 SNPs, resulting in an average SNP density of 993.44 and 1045.91 SNP/Mb, respectively ( Table 2 ; Figure 1 ). The percentage of the SNPs in each chromosome varied from 1.4% on chromosome D04 to 11.4% on chromosome A08 ( Figure 1 ). Most of the SNPs were located in intergenic regions (84.38%), whereas the exonic and intronic genome regions contained only 0.89% and 3.03% of SNPs, respectively ( Supplementary Table S1 ). In addition, SNPs in the coding regions (coding sequences, CDSs) included 33.26% synonymous mutations and 64.13% nonsynonymous mutations.

The three SV-related traits (GI, GP, and GR) of natural population accessions were measured in three environments. The values followed a normal distribution for the GI and GP but showed a skewed distribution for GR based on Shapiro–Wilk tests ( Supplementary Table S2 ). The frequency histograms of SV-related traits are shown in Figures 2A–I . The lowest average GI was 55.23 in the E1 environment, and the highest average GI was 58.92 in the E2 environment, with a coefficient of variation (CV) ranging from 6.52% to 12.02% ( Supplementary Table S2 ). For GP, the E1 environment had the lowest average value of 70.92%, while the E2 environment had the highest average value of 79.61%; the CV in the E1 environment (11.67%) was higher than that in the E2 environment (9.21%) and the E3 environment (9.37%) ( Supplementary Table S2 ). For GR, the lowest average value was 87.37% in the E1 environment, and the highest average value was 93.25% in the E2 environment, with a CV ranging from 3.71% to 10.48% ( Supplementary Table S2 ). Two-way analysis of variance (ANOVA) showed that genotype (G) and the genotype-by-environment interaction (G × E) had significant effects on the GI, GP, and GR (P < 0.001) ( Supplementary Table S3 ). Furthermore, the heritability of these three SV-related traits ranged from 74.23% (GR) to 81.75% (GP), whereas that of GI was 76.03% ( Supplementary Table S3 ). These results suggested that SV-related traits have extensive PV in the GWAS panel, which is suitable for further GWAS.

A total of 292 significant SNPs for three SV-related traits were identified on 11 chromosomes using the linear mixed model ( Figure 3 ; Supplementary Table S4 ; Supplementary Figures S1 - S3 ). Only 11 SNPs were identified in the At subgenome, whereas 281 SNPs were localized to the Dt subgenome. Among them, chromosome D03 had the highest number of SNPs (281), with a total of 254, and the range of -log₁₀(p) values was from 6.00 to 8.27. Furthermore, 26 stable SNPs were identified in a minimum of two environments (including for the best linear unbiased predictor, BLUP) or two traits, which were declared as six stable QTLs, focusing on chromosomes A09, A10, and D03. Notably, a QTL region (qGR-A09-1) located on chromosome A09 showed a strong SNP cluster associated with GR, which had a PVE of 6.76–8.56% and -log₁₀(P) ranging from 6.19 to 7.74. qGP-A10-1 on chromosome A10 had only one SNP that explained 8.15% of the observed PVE, with a LOD score of 7.39. Four QTLs on chromosome D03 (qGR/GI-D03-1, qGI/GR-D03-2, qGI/GP/GR-D03-3, and qGI/GP/GR-D03-4) were identified in two, three, three, and four environments, explaining 6.61–7.39%, 6.72–7.79%, 6.65–8.43%, and 6.61–8.90% of the observed PVE, respectively. Interestingly, a stable QTL (qGI/GP/GR-D03-3) region was revealed on chromosome D03 from 31.68 to 32.61 Mb and was flanked by regions associated with the GI, GP, and GR in the E1, E3, and BLUP environments. Thus, the QTLs qGR-A09-1 and qGI/GP/GR-D03-3 could be treated as major QTLs for further dissection.

In this study, a novel QTL, qGR-A09-1, exhibited a significant SNP cluster (rsA09_7745467, rsA09_7791621, rsA09_7878527, rsA09_7908017, rsA09_7954329, rsA09_7954353, and rsA09_7962794) occupying a physical region of 0.2 Mb on chromosome A09 ( Figure 4A ). Meanwhile, 22 genes were annotated in this QTL region based on the G. hirsutum reference genome (Wang et al., 2019), except for Ghir_A09G002720 and Ghir_A09G002760, which did not have annotation information ( Supplementary Table S5 ). We further conducted LD analysis on the significant SNP rsA09_7962794, and LD blocks were found in this region ( Figure 4A ). In this QTL interval, rsA09_7962794 on chromosome A09 showed a strong association with GR, with 7.95% of the PVE downstream of Ghir_A09G002730 ( Table 3 ). rsA09_7962794 had two haplotypes, GG and AA, which resulted in the accessions carrying the AA genotype having a significantly higher GR than those carrying the GG haplotype in three environments (P < 0.01) ( Figure 4B ). In addition, to gain a further understanding of the genetic characteristics of rsA09_7962794 in relation to geographic distribution, the 355 upland cotton accessions were divided into four groups: NIR, YZRR, YRR, and NSER. Interestingly, YRR and NSER showed an extraordinarily low frequency of the nonfavorable haplotype (GG), while the accessions obtained from YZRR and NIR had a relatively high frequency of the favorable haplotype (AA) (>75%) ( Figure 4C ). Furthermore, the genetic diversity of Ghir_A09G002730 decreased following the breeding period. Cotton accessions released before the 1980s showed greater diversity than accessions bred from the 1980s to the 2000s, while accessions bred after the 2000s showed the lowest diversity ( Figure 4D ). Specifically, Ghir_A09G002730 belongs to the pentatricopeptide repeat (PPR) superfamily protein family and has higher expression levels during the seed germination stage from 0 to 10 h than other genes ( Figure 4E ). The qRT-PCR analysis also showed that Ghir_A09G002730 had higher expression levels in the accessions (‘Liaomian27’ and ‘Xinluzhong35’) carrying the AA allele than in accessions (‘PB12-1-8’ and ‘Xiazao2’) with GG allele during the seed germination stage ( Supplementary Figure S4 ).

As mentioned above, another distinct SNP enrichment QTL region, qGI/GP/GR-D03-3, was detected for the GI, GP, and GR across multiple environments, which could explain the relatively high PVE of 6.65–8.43%, indicating that a major gene in this genomic interval may improve seed germination ( Table 3 ). Interestingly, 12 associated SNPs were located within the most significant haplotype block, which was almost 920 kb long and contained five haplotypes ( Figures 5A, B ). A haplotype analysis revealed that qGI/GP/GR-D03-3 had two major haplotypes according to SNP location. Comparatively, Hap1 had a higher GP than Hap1 ( Figures 5C, D ). In total, 46 candidate genes contained in the qGI/GP/GR-D03-3 region on chromosome D03 were identified. Among them, Ghir_D03G009280 was annotated as auxin response factor 9 (ARF9) in Arabidopsis ( Supplementary Table S6 ), and its homologs played a crucial role in seed dormancy. The RNA-seq and qRT-PCR assays also showed that Ghir_D03G009280 had higher expression levels during the seed germination stage, suggesting a positive regulatory effect ( Figure 5E ; Supplementary Figure S5 ).

SV is an indispensable indicator of seed quality, which directly affects the rapid and uniform germination of seeds and the robust growth of seedlings and affects the tolerance of plants to abiotic stress in the early stage of seedling growth (Qun et al., 2007; Fujino et al., 2008). In recent years, mechanical DS of cotton has been widely used due to its cost-saving and labor-saving advantages, leading to rapid and uniform seed germination becoming necessary conditions for high yield and mechanization in the cotton industry. However, seeds with low SV make it difficult for mechanical DS to achieve full seeding, which leads to problems such as subsequent filling of the gaps with seedlings and final singling of seedlings (Xie et al., 2014)—for example, Xinjiang Province is the major cotton-growing area in China and experiences serious saline–alkali stress (He et al., 2023). A high SV of cotton varieties will improve seed germination in the field and thus increase the yield. In addition, cotton breeding without plastic film in Xinjiang Province to eliminate “white pollution” has become popular. The germination rate and seedling emergence rate of seeds have higher requirements for cotton without plastic film (CWPF). CWPF needs to quickly establish robust seedlings after seed germination to resist the invasion of diseases, insect pests, adverse environments, and other factors in the field. Importantly, SV is the result of genetic and environmental factors and is thus often difficult to effectively select in conventional breeding (Dai et al., 2022). This study utilized high-throughput sequencing to generate widely distributed SNP markers that cover the whole genome ( Figure 1 ), and over 200,000,000 high-quality SNPs were detected in a diverse set of 355 cotton accessions. Combining phenotype data from multiple environments for GWAS analysis can be used to effectively identify genetic loci and candidate genes that improve SV in upland cotton, providing an effective way to improve cotton yield in Xinjiang when using the MAS method.

To ensure the accuracy of the GWAS results, phenotypic identification in multiple environments was conducted with at least three replicates per environment. The three SV-related traits (GI, GP, and GR) were measured for seeds collected from three locations: E1, E2, and E3. Among them, GR and GP did not show an absolute normal distribution, which was also found in previous studies (Dai et al., 2022; Si et al., 2022), indicating a complex genetic basis for these SV-related traits. Through phenotypic correlation analysis, it was found that there were significant positive correlations between the three traits. The GI showed a strong correlation with GR and GP (0.71 and 0.76, respectively) ( Figure 2J ). The highest GI was accompanied by the highest GP and GR, which is consistent with previous findings (Si et al., 2022). Furthermore, according to the measurement results for each trait, the CV of SV-related traits in upland cotton is affected by the environment ( Supplementary Table S2 ), resulting in different variations in the seeds of each accession harvested in different planting locations and years—for example, the CV of the GI and GR in E1 showed a larger range of variation than that in E2 and E3. Previous studies have shown that the environment in the planting area has a great influence on the growth and development of seeds (Fenner, 1992). It is speculated that the E2 and E3 (Sanya City, Hainan Province) environments with tropical climates are more suitable environments for seed growth, and the performance of the seeds may be relatively stable. In contrast, the E1 environment (Huanggang City, Hubei Province) has high precipitation and temperature during the seed maturation period, which can affect the success of pollination.

In the past two decades, GWAS has become a powerful and widely used tool for analyzing the genetic mechanisms underlying complex quantitative traits in crops (Tibbs Cortes et al., 2021). At present, most research on SV mainly focuses on the mechanism under stress in upland cotton (Sun et al., 2018; Yuan et al., 2019; Zheng et al., 2021), while genetic analysis of SV-related traits associated with normal seed germination is less common (Si et al., 2022). In this study, a GWAS panel was used to measure three SV-related traits of seeds harvested in three environments. In total, six significant QTLs were stably identified on three different cotton chromosomes ( Table 3 ), including 26 SNPs. Numerous studies have reported that several pathways are involved in regulating SV in plants, such as phytohormone signaling (GA, ABA, and auxin), amino acid metabolism, and the reactive oxygen pathway, which play a crucial role in the seed germination process and have a significant effect on the molecular mechanisms related to SV (Reed et al., 2022). It has been reported that high concentrations of ABA promote dormancy and inhibit seed germination, while high concentrations of GA promote seed germination by reversing dormancy, leading to an endogenous balance of the ABA/GA ratio but not the absolute hormone contents (Finch-Savage and Leubner-Metzger, 2006; Chen H. et al., 2020). Ghir_A09G002650 was annotated on chromosome A09, belonging to the GA-regulated family of proteins and encoding a protein containing the GASA domain, which is most closely related to the known homolog GASA14 in Arabidopsis. GASA14 regulates the increase in plant growth through GA induction and DELLA-dependent signal transduction, which could increase resistance to abiotic stress by reducing the accumulation of ROS (Sun et al., 2013). Thus, it is speculated that Ghir_A09G002650 has the potential to improve the SV of cotton under stress. MYB-type and bHLH-type transcription factors have been reported to be involved in the regulation of seed germination signaling in plants (Penfield et al., 2005; Reyes and Chua, 2007; Kim et al., 2015; Wang X. et al., 2022; Xu et al., 2022). Specifically, Ghir_D03G006550 is in the qGI/GR-D03-2 region and is homologous to MYB52. It has been previously shown that its shared common targets with ERF4 regulate the development of the seed coat in Arabidopsis (Ding et al., 2021). Ghir_D03G010510 encoded bHLH-type family proteins in the QTL region of qGI/GP/GR-D03-4, sharing 35.52% sequence identity with the PIF8 protein in Arabidopsis, which binds to promoter regions of AtPIF6. The expression level of AtPIF6 during seed development plays a crucial role in establishing primary seed dormancy levels (Peters et al., 2010).

Notably, Ghir_A09G002730 and Ghir_D03G009280 were detected in two distinct enriched regions located on chromosome A09 (qGR-A09-1) and chromosome D03 (qGI/GP/GR-D03-3) ( Figure 3 ). Interestingly, Ghir_A09G002730, within the strong-LD region at 21.9 kb upstream of rsA09_7962794 and highly expressed during the development of seed germination ( Figures 4A, E ), encodes a PPR superfamily protein in Arabidopsis. SOAR1 belongs to the PPR protein family and acts as a core negative regulator downstream of ABAR and upstream of ABI5, participating in ABA signaling regulation of seed germination and seedling growth processes (Ma et al., 2020). We also discovered that cotton accessions carrying rsA09_7962794-A with a higher GR had a much higher allele frequency for Ghir_A09G002730 in YZRR and NIR than in YRR and NSER ( Figures 4B, C ). It is possible that the planting mode of seedling raising and transplanting in YZRR and mechanized planting in the NIR all employed single-seed sowing, which increased the selection frequency of rsA09_7962794-A. In addition, we compared the genetic diversity of the region on chromosome A09 containing Ghir_A09G002730 in different breeding periods, and it was found that cultivars bred after the 2000s had lower genetic diversity than cultivars from other stages, implying that with the continuous increase in cotton SV during the breeding process, this gene was associated with artificial selection ( Figure 4D ). Therefore, it is reasonable to postulate that Ghir_A09G002730 is a new candidate gene influencing SV in cotton. Ghir_D03G009280 caught our attention based on the gene annotation of cotton. This gene encodes an auxin response factor. Recent studies have shown that ARF16 interacts with ABI5 and positively regulates the ABA response during seed germination (Mei et al., 2023). Furthermore, Ghir D03G009280, tightly linked with haplotype Hap1, showed a significant association with GP ( Figure 5C ), and materials carrying the Hap1 haplotype had longer roots ( Figure 5D ). The RNA-seq analysis showed a high expression level of this gene during seed germination ( Figure 5E ). From the above-mentioned results, we inferred that Ghir_A09G002730 and Ghir_D03G009280 were two major candidate genes that may play an important role in cotton SV.