Plant materials used in this study were developed from reciprocal crosses between Gossypium hirsutum acc. Acala Maxxa and G. barbadense acc. Pima S6 (both inbred lines). These genotypes have been extensively used to produce several molecular tools and resources including BAC libraries and Illumina genome sequences. Reciprocal advanced backcross populations were developed by first crossing the parents in a two-way cross (Acala Maxxa (♀) × Pima S6 (♂) – GH background; and Pima S6 (♀) × Acala Maxxa (♂) – GB background), then independently backcrossing F₁ plants to the respective female parents to create 300 to 400 BC₁ progenies for each cross. The backcrossing scheme utilized single seed descent to generate each subsequent generation (Adhikari et al., 2023).

After five generation of backcrossing, 179 BC₅F₁ plants from the GH background and 190 BC₅F₁ plants from the GB background were self-pollinated and a total of 369 BC₅F₂ families were grown at Iron Horse Farm (IHF), Watkinsville, Georgia in 2019 and 2021 and at Southwest Georgia Research Station, Plains, Georgia in 2021 under cultural conditions consistent with commercial irrigated cotton production. Individual BC₅F₂ plants that contained only one introgressed segment from the donor parent were deemed as NILs. Selfed seeds of these NILs (BC₅F₃ seeds) were grown at IHF and Plains in 2021 under cultural conditions consistent with commercial irrigated cotton production.

The genomic composition of the BC₅F₁ plants was inferred based on genotyping by sequencing (GBS). DNA was extracted from the parents and 369 BC₅F₁ plants using a scaled-down version of a published CTAB protocol (Paterson et al., 1993). A total of five multiplexed GBS libraries were constructed according to (Andolfatto et al., 2011) wherein the DNA were double digested with HinP1I-HaeIII enzymes. The libraries were sequenced on Illumina MiSeq (in-house) with 75 bp single end reads (SE75). The TASSEL5 GBSV2 pipeline was used for sequence data processing and genotype calling (Glaubitz et al., 2014). Reads were aligned to G. hirsutum acc. TM-1 (Zhang et al., 2015) using Burrow-Wheeler Alignment (bwa) and exported to variant call format (VCF). To minimize sequencing errors, only the first 64 base pairs were used to map reads to the reference genome. Filtering of the VCF was done for bi-allelic SNPs using Fisher’s exact test with a threshold P value < 0.001, considering that true variants should represent biallelic homozygous state for inbred accessions. Genotypes for lines in G. hirsutum background were called together with one another; and those for lines in G. barbadense background cross were also called together. The SNPs were filtered for MAF > 0.01, missing scores<30% and heterozygosity<10% at the population level. The retained SNPs were imputed using the Fast Inbred Line Library Imputation (FILLIN) pipeline available in TASSEL5 GBSv2 (Swarts et al., 2014).

The genomic composition of BC₅F₂ plants were inferred based on targeted microsatellite (SSR) genotyping of the introgressed chromosomal segments identified in their respective BC₅F₁ parents. At least two (and at most four) SSR markers were used to verify most of the introgressed regions while for small introgressions only one SSR marker was deployed. A total of 852 polymorphic SSR markers ( Supplementary Table 1 ) spanning the introgressed regions were derived from several published genetic maps of crosses between G. barbadense and G. hirsutum stored in the CottonGen SSR database (https://www.cottongen.org/data/download/marker). A total of 47 candidate SSRs were monomorphic in our lines and discarded, as were 23 with ambiguous bands. Among the 8364 BC₅F₂ individuals planted in 2019, the remaining 782 SSR markers were used to genotype 5315 plants (from BC₅F₁ parents carrying 2 to 5 introgressions) for the presence (or absence) and dosage (homozygous vs heterozygous) of the respective introgression/s. Individual BC₅F₂ plants that were verified by SSR markers to carry only one introgression (homozygous or heterozygous) from the donor parent were deemed as NILs.

Two replications of each BC₅F₂ family and NIL were planted in a randomized complete block design (RCBD) in three environments (IHF-2019, IHF-2021 and Plains-2021). Six replications each of the two parents were included in all three environments. Fiber samples were collected by harvesting 25 bolls from each plot, ginned in a laboratory gin, and evaluated by HVI (Cotton Incorporated Textile Service Laboratory, Cary, NC) to extract values for fiber elongation (ELO).

All statistical analyses were conducted in R programming language. Single marker analyses were done in R/qtl (Broman and Sen, 2009). The significance threshold was set to LOD of 3, to mitigate the multiple-comparison problem. Filtration of significant markers adopts the method proposed by Szalma et al. (2007). If several markers on the same introgressed segment show significant association with phenotype, the most significant one was reported. For the co-segregation of multiple introgressions, the QTL location is examined as follows. First, if multiple families show significance for the trait and carry overlapping introgression, the introgression is considered to carry QTL. Second, if the co-segregation of introgressions is in single families, the most significant introgression is considered to carry QTL.

Phenotypic variance explained by each locus was reported by taking the most significant marker as independent variable and phenotypic value as dependent variable in R (R Core Team, 2016). Additive effects were estimated by half the difference of phenotypic values between the lines carrying the homozygous introgression and lines not carrying the introgression. Dominance effects were estimated by the difference of phenotypic values between the lines carrying the heterozygous introgression and the remaining lines that do not carry the introgression. If multiple or overlapping introgressions were present at both homozygous and heterozygous state, the estimation of additive effects utilized the lines carrying the introgression at homozygous state only and, the estimation of dominance utilized the lines carrying the introgression at heterozygous state only.

Gene actions for QTLs were determined by calculating the degree of dominance (absolute values) for every QTL that has both additive and dominance effects. The degree of dominance is the ratio of dominance effect to additive effect (d/a) of the QTL and based on this ratio, gene action of the QTLs can be categorized as (i) additive (0<d/a< 0.2) (ii) partially dominant (0.2< d/a<0.8) (iii) dominant (0.8< d/a< 1.2 and (iv) over-dominant (d/a > 1.2) (Edwards et al., 1987). QTLs with dominant and over-dominant effects are considered to have heterotic effect or heterozygous advantage.

‘Common QTL’ refers to cases in which the same marker is detected in two reciprocal populations, or two different markers are statistically significant for the same/overlapping QTL(s)/introgression(s) in each of two populations. The hypothesis that two sets of QTLs are randomly distributed across the entire genome can be falsified using the hypergeometric probability function (Paterson et al., 1995; Feltus et al., 2006). In the model, p is the probability of non-random correspondence of QTLs being compared; n is the number of comparable intervals which is calculated by dividing the total genome size by average introgression size in both populations; m is the number of common QTLs; l is the number of QTLs in the GH background; s is the number of QTLs in the GB background. The same model was also adopted to detect correspondence between QTLs reported in this study with those previously published. In this case, l is the total number of QTLs identified in the larger sample (study reporting higher number of QTLs) and s is the number of QTLs identified in the smaller sample.

In silico annotation was performed on the identified QTLs to look for candidate genes related to fiber elongation in cotton. For each QTL identified in the study, the genomic region spanning 1 Mb on each side of the most significantly associated marker was used for in silico analysis. The DNA sequence from this tightly linked region was used to look for G. hirsutum genes in the CottonGen database and these genes were then analyzed for biological functions, with particular focus on fiber growth and development.

Detailed description of the genomic composition of the reciprocal advanced backcross lines (BC₅F₁ populations) is published (Adhikari et al., 2023, 2025). The NIL population in the Acala Maxxa (G. hirsutum) background consisted of 397 individuals with only one chromosomal segment introgressed from the donor parent (G. barbadense ‘Pima S6’), while that in the Pima S6 background consisted of 423 individuals. In total, these lines covered 78.72% and 71.48% of the donor genome respectively in the Acala Maxxa and Pima S6 backgrounds. Single introgressed segments were identified for all chromosomes except chromosome 24 in the Acala Maxxa background and chromosome 25 in the Pima S6 background. Individual lines contained an average of 0.904% of the donor genome, ranging from 0.07% (1.37 Mb) to 2.92% (56.65 Mb) in the Acala Maxxa background while in the Pima S6 background, individual lines contained an average of 0.97% of the donor genome, ranging from 0.15% to 3.07%.

The phenotypic performance of the two parents, reciprocal backcross populations, and reciprocal NIL populations for fiber elongation is shown in Figure 1 . The distribution of fiber elongation was approximately normal (Shapiro and Wilk test; p > 0.05) and typical of quantitative inheritance. Pima S6 outperformed Acala Maxxa in all three environments ( Supplementary Tables 2 , 3 ; Figure 1 ). Both advanced backcross and NIL populations in the GB background performed better than in the GH background. Transgressive segregation is seen for all populations across all environments tested. Transgressive segregants outperforming both parents were identified in both population types and both backgrounds ( Figure 1 ).

To identify the effect of genotypes and environment in the overall performance of the advanced backcross populations and the NILs, we conducted analysis of variance (ANOVA) treating all variables as fixed factors. Results showed significant effects of both genotype (GEN), and genotype-by-environment (GXE) ( Supplementary Table 4 ). GEN captured the most variation in both population types and GXE also captured significant variation, precluding the use of combined phenotypic values in identification of QTLs. Thus, marker-trait association and identification of fiber elongation QTLs was performed separately for each environment tested.

A total of 36 marker-trait associations were identified for fiber elongation. Phenotypic variances explained by these QTLs ranged from 3.01% to 18.41% ( Tables 1 , 2 ). Among the 36 QTLs identified, 11 were of large effect, explaining > 10% of the total phenotypic variation while the remaining 25 were small effect QTLs (explaining<10% of total phenotypic variation).

In the GH background, a total of 21 QTLs were identified for ELO in the three environments tested, explaining 3.01% to 18.41% of total phenotypic variation ( Table 1 ). These QTLs were distributed over 15 chromosomes, with chromosome 12 carrying the most QTLs. A total of 11 of these QTLs were small effect and the remaining 10 were large effect. Among the 21 QTLs identified, eight were identified at IHF in 2019, nine at IHF in 2021 and four at Plains in 2021. All 21 QTLs identified in the GH background increased fiber elongation consistent with the parental phenotypes. The ELO QTL on chromosome 12 was identified in all three environments tested and is a potential candidate for fine mapping and QTL identification studies. Another QTL on chromosome 15 was identified at IHF in 2019 and at Plains in 2021, making it a priority for further study of its effect on fiber elongation.

In the GB background, a total of 15 QTLs were identified for ELO, explaining 3.67% to 12.03% of the total phenotypic variation ( Table 2 ). Among the 15 QTLs identified for ELO, two were identified at IHF in 2019, seven at IHF in 2021 and six at Plains in 2021. Of the 15 QTLs, only on chromosome 3 was a major effect QTL and the remaining 14 were small effect QTLs. Five of these QTLs decreased fiber elongation while 10 increased fiber elongation. One QTL on chromosome 6 and one on chromosome 9 were each identified in two environments (IHF 2021 and Plains 2021).

Among the 36 QTLs identified for fiber elongation, 22 were identified in the At subgenome and 14 were identified in the Dt subgenome. In both backgrounds, At subgenome carried more QTLs for fiber elongation than Dt subgenome (12 vs 10 in the GH background and 9 vs 5 in the GB background). QTLs for fiber elongation were identified in all chromosomes except chromosomes 11, 13 and 19. In the GH background, fiber elongation QTLs were identified in 15 of 26 chromosomes while in the GB background, they were identified in 11 chromosomes.

A total of 11 QTLs showed transgressive effects for fiber elongation, i.e. which were opposite of what would be predicted by comparison to the parental phenotypes ( Table 3 ). Five of the 11 transgressive QTLs were identified in the Acala Maxxa background while six were identified in the Pima S6 background. Except for EL08.1 which showed fiber elongation values lower than both the parents, all other transgressive QTLs outperformed the parents in both backgrounds. Two QTLs in Acala Maxxa background and one in Pima S6 background had average fiber elongation values above 7.0, far exceeding the performance of the better parent (Pima S6), thus, indicating the potential for their use in breeding superior cultivars.

Among the 36 QTLs identified in the study, three showed heterotic effects (d/a >|1.2|, see methods section for details). All three QTLs showed overdominant gene effects (d/a: 1.49 for EL07.1, d/a: 4.02 for EL05.1 and d/a: -5.46 for EL10.1; Tables 1 , 2 ). It is well known that the F1 hybrids between these two species have striking fiber quality, and are directly used as cultivars in regions such as India where hybrid seed production costs are low. Further evaluation of salient lines carrying alleles such as these to validate their effects in multi-environment and multi-location trials may elucidate the genetic underpinnings of such hybrid advantages.

NILs and advanced backcross populations in both backgrounds showed average phenotypes consistent with their recurrent parent ( Figure 1 ). Albeit the average performance of these two populations behaved like their recurrent parents, which is expected given the exceptionally large proportion of their genome coming from the recurrent parent, a large amount of variation was observable (Adhikari et al., 2023, 2025). The presence of transgressive segregants in both directions suggests that the chromatin segments introgressed from the donor parents have effects that could significantly alter the performance of individual lines.

The total number of QTLs identified for the three fiber quality traits differed little in the two backgrounds (21 in GH vs 15 in GB). While identification of similar total numbers of QTLs in the reciprocal backgrounds might suggest the involvement of similar number of genes in controlling these traits, many other findings indicate the effect of species background on effects, locations, and patterns of these QTLs. For example, all 21 QTLs identified in the GH background had positive effects, i.e. introgression from GB into GH always increased fiber elongation. While GB is known to have superior fiber elongation than GH and has the potential to transfer favorable alleles to GH background, it is an interesting observation to see that all the QTLs conferred by PimaS6 to the GH background increased fiber elongation. Observations like these, where alleles contributing favorable alleles are spread throughout the genome, would be suggestive of utilizing a comprehensive approach like genomic selection (rather than targeted selection) when breeding for such traits that have favorable alleles spread throughout the genome.

GH is usually not superior to GB in terms of fiber quality, and for this, it has rarely been used as a donor parent in trait mapping studies (Liu et al., 2023). The reciprocal nature of our populations allowed us not only to study GH as a recipient parent of favorable alleles from GB, but also to explore GH itself as potential source of novel alleles to the GB background. Our results suggest that GH can also contribute favorable alleles to GB. Of the 15 QTLs identified in the GB background, nine (more than half) had positive effects in increasing fiber elongation.

A perplexing observation, however, is the complete lack of reciprocity of QTLs identified with opposite phenotypic effects at corresponding locations in the reciprocal genetic backgrounds. In principle, one would expect the majority of QTLs to show such reciprocity if alternative alleles at a QTL show additive or dominant-recessive effects. Limited correspondence of identified QTLs in the two backgrounds could be a result of several factors. From first principles, noting that these lines covered 78.72% and 71.48% of the donor genome respectively in the Acala Maxxa and Pima S6 backgrounds, only a little more than half the genome is actually tested in both backgrounds. Further, the small phenotypic effects of most of the identified QTL, increases the likelihood that one or both members of a reciprocal pair elude detection (Broman, 2001). Another intriguing factor that could account for some failures to identify correspondence of QTLs in the reciprocal backgrounds, especially in advanced backcross lines with multiple introgressed chromosomal segments, is epistasis. Interaction between introgressed loci might result in underestimation of their effects which might have resulted in some QTLs failing to reach the biometric thresholds required to declare them as QTLs per se. The widespread observation that fiber quality parameters generally have high heritability (Fang et al., 2014) suggests a limited role of epistasis, but it could contribute to failures to identify reciprocal QTLs with relatively small effects (Chandnani et al., 2018). This might be the case here as most of the QTLs detected in this study show low genetic contribution to the total phenotypic variation explained by the phenotype.

While it is no longer controversial that an appreciable portion of variation in fiber quality derives from the D subgenome (Jiang et al., 1998) from an ancestor that does not produce spinnable fiber, the relative roles of the two cotton subgenomes vary by trait and population. With regard to the distribution of QTL in A and D subgenomes, many prior studies have concluded that more QTL for growing period, yield and fiber quality were distributed in the Dt subgenome than in the At subgenome (Rong et al., 2007; Said et al., 2013; Wang et al., 2013; Adhikari et al., 2017; Chandnani et al., 2018). However, Shen et al. (2005) and Lin et al. (2005) reported that more QTL for fiber quality were located in the At than the Dt subgenome. In the current study, we identified more QTLs in the At than the Dt subgenome. For NILs and advanced backcross populations in the GH background, the subgenomic affinities of QTLs were not significantly different (p>0.05) but nominally agreed with most previous findings, with the At subgenome harboring three more QTLs (12 vs 10) than the Dt subgenome. In the reciprocal background, however, a significantly higher (p<0.05) number of QTLs were identified in the At subgenome (9 vs 5).

While different loci associated with the same trait under various environments might suggest interaction between genotype and environment, QTLs being detected across environments might indicate environmental stability. Although most QTLs identified in this study were single-environment QTLs, seven genomic regions (on chromosomes 2, 6, 9, 12, 15 and 18) were consistently associated with ELO in two different environments ( Tables 1 , 2 , Figure 2 ). Since NILs facilitate pinpointing genomic regions affecting phenotypes to a rather small region around the causal mutation without accumulating a lot of extraneous variation, the identification of these seven stable regions for QTL in the cotton genome warrants more indepth look into these regions for genes controlling the trait. While this investigation is not in the scope of our current study, it will certainly help the scientific community to streamline efforts in these regions to identify and verify candidate genes controlling fiber elongation.

The genetic composition of our experimental populations provides a platform to identify novel small-effect QTLs in addition to major QTLs. Since NILs serve as a resource to not only identify marker trait associations but also an important tool to verify QTLs previously identified, mostly those using early generation populations, the correspondence identified here could be used as a means of validation of previously published QTLs. Given that there are hundreds of studies reporting fiber quality QTLs and owing to the genetic structure of our experimental populations, we mostly compared our results with populations of similar genomic composition (Chee et al., 2005b; Draye et al., 2005; Yu et al., 2013; Brown et al., 2019). We also performed elaborate statistical comparisons with other previous reports on the correspondence of QTLs identified in our study.

In a comprehensive study of biometric parameters of QTLs affecting fiber elongation using a backcross-self approach, (Chee et al., 2005a, b) identified a fiber elongation QTL (EL05.1) on chromosome 5 with the nearest locus being RFLP marker pAR206b, in the same introgressed segment where we identified the most significant marker for our ELO QTL (qEL05.1) in the GH background, both GB alleles increasing fiber elongation. In the same study, another fiber elongation QTL EL23.1 was identified on chromosome 23 with the nearest marker being an RFLP marker pAR209, near which we also identified a fiber elongation QTL (EL23.1) – albeit with a different phenotypic effect. While the fiber elongation QTLs identified on this chromosome in the previous study decreased fiber elongation, the one identified in our study increased fiber elongation. In both studies, the donor parent was the same (Pima S6) while the recipient parents were different (G. hirsutum cv. Tamcot 2111 in the previous study), implicating variation among cultivars within the same species as a possible cause of this difference.

While correspondence of individual marker-trait associations may reflect, for example, a gene for which the recurrent parent, Acala Maxxa or Pima S6, has a rare allele, non-random patterns of association across an entire genome can reflect other properties such as convergent domestication (Paterson et al., 1995). To investigate such genome wide non-random patterns of association, we also conducted a comprehensive study on the correspondence of QTLs identified here, with previous reports irrespective of the population type. Information on a total of 830 QTLs related to fiber elongation in 103 previous reports were downloaded from CottonGen (https://www.cottongen.org). The hypergeometric probability distribution was used to analyze QTL correspondence, providing a means to infer statistically whether QTLs for a trait are randomly distributed between two populations or environments. Correspondence was identified for 15 of the 31 previously reported studies for the QTLs reported in Acala Maxxa background and for 10 of the 31 previously reported studies in Pima S6 background ( Table 4 ). Common QTLs were identified in most of the 31 reports listed, however, only 15 (or 10) significant P values based on the hypergeometric distribution suggest that across the genome, this correspondence is not sufficient to infer a non-random distribution of QTLs between published studies.

While it is a common practice to compare and contrast QTLs mapped in different populations to verify underlying loci, it is also important to highlight the advantages of some mapping populations over others, especially in terms of the quality of the estimates that they can offer. With reduced background noise, advanced backcross populations such as NILs offer more precise estimates of the effects of underlying QTLs. While the total variations attributable to the QTLs obtained from our NIL populations ( Tables 1 , 2 , Figure 2 ) are not as pronounced as those from F2, F2:3, RILs or other early generation populations reported elsewhere (which is expected given the nature of our populations), the effects of these QTLs individually were significantly larger than those reported by early generation studies of population types that carry a lot of segregating genomic regions ( Table 5 ).

Since the parents used, type of crosses and population sizes, precision of reported location of the QTLs, and type of diagnostic markers for reporting QTLs in prior studies vary greatly, the distribution of the additive allele effects of the QTLs were compared to ours using a two-sample t-test. Absolute values of the effects were used to perform the t-test to remove/reduce cancellations of effect sizes due to the direction of the QTL source (sign convention). The QTLs reported in the current study have significantly larger average additive effects that those from other population types ( Table 5 ). While, in general, fiber elongation QTLs from other population types contributed 0.009 to 0.37 percent additive effects with maximum observed values of 0.6 percent, the average additive effect of QTLs reported here average 1.07 percent with a maximum effect of 2.19 percent. Indeed, the minimum effect of 0.37 percent found in the NILs is greater than most averages or even maximums reported in other studies. These results suggest that NILs offer a much greater advantage in precisely identifying QTL effects, especially relatively subtle effects that might be masked by either genetic or non-genetic variation in other population structures.

Availability of cotton reference genomes (Paterson et al., 2012; Zhang et al., 2015) enabled us to scrutinize physical regions surrounding the identified QTLs for genes/gene families known or suspected to affect fiber quality parameters in cotton. This investigation was limited to tightly linked regions i.e., 1 Mb on both sides of the SNP marker that is most significantly associated with the QTLs. On chromosome 5, the nearest gene of interest to an ELO QTL (EL05.1) was a C2HC zinc finger superfamily protein, a member of a gene family whose gene expression level at 11 days post-anthesis has been reported to be highly correlated with UHM and which has been implicated in secondary cell wall thickening (Al-Ghazi et al., 2009). Similarly, on chromosome 10, the nearest gene of interest to another fiber elongation QTL (EL10.2) was a glucosyltransferase like UDP-Gp gene. The enzymatic activity of this gene increases during the period of development when cotton fiber is synthesizing massive amount of cellulose (Taliercio and Kloth, 2004; Xu et al., 2022), which like protein is involved in cellulose biosynthesis during secondary cell wall formation stage. While it is premature to suggest the active roles of these genes towards fiber trait phenotypes, these genes are potential candidates for these roles given that they are so close to the most significant marker for fiber quality QTLs that have been repeatedly identified in multiple environments and studies. Improved genomic resources together with these novel populations are expected to accelerate candidate gene identification and validation for many such loci.