Recombinant inbred lines derived from the intra-specific cross between G. hirsutum var. MCU5 and TCH1218 were used in this study. The F₁ was made at Tamil Nadu Agricultural University (TNAU), Coimbatore, in the year 2006 and forwarded to F₁₁ using the single seed descent method. In total, 181 RILs from the F₁₀ generation were employed in the construction of the genetic map. Data on PH (cm), number of sympodial branches, number of bolls, and boll weight (g) obtained from the phenotypic evaluation of F₇ to F₁₁ of those same RILs were used for QTL analysis. The female parent MCU5 is a multiple cross derivative cultivated widely in South India as it possesses medium staple fiber (29 mm) and can be spun up to 70s. Besides, it has a 34% ginning outturn and produces 1,850 kg seed cotton yield per ha; however, it cannot withstand water stress at the flowering phase (Boopathi et al., 2014). On the other hand, the male parent, TCH1218, has relatively better drought tolerance but low fiber yield and quality characters than MCU5 (Boopathi et al., 2014). Our historical breeding effort has shown that MCU 5 has long-staple cotton with a good yield, while TCH 1218 is a good combiner.

All the plant materials (181 RILs and two parental lines) were evaluated in randomized block design with two replications at different locations and seasons (the details of testing environments are provided in Table 1). The plant-to-plant distance was 45 cm, while the row-to-row distance was 90 cm, and 13 plants were maintained in each row. Regular crop husbandry measures were followed to ensure a healthy crop. To collect phenotypic data, three plants were selected randomly per replication on 115 days after sowing (DAS), and information was collected on different yield parameters from each experiment as detailed in Table 1. RILs were evaluated under field water stress conditions in two different water-stress environments: one with managed water stress (by withholding irrigation water after 45th day) at Maize Research Station (MRS), TNAU, Vagarai, and another under purely rainfed conditions at Cotton Research Station (CRS), TNAU, Veppanthattai and Agricultural Research Station (ARS), TNAU, Aruppukottai. Regular agricultural practices were used to control weeds, diseases, pests, and fertilizers were applied at sufficient levels to ensure that yield potential was not limited by any factor other than water. At MRS (T2 in Table 1), water stress was imposed by withholding water from 45th DAS, whereas the irrigated control plots received water at routine intervals. At CRS T3 and T7 (Table 1) and ARS T6 (Table 1), the RILs were grown purely under rainfed conditions, and there was no rain after 51 and 57 DAS at CRS during 2012 and 2018, respectively, and 58 DAS at ARS during 2018.

Standard descriptive statistical analysis to examine the significance of the difference in the investigated traits between two parents and among population and the estimation of Pearson’s correlation coefficient was done using Calculator Soup^®4 and Minitab^® 19⁵, respectively.

Genomic DNA was extracted from the young leaves of mapping population and parents using a modified cetyltrimethylammonium bromide (CTAB) method (Shukla et al., 2021). SNP genotyping was performed using the Gene Titan Multi-channel instrument (Thermo Scientific) facility at M/s Imperial life Sciences, New Delhi, with the Cotton 50K SNP array chip. DNA samples were first processed to the amplification, fragmentation, precipitation, and re-suspension and then hybridized to the chip using the Affymetrix reagents. The arrays were scanned with the Gene Titan MC and automated allele calling, and the quality assessment of called genotypes was done with Genotyping Console™ Software (GTC) with a new Axiom Genotyping Algorithm v1 (Axiom GT1). Consequently, the raw hybridization intensity data were processed for clustering and genotype calling with Affymetrix^® GTC (v4.2), and those data with Dish Quality Control (DQC) value <0.82 and call rate <0.97 were excluded from the further genotyping analysis. To finish, GTC was further processed using APT (v1.19.0) and classified those SNPs into six major classes. The raw genotyping data were again filtered by selecting only those co-dominantly segregating SNPs with <20% missing genotyping data (Hulse-Kemp et al., 2015). Subsequently, the genotyping data were transformed into a mapping data format. The probe sequences (10–20 bp overlapping to SNP) were also mapped on reference genome G. hirsutum cultivar TM-1 (ASM98774v1), and the SNPs were assigned to the particular chromosomes. The same data were used as anchoring information while preparing the input file for the mapping software while performing chromosome assignments.

A total of 522 polymorphic and segregating SNP markers were used to construct the genetic map using ‘IcImapping v4.0’ (Meng et al., 2015), and the genetic map distance (centiMorgan, cM) was calculated using the Kosambi mapping function. The threshold for logarithm of odds (LOD) score was fixed at five, revealing 26 linkage groups. The exact ordering of the SNPs across chromosomes was done using RECORD, and Rippling was also executed for fine-tuning the order of markers by the sum of adjacent recombination frequencies (COUNT) with a window size of five. If there were more than three consecutive adjacent markers in the genetic map with a significance level 0.001 < p < 0.05, it was taken as segregation distortion, and its distribution on the map was also analyzed (Coulton et al., 2020).

Composite interval mapping (CIM) with forward regression, in a window size of 10 cM and 5 background control markers at a walking speed of 2 cM, was employed for QTL analysis using Windows QTL Cartographer 2.5 (Wang et al., 2012). The threshold LOD of each trait QTL was determined with the trait threshold calculated with a permutation of 10,000, and multiple QTL peaks detected within 10 cM of each other were regarded as a single QTL. To calculate the effects of additive QTL in multiple environments, the MET functional module of QTL IciMapping v4.0 was employed by combined mapping analysis under multi-environment with 3.00 permutation tests’ LOD cut-off values.

All the SNP markers used to construct the genetic map reported in this study were aligned to the physical map of the upland cotton genome (TM-1 Genome NAU-NBI Assembly v1.1 and Annotation v1.1 database) through Basic Local Alignment Search Tool (BLAST). CIRCOS 0.66 with default parameters was employed to compare each investigated SNP marker’s genetic and physical positions collinearity. The recombination hotspot (RH) was estimated by inferring the recombination rate of investigated SNPs. If the genetic distance between adjacent SNPs was higher than 20 cM/Megabase, the genomic region between those markers was viewed as RH (Zhang J. et al., 2015).

Quantitative trait loci co-localized at least in two trials with ≥10.0 R² for at least one trait were depicted as major and stable QTL, and they were selected to scan for candidate genes. In general, a line with at least 10% higher variation than the parents for the target traits is considered an improved line; hence ≥10.0 R² is depicted as a major QTL in this study. The sequences of SNP markers flanking the confidence intervals (CIs) of the QTL were aligned back to the physical sequence of the upland cotton genome database (Zhang J. et al., 2015). Based on the position of these flanking markers, all the genes within the target QTL were identified as candidate genes. Gene ontology (GO) grouping was employed to categorize the identified candidate genes, and pathways correlated to the candidate genes were also discovered using Kyoto Encyclopaedia of Genes and Genomes (KEGG) analysis (Charmpi and Ycart, 2015).

The differences noticed for the investigated phenotypic traits between parents and the variation among the investigated RILs are provided in Table 2. As expected, the elite parent MCU5 performed well under the irrigated conditions, whereas TCH1218 had produced substantially good yield under both managed water stress and rainfed conditions by having higher values in the relative performance of investigated traits (significant at the 0.05 probability level) (Table 2). Transgressive segregation was observed in the RIL population for all the investigated traits under irrigated and water-stressed conditions. However, PH and boll weight have shown substantially skewed distribution as their skewness values were greater than +1 or lower than −1 (Table 2). All possible kinds of distribution, such as fairly symmetrical (skewness −0.5 and 0.5), moderately skewed skewness (skewness −1 and −0.5 or between 0.5 and 1), and highly skewed (skewness less than −1 or greater than 1) distribution of the traits examined in this study were noticed (Table 2). As the calculated kurtosis values for all the investigated traits were greater than zero (Table 2), it can be concluded that each trait distribution has a heavier tail (leptokurtic distribution). Substantial variations (minimum and maximum values of each trait) due to genotype differences were reported for all the investigated traits among the RILs. However, the relative proportions of variance varied from one trait to another, and there was low to moderate heritability noticed for the investigated traits (Table 2). The correlation between investigated traits measured from all the trials (T1–T8) was also evaluated (Table 3). Although not in all the trials, the majority of correlations have shown a highly significant positive correlation between yield (BN and BW) and growth parameters (PH and SYM) both under irrigated and water-stressed environments, indicating that genetic improvement BN and BW would likely be accompanied by the improvement of PH and SYM. However, significant negative correlations were also observed (Table 3), which showed that the simultaneous improvement of investigated traits could not be done for all the cotton-growing regions. In general, the correlation between the investigated traits increased dramatically under rainfed conditions (T3 and T7); however, mixed trends were noticed in one rainfed environment, T6 (Table 3). On the other hand, under managed water stress conditions (T2), a reduced level of correlation between the measured traits was noticed when compared with the traits measured from irrigated conditions (T1) (Table 3).

Affymetrix’s Axiom custom designed Cotton 50K SNP array (unpublished) was used to the SNP typing of both parents and 181 RILs. This resulted a total of 2,457 polymorphic SNPs between the two parents and those polymorphic markers were analyzed for their segregation in the mapping population. The initial attempt of binning of these 2,457 loci by excluding those SNPs with more than 20% missing data (totally there were 17 such markers) resulted in 2,440 polymorphic SNPs, which were further analyzed with χ² test to determine if the allele frequency was deviated from the expected segregation ratio (1:1). Among the selected 2,440 polymorphic loci, 1,913 loci showed segregation distortion (p < 0.05), and only 527 loci followed the perfect Mendelian ratio. It was also noticed that among 1,913 distorted loci, 1,186 were in favor of TCH1218 alleles, and 727 loci were biased toward MCU5 alleles. Further, it has been found that among the 527 loci that were segregated in the RILs, and five were found to be duplicated.

Therefore, only 522 polymorphic SNPs were taken for the genetic mapping, out of which 504 were mapped on 26 AD linkage groups (LGs) or cotton chromosomes (as 18 SNPs were unlinked), and the resultant genetic map represented the genetic span length of 4,416 cM. The At sub-genome had 204 SNP markers covering 1,830.03 cM with an average genetic distance of 8.97 cM between adjacent loci, while the Dt sub-genome possessed 300 SNP markers that span 2,585.97 cM with an average of 8.62 cM between successive loci. Overall, the average marker distance among the 26 chromosomes was 8.76 cM (Table 4). The longest and most dense chromosome found in this study was chromosome 18 (395.48 cM) with 79 markers, and the shortest and most sparse was chromosome 02 (85.03 cM) with 2 markers (Table 4 and Figure 1). In the genetic map of chromosome 2, there were only two markers mapped, and 19 markers showing segregated distorted (SD) were removed from the map. Thus, there is a huge gap in chromosome 2. Besides chromosome 2, other small gaps were recorded in this genetic map due to the removing of SD markers. Therefore, the marker interval in each chromosome was ranged from 7.9 to 85.0 cM (Table 4 and Figure 1).

The quality of the genetic linkage map was evaluated using parameters, such as the segregation distortion of mapped markers, gaps in the map, and collinearity between the linkage map constructed in this study and the reference physical map of TM-1⁶. All the 504 polymorphic SNPs found in this study were aligned to the G. hirsutum var. TM-1 reference genome using BWA in Galaxy and a high level of collinearity was shown by the SNP genetic map reported in this study with the physical map of the TM-1. It was revealed that there were no chromosomal rearrangement in the mapping population with context to the parents. The CIRCOS plot further confirmed the accuracy of the genetic maps (Figure 2). All the synteny blocks of each chromosome corresponded to the physical map and no single synteny block mismatched with the corresponding chromosome. Thus, it can be inferred that there were no chromosomal rearrangements. The sequence-based reference genome of 1,027.9 Mb corresponded to our SNP based genetic map of 4,416 cM and our genetic map represented 53.13% of the total length of the sequence-based physical map (Table 5). A good collinearity was revealed by all the linkage groups with the physical reference map (Figure 2). The coverage of individual chromosome in the constructed genetic linkage map of our study ranged from 0.64 to 99.11% of the physical reference map (Table 5). Chromosome 26 was least covered whereas, chromosome 3 was comprehensively covered. Though least covered chromosomes in the genetic map would provide little information during QTL analysis, it would be useful to have such preliminary data to proceed further.

As the variation found in the investigated traits was normally distributed (Table 2), it was predicted that multiple genes might control these traits. There was low to moderate heritability noticed for the investigated traits (Table 2). Thus, higher similarity between parents and offspring for the investigated trait was not expected making the QTL study more complex, and interesting for significant environmental influence. In such context, the identification of co-localized QTLs across the seasons and environments was the only means of solid proof for the presence of QTL linked to the investigated trait. A total of 53 QTLs for four investigated traits identified in this study, and among them, 21 and 32 were found to be At sub-genome- and Dt sub-genome-specific, respectively (Supplementary Table 1). The LOD values of all the identified QTLs ranged from 2.0 to 18.5. There were eight water stress-specific QTL further linked to different investigated traits. Though there are 8 QTL identified for PH (6 under irrigated conditions and two under water stress conditions), only two QTL had shown more than 10% phenotypic variation (Supplementary Table 1). Among them, qPH-c11-1 explained 38.41% phenotypic variation (PVE) and was detected under irrigated conditions; whereas qPH_S-c16-2 was detected under rainfed conditions had 19.39% PVE. A maximum number of 36 QTL were detected for the number of sympodial branches (SYM), and among them, 25 QTL were detected under irrigated conditions, and 11 were detected under managed water stress/rainfed conditions. In these, QTL detected on chromosome 18 (qSYM-c18-2) had maximum phenotypic variation under irrigated conditions (87.83%) with QTL LOD of 8.6 and a QTL observed on chromosome 21 (qSYM_S-c22-2) had accounted for a maximum of phenotypic variation under rainfed conditions (88.01%) for SYM with QTL LOD value of 3.6. Seven QTL were also identified for boll number (BN) in which all the QTL were found under irrigated conditions. Among, qBN-c26-1 exhibited 83.61% phenotypic variation. Only two QTL were found to be linked with boll weight (BW) under irrigated conditions, and no QTL was observed for this trait under rainfed conditions. The maximum PVE (81.34%) for this trait was reported for qBW-c05-2.

Twenty QTL had a positive additive effect for the investigated traits suggesting the independent effects of alleles on the trait (Supplementary Table 1). Several QTL were found to have a dominant positive effect (Supplementary Table 1). Among the 53 QTL, 25 QTL were derived from MCU5, which exhibited a positive additive effect, and 28 QTL were derived from TCH1218, which exhibited a negative additive effect. As expected, all the QTL detected under water stress conditions were contributed by the tolerant drought parent, TCH1218 (Supplementary Table 1) except qPH_S-c16-1, qSYM_S-c18-1, and qSYM_S-c18-2, which were contributed by drought susceptible parent, MCU5.

It was found that nine QTL that were not only co-localized for more than one investigated trait but also stable (identified at least in two seasons/locations) and major QTL with >10% PVE for at least one trait that colocalized in that region (Table 6). In addition, the maximum number of seven QTL was clustered on chromosome 22 and a minimum number of two QTL was found on chromosomes 5, 21, and 24 (Table 6 and Figure 1).

Further, joint multi-environmental QTL mapping effort to estimate the impact of QTL × Environment interaction resulted in the identification of 11 additive QTL (Table 7). Eight QTL associated with SYM were identified with LOD ranging from 3.02 to 5.70. The contribution rates of interaction among eight additive and the environment QTL ranged from 0.68 to 3.33%. Two QTL associated with PH were identified with LOD ranging from 3.13 to 3.53. The contribution rates of interaction among two additives and the environment QTLs ranged from 2.16 to 4.20%. A single QTL associated with BN was identified with LOD 4.36 and the contribution rates of interaction was 4.08% (Table 7). However, no QTL for BW was detected in this joint multi-environmental QTL analysis.

Large array of candidate genes were identified within the QTL reported in this study (Supplementary Table 2), and they were compared and annotated with Xu et al. (2020), where meta-QTL analysis along with transcriptomic approach utilized for the identification of candidate genes related to fiber quality in upland cotton. Though there were genes specific to abiotic stress responses and productive traits, large numbers of genes identified in this study warrant the use of additional markers to fine map these QTL and identify precise genes involved for the target traits. For example, for the PH QTL under irrigated conditions, qPH-c11-1, there were 32 genes, whereas under water stress conditions, the QTL identified for PH, qPH_S-c16-2, possessed 94 candidate genes (Supplementary Table 2). Similarly, the QTL identified for a number of sympodial branches underwater-stressed environment, qSYM_S-c24-1 harbored 169 genes and QTL for BN under irrigated conditions, qBN-c26-1 contained 733 genes, and the BW QTL under irrigated conditions, qBW-c05-2 had 2,074 genes.

The GO analysis used all the identified candidate genes to identify potential biological functions and grouped under three main GO categories: biological process, molecular function, and cellular component (Supplementary Table 3; Sheet “GO”). Within the biological process category, there were 12 sub-categories, such as abiotic stress-specific, such as response to freezing, and response to biotic stimulus sub-category possessed the maximum number of 15 candidate genes. ADP binding with 21 genes was the main sub-category in the molecular function category with 14 principal sub-categories. Finally, in the cellular component category, there was one sub-category (vacuolar membrane) with 11 genes (Supplementary Table 3; Sheet “GO”). On the other hand, KEGG analysis resulted in only one category, namely, plant-pathogen interaction with four candidate genes (Supplementary Table 3; Sheet “KEGG”).

The lack of large numbers and the segregation distortion of polymorphic markers are the main leading cause for the availability of a limited number of high-density linkage maps in cotton (Zhang T. et al., 2015). As a result, cotton genetic maps with SSR or AFLP markers often shows large gaps, poor marker density, and low marker coverage, whereas SNPs have shown their utility in constructing a high-dense genetic map as they are distributed throughout the genome (Cai et al., 2017; Liu et al., 2018). Among various strategies used in SNP genotyping (such as, genotyping by sequencing and whole-genome resequencing), SNP arrays have been shown to be simple and useful in the development of dense genetic maps in cotton (Li et al., 2016; Tan et al., 2018; Kumar et al., 2019). For example, in the present study 2,457 polymorphic markers between the parents of an intraspecific cross of upland cotton were identified using Cotton SNP50K array.

Progenies resulting from various cotton cross combinations frequently have unequal allelic distribution (Shao et al., 2014). Because of marker and population types, preferred fertilization, gametic combinations, genetic drift, and environmental variables, segregation distortion can vary across and within species (Shen et al., 2007). Despite the fact that a considerable number of polymorphic markers were found in this study, a large proportion of them (1,953, or 79.48% of total polymorphic loci identified between the parents) were removed due to segregation distortion. As a result, the number of loci in this study that can be adequately mapped has been drastically decreased. Cotton has also shown a low degree of mapping effectiveness for discovered polymorphic markers (Lacape et al., 2010; Yu et al., 2013). On the other hand, the constructed genetic map in present study did not have any segregation distortion region.

Thus, after a stringent screening of missing data, segregation distortion, and similar loci, a total of 504 markers were mapped only on 26 chromosomes. As shown in Table 4, the highest number of markers (300 SNPs) was distributed on Dt sub-genome when compared with At sub-genome (204 SNPs). It was inferred that the genetic map constructed in this effort had highest collinearity with the physical map (all the 504 polymorphic SNPs were aligned perfectly; Figure 2) and had covered 53.13% of the physical map (TM-1 Genome NAU-NBI Assembly v1.1 and Annotation v1.1 database; Table 5), which highlighted that there was no chromosomal rearrangement in the intraspecific genetic map constructed in this study.

Several small effect QTL are involved in confirming the drought tolerance in cotton, and each loci represents hundreds of genes, which are the genetic basis for expressing an extensive array of water stress responses in the form of morpho-physiological traits (Abdelraheem et al., 2021). As only a few studies have focused on the QTL mapping of drought tolerance in cotton (reviewed in Mahmood et al., 2020), the molecular breeding of cotton for the genetic improvement of water stress tolerance is considered as a challenging task. Further, such QTL have large intervals (Mahmood et al., 2020; Abdelraheem et al., 2021).

This study found large intervals in the identified QTL, mainly attributed to the poor polymorphism rate detected between the two parents used in this study. Further, due to the segregation distortion, the number of mappable markers was reduced drastically from 2,457 to 504, which resulted in due to the segregation distortion. This drastic reduction of mappable markers resulted in gaps in the genetic map constructed in this study, reporting in G. hirsutum (Wang et al., 2006; Wu et al., 2009; Sun et al., 2012) and even in interspecific crosses of Gossypium spp. (He et al., 2007; Lacape et al., 2010).

Though there are different strategies (such as, genotyping by sequencing) to increase the efficiency of genetic mapping and QTL analysis through fine mapping, it would be desirable to employ high-density SNP chips as it enables speedy and automated detection of a high rate of polymorphisms across cotton accessions within a short span of time. Therefore, the development of new SNP chips by using an additional set of SNPs, detected on the water-stress responsive candidate genes identified in this study, would further enhance the efficiency of cotton molecular breeding for abiotic stress resistance.

Compared with the published reports on QTL linked to the drought tolerance and fiber yield and quality traits (Mahmood et al., 2020), a relatively low number of QTL were detected in this study. It indicated that the traits examined in this study were under strong genotype by environment (GxE) interaction and low heritability. Despite this, nine common QTL were detected across the seasons and locations, which pointed out that the collected phenotypic data are reliable even though there were GxE interactions. Further, a significant correlation among those traits and a substantial amount of heritability and additive effect estimates for most of the QTL identified in this study showed that the selection of investigated traits for drought tolerance improvement in cotton under filed water stress conditions would be more efficient. Similar evidence for selecting productive traits under abiotic stress environments in cotton has already been reported (Abdelraheem et al., 2018, 2021).

Surprisingly, both parents contributed to the additive effects of water stress tolerance alleles (Table 6). This demonstrated that even the inferior parent (MCU5) might contribute to the establishment of drought tolerance, a feature that has previously been demonstrated (Paterson et al., 2003; Ulloa et al., 2005; Shen et al., 2007). It is imperative to identify the consensus and hot spot QTL across the testing environments (as shown below), which will facilitate the introgression of causal genomic regions that impart at least a minimum increase in the productivity under water stress. Though this study had identified several QTL (Supplementary Table 1), only those QTL that have shown their potential in molecular breeding have been selected (Table 6) for further analysis. This study had identified a total of 11 QTL in joint multi-environmental analysis. Compared with individual analysis (Supplementary Table 1), the joint multi-environmental analysis (Table 7) identified a smaller number of QTL due to both QTL additive and QTL-environment interaction effects. However, the former analysis estimates the only additive effect of QTL but not the environmental influence. Thus, this effort has helped to get additional information on the influence of the environment on the expression of QTL under different water-limited environments.

Further, this study attempted to identify hotspot QTL (which is defined here as a cotton genomic region), where QTL were discovered for more than one trait that was investigated across the seasons and locations expressing with more than 10% phenotypic variations (in this study) and genetic backgrounds (elsewhere) were colocalized in the same genomic segment. Nine hotspot QTL were identified in this study that were located on chromosomes 5, 7, 9, 11, 15, 18, 21, 22, 24, and 26 (Table 6). Among them, a region on chromosome 22 flanked by AX-182254626–AX-182264770 with a span length of 89.4 cM, is co-localized for a maximum number of major and stable QTL (7) linked to the number of sympodial branches both under irrigated and water stress conditions (Figure 1). This region was specifically detected for fruit branch node number by Hulse-Kemp et al. (2015). Another notable region was on chromosome 18 that span a length of 61.4 cM flanked by the SNPs, AX-182256182 and AX-182255051. This region was associated with five QTL linked to the number of sympodial branches both under irrigated and water-stressed conditions in this study. Similarly, a QTL cluster was found on chromosome 5 for BW and BN in this study and similar kind of correlation of this chromosomal region with several productive traits has been reported earlier (Abdelraheem et al., 2018).

Except for chromosomes 5 and 26, almost all of the hotspots discovered in this study were connected to the number of sympodial branches. Cotton bolls can form at the nodes of sympodia that develop from the monopodia, as well as at the main or secondary vegetative axis sites. Orderly development of bolls on sympodia has long been favored (McClelland, 1916), and plant design with 0% monopodial growth would be more ideal, favoring efficient mechanical harvesting and improving ultimate yield. Since only major, stable, and consistent QTL were reported in this study, it is likely that all of the QTL listed in Table 6 could be potential assets for consensus mapping across genetic backgrounds for their validation and identification of most useful QTL hotspots, map-based cloning, and MAS for productive traits in upland cotton with high predictability. Among the three PH QTL chosen in this investigation (Table 6), a QTL on chromosome 11, qPH-c11-1, has demonstrated the most phenotypic variance (38.41%), and this area has been proven to host PH QTL under irrigated (Jia et al., 2016) and salt stress environments (Abdelraheem et al., 2021). Optimum plant height in cotton is the most desirable trait in the terms of mechanical harvesting, besides its direct relation with biomass that influences the final yield (Shang et al., 2015).

It would be an elaborated list to identify candidate genes for the QTL reported in this study, as it has been estimated that ∼4,500 cM of G. hirsutum genome consisted of 70,478 predicted protein-coding genes (Zhang T. et al., 2015). Similarly, the large array of candidate genes was found in this study (Supplementary Table 2). Despite this extensive list, identifying candidate genes in the target hotspot QTL can open up new avenues in understanding the molecular basis of drought tolerance in cotton.

For example, qSYM-c22-1, a QTL on chromosome 22 flanked by AX-182254626 and AX-182258142, was found to be co-localized for the number of sympodial branches both under irrigated and water-stressed conditions, and it harbors 1,046 genes (Supplementary Table 2). Among them, the notables are GDSL esterases/lipases, involved in the regulation of plant development by the synthesis of secondary metabolites in response to the biotic and abiotic stresses (Chepyshko et al., 2012), zinc finger protein CONSTANS, regulating flowering under normal and stress conditions (Putterill et al., 1995), abscisic acid receptor PYL8, expressed in response to dark-induced leaf senescence (Lee et al., 2015), choline mono oxygenase, conferring abiotic stress tolerance by synthesizing glycine betaine (Bao et al., 2011), potassium transporter 2, copper transport protein ATX1, and peroxidase 52, contributing significantly for abiotic stress resistance (Zhang et al., 2021) and several network of reactive oxygen species (ROS) genes, those have played critical role in both abiotic stress response and fiber development (Xu et al., 2019). Similar kinds of candidate genes were identified for another hotspot QTL (Supplementary Table 2). Thus, this study provides a promising lead to precise the QTL mapping consisting of functional SNPs derived from those candidate genes identified in the target QTL by developing a new SNP chip.