Cotton is one of the most important cash crops used for natural textile fibers and oilseeds. It is cultivated worldwide as the primary raw material for the textile industry (Chen et al., 2007). It is derived from the Arabic word “quotn,” “gutum,” or “kutum” meaning those crops which possess spinnable fibers on their seed coat (Lee, 1984). Among the plant species commonly used by people, fiber-yielding crops hold the second position after food crops. Cotton serves as an excellent model for studying various aspects, such as plant cell elongation, evolution, cellulose biosynthesis, and plant polyploidization (Qin and Zhu, 2011; Huang et al., 2021). It belongs to the genus Gossypium in the family Malvaceae and is one of the largest genera in the Gossypieae tribe owing to its high species diversity. It encompasses approximately 44 diploid (2n = 2x = 26) and seven tetraploid (2n = 4x = 52) species (Grover et al., 2015; Ditta et al., 2018).

Cotton is cultivated in more than 123 countries, encompassing regions ranging from arid to semi-arid areas in the tropics and subtropics. Among the four cultivated species, upland cotton (Gossypium hirsutum) constitutes ~95% of the global cotton production because of its high adaptability to various environments and significant yield (Chen et al., 2007). Thus, compared with other cultivated species, most breeding and improvement programs related to fiber yield, fiber quality, insect resistance, and drought tolerance, are mainly focused on upland cotton. Genetic improvement for high cotton fiber yield has always remained the primary focus of cotton breeders to increase their productivity. With progress and advancement in the textile industry, the demand for high fiber yield and quality is growing exponentially. Although fiber quality is an essential trait in cotton breeding programs (Rong et al., 2007; Said et al., 2013), enhancing cotton fiber yield using Indian cotton varieties is still a primary goal. In comparison to the rest of the world, India has the largest land area under cotton cultivation (12,150 thousand hectares), yet its productivity is greatly compromised (only 457 kg/ha) by 2022 (US Department of Agriculture, https://usda.library.cornell.edu/). In India, most cotton cultivars are released using conventional breeding techniques based on morphological traits that are affected by agronomic practices and environmental changes. Consequently, the diversity and quality of elite cotton are reduced daily because of the narrow genetic background, domestication, and selection of elite cotton cultivars (Iqbal et al., 2001; Rungis et al., 2005; Abdurakhmonov et al., 2012); thus, genetic improvement in cotton fiber yield and other quality traits is still challenging.

Cotton fiber yield is a multifaceted quantitative trait influenced by various traits, including boll weight (BW), boll number (BN), seed cotton weight (SCW), lint percentage (LP), lint index (LI), seed index (SI), first fruit branch position (FFBP), plant height (PH), flowering period (FP), fruit spur branch number (FSBN), and many others (Li F. et al., 2018; Sun et al., 2018). Such quantitative traits are governed by both quantitative trait loci (QTLs) and environmental factors, which are difficult to concurrently improve using traditional breeding techniques. However, the advancement of applied genomics research has introduced the use of QTL-linked or QTL-associated molecular markers in marker-assisted selection (MAS) and/or genomic selection programs. These emerging techniques offer promising avenues for enhancing the efficiency of cotton breeding and simultaneously targeting specific traits. QTL mapping has been extensively employed to analyze the genetic variations underlying complex traits in cotton, including fiber quality and yield component traits (Rong et al., 2007; Said et al., 2013; Said et al., 2015). Over the last 20 years, biparental linkage mapping in upland cotton has led to the identification of numerous QTLs associated with cotton yield (Zhang et al., 2005; Abdurakhmonov et al., 2007; Shen et al., 2007; Wan et al., 2007; Liu et al., 2012; Yu et al., 2013a; Yu et al., 2013b; Wang et al., 2014; Liu D. et al., 2015). Nevertheless, biparental QTL mapping encounters two main challenges: limited allelic diversity and restricted genomic resolution due to the relatively low number of recombination events that occur throughout the development of the mapping population (Jannink et al., 2001; Flint-Garcia et al., 2005; Hall et al., 2010).

To overcome these limitations of bi-parental QTL mapping, association mapping offers an alternative approach for mapping QTL, which relies on linkage disequilibrium (LD). It determines whether particular alleles in a population are more frequently associated with specific phenotypes than expected, thereby providing insights into the genetic basis of complex traits (Flint-Garcia et al., 2003). In addition, this analysis can be applied to a large natural population, allowing the identification of traces of genetic crossovers and loci responsible for traits at a much higher resolution than previously reported techniques. In recent years, genome-wide association studies (GWASs) have emerged as a more precise and cost-effective method for identifying significant QTLs or genes linked to complex traits than linkage mapping. The first association mapping study in cotton focusing on fiber quality was reported in G. arboretum (Kantartzi and Stewart, 2008). Subsequently, association mapping has been extensively applied in upland cotton to investigate crucial economic traits, including fiber quality (Abdurakhmonov et al., 2008; Nie et al., 2016; Liu et al., 2020; Song et al., 2021), seed oil content (Liu G. et al., 2015; Yuan et al., 2018; Zhao et al., 2019), fiber yield (Ademe et al., 2017; Sun et al., 2018; GUO et al., 2021; Niu et al., 2023), and biotic and abiotic stresses (Fang et al., 2010; Bardak et al., 2021; Xu et al., 2021; Zhang et al., 2021; Zhao et al., 2021), as well as those associated with epistasis and environmental interactions (Jia et al., 2014).

To date, studies on association mapping using Indian germplasm are scarce and were first reported by Handi et al. (2017), who used a cotton 63 K SNP chip to associate fiber yield and fiber quality traits in 201 upland cotton germplasm lines. Owing to their limitations, SNP chip-based studies only fetch SNPs at a specific location in the genome (where relevant information is gathered), and genotyping-by-sequencing (GBS) methods, on the other hand, proved to be an attractive approach for discovering and genotyping high-density SNPs. Although other marker classes can also be used for this study (for example, Kumar et al. (2021) used SSRs), SNPs demonstrated higher resolution in revealing genetic relatedness as well as delineating population structure in crops. Additionally, since SNP markers are abundant and primarily derived from genes, genetic diversity studies using these markers can reveal the functional variation that can be used in association mapping studies for specific traits (Singh et al., 2013).

Therefore, the specific goals of the present study were to (a) examine the phenotypic variability within the upland cotton germplasm released from different cotton growing belts in India; (b) explore the QTNs (qualitative trait nucleotides) or SNP markers underlying fiber yield traits; (c) compare identified QTLs with previously reported QTLs for yield traits, if available; and (d) identify and validate through expressional profiling of putative candidate genes found within the genomic regions controlling yield traits. Our results using 117 diverse Indian germplasms and 2,41,086 SNPs identified many stable and novel QTL/genes that may offer crucial information on the genetic control of fiber yield traits in cotton. This information will aid in improving cotton yield and the development of elite Indian cotton varieties through marker-assisted breeding programs.

Cotton is one of the most important natural fibers and is a raw material for many textile industries. Its yield can significantly influence the textile industry because it is directly related to overall supply and demand dynamics, production capacity, pricing, innovation, global trade, and sustainability efforts within the industry. Cotton yield is a complex quantitative trait governed by multiple yield-related descriptors that are difficult to improve simultaneously using traditional breeding methods (Wang et al., 2019). Thus, the identification and characterization of genetic factors for targeted traits and their manipulation through conventional breeding coupled with genomic tools have consistently been the primary goals for cotton breeders, enabling a more precise selection of genotypes in the pursuit of developing varieties with higher yield potential. In the present study, GWAS was performed based on the natural population of 117 Indian upland cotton germplasms and genotyping-by-sequencing methods to deploy markers for six yield-related traits. Association mapping study is an effective tool for mapping complex quantitative traits to identify of key genes associated with such traits in many plants (Huang et al., 2012; Horton et al., 2014; Crowell et al., 2016; Varshney et al., 2019b; Wang et al., 2020). It is an analytical method and its ability to resolve associated loci relies on factors such as the size of the experimental population, marker density, and the selection of appropriate statistical tools or GWAS models for marker-trait identification (Liu et al., 2016). Population size has a great impact on association mapping studies, and it has been reported that a relatively large population size ensures sufficient genetic variation, thereby influencing the detection of significant QTLs (Su et al., 2016; Huang et al., 2017). Population sizes ranging from 95 to 800 have been previously studied in upland cotton in association studies (Gapare et al., 2017; Sun Z. et al., 2017; Dong et al., 2019; Kumar et al., 2021). Although our population size of 117 released varieties from different agro-climatic zones was not sufficiently large, they are suitable for mapping analysis because they stabilize after many years of genetic recombination (Nordborg et al., 2002). Along with population size, phenotypic variances in a single environment affect the accuracy of phenotypic data, thereby affecting the reliability of the mapping study. To rectify this error, the use of multi-environment and unbiased prediction data in association analysis has been proven effective in previous studies (Roorkiwal et al., 2018; Tomar et al., 2021). Phenotyping from multiple field locations effectively eliminates the influence of the environment and aids in the interpretation of environment-specific as well as general QTL (Gutiérrez et al., 2015). In the present study, we used two sites for the evaluation of six yield-related phenotypic traits: Punjab Agricultural University, Punjab (E1), and Tamil Nadu Agricultural University (E2), which are the two main cotton-growing belts in India. In addition, there is a large difference in the geographical position as well as the climatic conditions between these two zones (https://www.cicr.org.in/pdf/long_staple), and all six traits exhibited abundant phenotypic variation ranging from 6.32% to 24.7% and 6.17% to 20.1% in E1 and E2, respectively. Our findings also showed the presence of medium to high heritability (>50%, except for LP in E2), illustrating the importance of these germplasm lines for genetic improvement and showing the same phenotypic expression if sown in the same environment (Bhat et al., 2022). High and stable heritability has always been proven to be useful in determining the strong association between markers and traits of interest (Courtois et al., 2013). In addition, the correlation study among the traits showed a significant positive correlation between SCY and BN, and BW, and a maximum negative correlation between LP and SI ( Figure 2 ) which corroborates the findings of Dong et al. (2018) and Sun et al. (2018). Thus, the presence of significant variation among the genotypes in different environments along with the strong correlation led us to explore markers associated with six yield-related traits through GWAS analysis.

Association mapping using various SSR or SNP chips has been previously reported in upland cotton (Mei et al., 2013; Gapare et al., 2017; Handi et al., 2017; Huang et al., 2017; Kumar et al., 2021); however, the development of high-resolution sequencing technologies such as GBS, SLAF-seq, and RAD-seq offers significant advantages and has led to tremendous progress in the development of numerous SNP markers for genetic mapping in cotton (Islam et al., 2016; Su et al., 2016; Geng et al., 2020; Wang et al., 2021) and other crops (Arruda et al., 2016; Maldonado et al., 2019; Ravelombola et al., 2021). Genotyping by sequencing (GBS) provides a comprehensive view of the genome by generating genetic markers that spread across the entire genome. In contrast, SSR markers and SNP chips have limited coverage, and may not capture the full genetic variation present in a population. This increased marker density improves the resolution and efficiency of GWAS, enabling the detection of more genetic associations and fine mapping of genomic regions. In addition, SNP chips often rely on preselected SNPs that may not fully represent the genetic diversity of a particular population or species, conversely, GBS allows for unbiased genotyping, reduces ascertainment bias, and enables a more accurate assessment of genetic variation. GBS technologies have enabled the rapid sequencing and genotyping of breeding populations, allowing plant breeders to accomplish genomic diversity, GWAS, genomic selection (GS), and marker development in many species without prior knowledge of the species genomes (Poland and Rife, 2012). In our study, a substantial number of high-quality SNP markers (2,41,086) were identified using the GBS method with an average density of 1SNPs/11.36 kb. The number of markers reported in our study was relatively higher than that reported by Wang et al. (2021); however, it was almost consistent with the findings of Geng et al. (2020). The higher number of markers in the current study was sufficient to conduct GWAS analysis, offering potential advantages in uncovering additional prominent loci and candidate genes (Wang et al., 2018).

Upland cotton has an extensive and intricate history of domestication and breeding, with a narrow genetic background, and is mostly influenced by geographical isolation and gene flow (Huang et al., 2017). Furthermore, understanding the population structure and relatedness in the association mapping panel is crucial for elucidating the heterogeneity of the genetic architecture and controlling false associations (Lu et al., 2015). Therefore, it is crucial to consider the population structure and degree of relatedness among individuals in association-mapping studies. In our study, the association panel was categorized into two clades (CGP1 and CGP2) based on molecular analyses using STRUCTURE software but were not completely consistent according to the geographical origin, as reported in previous studies (Gapare et al., 2017; Sun Z. et al., 2017; Song et al., 2019). In addition, clustering analysis using Nei’s genetic distance, kinship matrix, and PCA calculation methods showed good consistency with population structure analysis. The absence of geographical correlation in our germplasm might be attributed to the relatively high level of gene exchange and interspersed introduction or crossbreeding among the germplasms across different geographical regions of India. Overall, these results highlighted that the germplasms were not highly structured and exhibited weak relatedness, but considering the greater continuous phenotypic variations among the six yield-related traits ( Table 1 ), the association population was further used in the GWAS analysis (Yano et al., 2016; Sun Z. et al., 2017; Song et al., 2019).

Several studies have indicated that effective control of the discovery of false positives resulting from population structure in crops may not be entirely achieved (Hamblin et al., 2011; Lipka et al., 2015). Therefore, to address these errors, we investigated various statistically robust models for the genome-wide association study using the GAPIT package in the ‘R’ programming. The strategy of utilizing two or more different models to identify significant markers in cotton has been reported in several studies (Badigannavar and Myers, 2015; Abdullaev et al., 2017; Sethi et al., 2017; Li C. et al., 2018; Su et al., 2018). In the present study, we implemented six different models; three univariate models (GLM, MLM, and CLMM) and three multivariate mixed models (MLMM, FarmCPU, and Blink). According to Segura et al. (2012), the use of a multivariate GWAS method addresses the issue of confounding between covariates and the test marker, and lowers the false discovery rate (FDR) compared to univariate GWAS when employing the same threshold, enabling the detection of a greater number of QTLs. By examining the Q–Q plots, we determined the most suitable models, and our findings demonstrated that the MLMM exhibited the best-fit model for our trait data, unlike the other models that exhibited early deviations from the expected line. The MLMM model has also shown successful results in other plants, such as Cannabis (Watts et al., 2021), tomato (Zhao et al., 2022), apricot (Omrani et al., 2019), and wheat (Mihalyov et al., 2017). Detected by the MLMM, 205 significant SNP from two environments and BLUP were found to be associated with six yield-related traits ( Supplementary Table S6 ). These markers were unevenly distributed on all chromosomes of cotton (except for Chr D13), suggesting that the regulation of these traits involves a complex gene network. This finding aligns with those of previous reports, highlighting the intricate genetic control of cotton yields (Rong et al., 2007; Sun et al., 2018). Using these significantly associated markers, QTL were determined within the ±200 kb upstream and downstream regions, as reported previously (Sun Z. et al., 2017; Song et al., 2019). In this study, 134 QTLs were derived from 205 significant markers ( Supplementary Table S7 ). Interestingly, some of these significant SNPs or QTLs were associated with multiple yield-related traits at pleiotropic loci. It has been reported that complexity and linkage between complex traits are common in many crops (Van Tienderen et al., 1996; Yan et al., 2011; Raman et al., 2019; Li A. et al., 2022). In upland cotton, four pleiotropic loci related to fiber yield, fiber quality, and flowering date have been reported by Wang et al. (2021). Similarly, four and nine pleiotropic loci have been reported for concurrent alterations in both lint yield and fiber quality traits in cotton (Geng et al., 2020; Li Y. et al., 2023). In our study, eight QTLs (six within the same genomic interval and two overlapping with adjacent QTLs) were found to be associated with six yield-related traits. QTL qGhBN_SCY_D6-1, qGhBN_SCY_D6-2, qGhBN_SCY_D6-3 was simultaneously associated with BN and SCY, having common significant markers NC_053442.1_36092628, NC_053442.1_37082185, NC_053442.1_37984052. Three other QTLs, qGhSI_LI_A5, qGhLI_SI_A13, and qGhLI_SI_D9, were concurrently associated with SI and LI, of which qGhLI_SI_A13 had a common marker (NC_053436.1_78553469), while the other two were overlapping QTLs. Similarly, qGhBW_SCY_A10 and qGhLP_BN_A8 were pleiotropic QTLs for BW/SCY and LP/BN, respectively. Additionally, correlation analysis of these yield-related traits showed a significant positive correlation (between SCY and BW, BN), and the maximum negative correlation between LP and SI indicated that these are the favorable loci that could be the genetic foundation for the correlation between these traits. Thus, the pleiotropic loci identified in this study provide information to explore the molecular mechanism that explains the simultaneous enhancement of six intricate yield-related traits in cotton, which will be helpful in the selection of germplasms in crop improvement programs for yield traits in Indian upland cotton.

In the past few decades, hundreds of QTLs or GWAS signals associated with cotton fiber yield have been identified in different intraspecific and interspecific populations through various linkage and association mapping methods (Said et al., 2015; Huang et al., 2017; Ma Z. et al., 2018). Of the 134 QTLs detected in the present study, 95 QTLs overlapped within the genomic interval or were adjacent to the genomic position of the QTLs and GWAS signals identified in a previous study ( Supplementary Table S7 ). Interestingly, some of the QTLs were trait-specific, as reported earlier within the same genomic interval. For instance, QTLs (qGhBN_D5-1, qGhBN_D6-1, qGhBN_D6-2, qGhBN_D6-3, and qGhBN_D6-4) for BN overlap within TM58714_TM58742_TM58749 (Zhu et al., 2021), and qBPP-1-43.5 (Zhang et al., 2011). Similarly, for BW, QTLs (qGhBW_A13-2, qGhBW_D4-2, qGhBW_D7-1, qGhBW_D9-3, qGhBW_D12-1) showed consistent results with TM47610_TM47614, TM56685_TM56686, TM63749 (Zhu et al., 2021), qBW-C16-1 (Wu et al., 2009) i15830Gh (Sun et al., 2018), and qBW-06A-c26-1 (Yu et al., 2013b) BW QTLs reported previously. LI (qGhLI_D6-2, qGhLI_D9-2) corresponds to QTLs qLI-D6-1 (Liu et al., 2012) and Br7_Lt%_23 (3,47+) (Lacape et al., 2013) LI QTLs reported in previous studies. For LP, three QTL (qGhLP_A6-2, qGhLP_D6-1, and qGhLP_D8-1) matched with the LP GWAS signals reported by Ma Z. et al. (2018) (A06_102555770), and QTL reported by Yu et al. (2013a) (F2:3-qLP-c25-1), and Chen et al. (2010) (qLP-D8-1/qLP-F2:3-JES-1a), respectively. Four QTLs (qGhSCY_A13-1, qGhSCY_D6-3, qGhSCY_D2-1, and qGhSCY_D2-2) were associated with SCY in the present and in the previously reported SCY QTLs: qSCY-07A-c13-1, qSCY-06A-c25-1 (Yu et al., 2013b), qSY-D2-1 (Wang et al., 2007), and TC-qSCY-c14-1 (Yu et al., 2013a), respectively. Interestingly, pleiotropic QTLs for BN and SCY (qGhBN_SCY_D6-1, qGhBN_SCY_D6-2, and qGhBN_SCY_D6-3) showed congruency for both BN (qBPP-1-43.5) and SCY (qSCY-06A-c25-1) traits, as reported previously (Zhang et al., 2011; Yu et al., 2013b). Pleiotropic QTL (qGhBW_SCY_A10 and qGhLI_SI_D9) were consistent with the BW association signal A10_99131954 (Ma Z. et al., 2018) and qSI-D9-1 (Shen et al., 2007) located on chromosomes A10 and D09, respectively. Therefore, our findings corroborate those of previous studies, validate the authenticity of current GWAS results, and increase confidence in the reliability of some QTLs/SNPs. These QTL/SNPs, exhibit stable inheritance and were consistently detected across diverse segregating populations with varying genetic backgrounds and through different mapping methods, and have a significant potential for future breeding programs aimed at improving cotton yield in India.

High-quality allelic loci are precious assets for agricultural breeding initiatives, and the identification of favorable alleles is an efficient approach to enhance the traits within crop plants (Su et al., 2016). Many recent studies have identified haplotype alleles for important traits, such as haplotype analysis of yield-related traits in soybean (Bhat et al., 2022), agronomically important traits in Arabidopsis (Lu et al., 2019), grain quality traits in rice (Wang et al., 2017), high thousand-kernel weight in wheat (Sun C. et al., 2017), grain yield, and flowering time under drought and heat stress conditions in maize (Yuan et al., 2019), which have shown considerable potential for the identification of traits and crop improvement. In cotton, Song et al. (2019) reported the haplotype analysis of two significantly associated SNPs with the lint percentage (LP) trait had a positive effect on LP, and these favorable alleles can be pyramided in a target line by marker-assisted selection. Similarly, two major haplotypes for fiber length and strength in cotton have been identified on chromosomes Dt11 and At07 (Sun Z. et al., 2017). In the present study, haplotype analysis of pleiotropic QTLs identified for different traits was conducted using a 10–20 kb window. The combination of favorable alleles identified within H003 (haplotype 3) was significantly higher for BN, BW, SCY, LI, and SI in the pleiotropic QTL qGhBN_SCY_D6-2, qGhSI_LI_A5, and qGhBW_SCY_A10, respectively. Eleven SNP combinations within haplogroup H001 were greater than other haplotype groups in the pleiotropic QTL qGhLI_SI_D9. Our results revealed that the combination of favorable alleles led to the identification of haplogroups that regulate a diverse range of phenotypes and significant phenotypic variation in yield-related traits in cotton. All identified diverse haplotypes can be further pyramided as a targeted line by marker-assisted breeding in cotton breeding programs. Thus, haplotype-based breeding strategies will aid in choosing favorable plant genotypes that carry advantageous haplotype alleles that have great potential for crop improvement (Varshney et al., 2019a).

Several genes associated with yield traits, such as Gh_D08G2376 (Huang et al., 2017), Gh_D12G2344 (Sun et al., 2018), Gh_D05G1124 (Song et al., 2019), Gh_A02G1268 (Su et al., 2016), Gh_D02G0025 (Ma Z. et al., 2018), AHP5 (Fang et al., 2017), Gh_A02G0111 (Niu et al., 2023), and Ghir_A08G009110 (Feng et al., 2022) have been previously reported in GWAS using different association mapping studies. In the present study, 2,509 candidate genes were identified within the confidence interval of the identified QTLs. Previous studies have suggested that genes preferentially expressed at different stages of fiber, ovule, and seed development may be involved in fiber yield and quality (Huang et al., 2017); therefore, we selected 40 highly active genes showing significant expression profiles ( Supplementary Table S8 ). For BN, three candidate genes were identified; the GhGELP isoform has been reported to have a direct function in ovule, fiber, and seed development in cotton (Ma R. et al., 2018). GhCBSX5 has reported to be involved in other biological processes (Ali et al., 2021); however, its role in fiber development in cotton has not yet been reported. GhCYP has no functional information in cotton; its isoform has been reported to play a significant role in drought tolerance in wheat (Zang et al., 2010). Gene annotation of the BW reported QTL intervals identified five active genes. GhCSGL3 plays a direct role in secondary cell wall biosynthesis leading to enhanced lint yield and quality in cotton (Li et al., 2015; Zhang et al., 2015), whereas GhAGXT2 and GhEP1-3 play a positive role in abiotic stress in cotton (Li L. et al., 2022). Two other genes, GhSMP1 and GhGRP4, were reported to be novel genes governing boll weight traits in the present study. Six candidate genes have been reported for LI traits of which two genes (GhKNAP2 and GhGA20OX1) were functionally validated to play a major role in the cotton fiber development process (Xiao et al., 2010; Gong et al., 2014). The other two genes (GhLEA and GhbZIP11) have been reported to be stress-responsive in cotton (Liang et al., 2016; Magwanga et al., 2018). In addition, the functional roles of GhCYP and GhMIF2 have not been reported in cotton, but its isoform has been reported to play a vital role in petal elongation in Gerbera hybrida (Han et al., 2017). Similarly, expressional analysis of candidate genes for LP traits resulted in the identification of seven active genes. Three genes (GhLEA-D-19, GhPDF1, and GhABCC15) have direct roles in cotton seed germination, fiber cell initiation, and elongation (Dure III and Galau, 1981; Zhu et al., 2003; Deng et al., 2012). The other four genes (GhGUX1, GhBBE18, GhMADS23, and GhChlADR1) were newly identified genes and might play a vital role during fiber development based on their expression profile; however, the functions of these in cotton remain to be further explored. Interestingly, for SCY highest number of genes was found (eleven) out of five GhPG, one GhPL, and one lncRNA reported encodes the highly significant marker (NC_053427.1_2009750) showing the highest phenotypic effects (~16.69%) could be the most promising genes found in the present study that could play a major role in fiber yield trait, which also corroborates with previously reported studies (Li Z. et al., 2023; Sun et al., 2020). Another gene, GhAP2/ERF/AIL5, has also been reported to play a key role in the growth and development of cotton plants (Zafar et al., 2022). GhPUP4, GhDCTPP1, and GhGGAT2 have no direct role in cotton; however, GhPUP4 has been reported to play a major role in the enhancement of grain size increase in rice (Xiao et al., 2019). SI has five active genes, (GhMYB22, GhAGP9, GhZAT10, GhACO3, and GhFLA7), which play an active role in fiber initiation elongation and freezing tolerance, as reported previously (Shi et al., 2006; Huang et al., 2013; Li P. et al., 2023). Additionally, pleiotropic QTLs (qGhBN_SCY_D6-3 and qGhSI_LI_A5) have three active genes (GhPPR, GhCHUP1, and GhSCPL42), whose functions in cotton fiber development have not been deciphered; however, their isoforms have been reported to play a functional role in defense mechanisms against abiotic and biotic stress (Wang et al., 2022). All 23 stable and 17 novel genes identified in the current study exhibited a high expression profile, which renders them promising candidate genes for future investigations and their functional validation would reveal their role in cotton yield improvement through functional genomics approaches. In conclusion, the present study unveiled a rich source of genetic elements, including SNPs, QTLs, and putative candidate genes associated with fiber yield traits in Indian upland cotton.

The present study has made significant strides in understanding the genetic structure and diversity of the Indian cotton germplasm, the identification of SNPs associated with fiber yield traits, and the subsequent identification of potential candidate genes. A weak population structure in the Indian cotton germplasm revealed two subgroups when the population structure was analyzed using a variety of approaches, including PCA, NJ tree, and kinship analysis. A low level of genetic relatedness among the genotypes was observed in the kinship matrix, which is required for breeding programs because it preserves the genetic diversity with wide adaptability. A total of 205 significant SNPs associated with six yield-related traits were identified through the GWAS and further delineated into 134 QTLs. Interestingly, several of these QTLs showed pleiotropic effects, indicating that these loci govern several traits and are thus advantageous candidate loci for introduction into breeding programs aimed at enhancing fiber yield. Additionally, 2,509 unique candidate genes were identified within the vicinity of these QTLs. The biological processes of the trait-associated genes were revealed using gene ontology enrichment analysis. While analyzing the public domain RNA-seq data, we identified 40 potential candidate genes across various cotton fiber developmental stages, several of which are known to be associated with fiber yield, while others need further functional validation to decipher their role in cotton yield. In conclusion, the present study unveiled a rich source of genetic elements, including SNPs, QTLs, and putative candidate genes associated with fiber yield traits in Indian upland cotton. These findings provide a solid foundation for further research on the functional roles of these genetic elements and their potential utilization in breeding programs to improve cotton fiber yield in India. To clearly determine the involvement of candidate genes governing cotton yield traits, further functional validation is necessary. It is anticipated that this effort will make a substantial contribution to the MAS breeding of high-fiber-yielding cotton varieties, thereby enhancing cotton farming productivity and the sustainability of the Indian cotton industry.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI SRA accession number: PRJNA982386.

The manuscript number allotted by the ethical committee is CSIR-NBRI_MS/2023/06/10.

SNJ and SVS conceptualized the idea of the project. NMB and HK performed field cultivation of cotton germplasms at TNAU, Coimbatore, Tamil Nadu, and PAU Faridkot, Punjab, respectively. BJ conducted the phenotyping and performed all the experiments, whereas GJT participated in phenotyping only. SS, BJ, SJ, and MI conducted all the computational and statistical data analyses. BJ wrote the first draft of the manuscript and SS helped into it. GJT, NMB, DA, SJ, DK, SVS, and SNJ revised the manuscript and provided various suggestions. BJ includes all the corrections and suggestions in the final manuscript. All authors contributed to the article and approved the submitted version.