QTL and candidate gene identification of the node of the first fruiting branch (NFFB) by QTL-seq in upland cotton (Gossypium hirsutum L.)

Background The node of the first fruiting branch (NFFB) is an important precocious trait in cotton. Many studies have been conducted on the localization of quantitative trait loci (QTLs) and genes related to fiber quality and yield, but there has been little attention to traits related to early maturity, especially the NFFB, in cotton. Results To identify the QTL associated with the NFFB in cotton, a BC4F2 population comprising 278 individual plants was constructed. The parents and two DNA bulks for high and low NFFB were whole genome sequenced, and 243.8 Gb of clean nucleotide data were generated. A total of 449,302 polymorphic SNPs and 135,353 Indels between two bulks were identified for QTL-seq. Seventeen QTLs were detected and localized on 11 chromosomes in the cotton genome, among which two QTLs (qNFFB-Dt2–1 and qNFFB-Dt3–3) were located in hotspots. Two candidate genes (GhAPL and GhHDA5) related to the NFFB were identified using quantitative real-time PCR (qRT-PCR) and virus-induced gene silencing (VIGS) experiments in this study. Both genes exhibited higher expression levels in the early-maturing cotton material RIL182 during flower bud differentiation, and the silencing of GhAPL and GhHDA5 delayed the flowering time and increased the NFFB compared to those of VA plants in cotton. Conclusions Our study preliminarily found that GhAPL and GhHDA5 are related to the early maturity in cotton. The findings provide a basis for the further functional verification of candidate genes related to the NFFB and contribute to the study of early maturity in cotton. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-08164-2.


Background
Cotton is one of the most economically important crops, and it plays an important role in China's agricultural economy. Upland cotton (Gossypium hirsutum L. AADD, 2n = 52) is one of 50 Gossypium species and the main natural fiber crop. It is widely planted and accounts for more than 95% of global cotton production [1]. China has a large population, and the available land for planting crops is decreasing. Improving the multiple crop index and land utilization rate is crucial for maintaining or increasing agricultural production. Early-maturing cotton has a short growth period and generally shows a dwarf and compact plant architecture [2]. Early maturity is an important target trait of cotton breeding, especially with the current adjustments to the cotton planting Open Access *Correspondence: ysx195311@163.com; w.wanghantao@163.com 1 State Key Laboratory of Cotton Biology, Institute of Cotton Research of Chinese Academy of Agricultural Sciences, Anyang 455000, Henan, China Full list of author information is available at the end of the article structure and the promotion of the mechanization of cotton production. Early maturity is a very complex agronomic trait that involves budding date, growth period, flowering time, NFFB, and the height of the NFFB [3]. These quantitative traits related to early maturity are regulated by both quantitative trait loci and the environment, which are reflected in different genetic models [4,5]. In order to improve the land competition between grain and cotton, cotton breeding for early maturity has attracted increasing attention.
The NFFB, an early-maturity trait, refers to the number of node that generates the first fruiting branch (cotyledon node not included) [3]. This trait is associated with the flowering time and provides an estimate of the relative photoperiodism [6]. The NFFB is an easily measured and dependable morphological indicator used to evaluate cotton early maturity [7,8]. A lower NFFB of a cotton variety indicates earlier maturation. The NFFB has direct and indirect effects on the yield percentage before frost through the bud stage [9]. Because of its high heritability and correlation with the final harvested yield, the NFFB has been used as a criterion for evaluating cotton early maturity [7].
Over the past 20 years, a large number of different quantitative traits related to fiber quality, yield, drought tolerance and disease resistance have been reported. However, inadequate attention has been focused on early maturity traits, especially the NFFB of cotton [10][11][12][13][14][15]. Previous studies have identified QTLs related to early maturity traits on most of the 26 cotton chromosomes [16][17][18][19][20][21]. More than 70 QTLs associated with the NFFB of cotton have been detected through linkage mapping or association analysis [2,6,[22][23][24]. These QTLs need further confirmation and may be useful for improving early maturity via molecular marker-assisted selection (MAS). There are significant negative correlations between early maturity and both fiber quality and yield [4]. It is difficult to improve maturity, yield and fiber quality simultaneously with traditional breeding methods. Therefore, it is essential to identify QTLs and genes contributing to the NFFB, and the use of markers closely linked to QTLs for marker-assisted breeding (MAB) is key to achieve the simultaneous improvement of early maturity and other traits.
NFFB is an important agricultural trait of early maturity in cotton, and it is controlled by multiple genes, with separate minor effects. The construction of genetic linkage maps using molecular markers is a reliable method to detect QTLs related to target traits. However, the construction of linkage maps is difficult and labor intensive, especially for large population. Bulked segregant analysis (BSA) is a rapid method to detect DNA markers that are closely linked to candidate genes for a specific phenotype [25,26]. QTL-seq has been used to identify a key QTL for grain length and weight in rice using a near isogenic F 2 population [27]. The genetic rules for Cf-10 in tomato were investigated by combining SNP-index and Indelindex analyses and local QTL mapping using the Kompetitive Allele-Specific PCR (KASP) genotyping approach [28]. In this study, we constructed a high-generation backcross BC 4 F 2 population to create a highly homozygous background genotype, and we selected extreme plants based on the NFFB of the offspring to identify related QTLs. Two DNA bulks of offspring (20 individuals each) showing extremely high or extremely low NFFB phenotypic values were constructed. QTLs related to the NFFB were identified using whole genome resequencing of the two DNA bulks, and two candidate genes associated with the NFFB were predicted and verified using qRT-PCR and VIGS. These results lay a foundation for analyses of the genetic mechanisms underlying cotton earliness and the breeding of improved early-maturing cotton cultivars in the future.

Identification of extreme phenotypic and construction of two DNA bulks
The average NFFB of the early-maturing cotton material RIL182 and the late-maturing cotton material G2005 were 4.5 and 7.65, respectively. Significant phenotypic differences between these two materials were observed (Fig. 1b). There were also significant differences in plant type between these two parents: RIL182 was compact, whereas G2005 was loose (Fig. 1c). The NFFB differed between the two materials (Fig. 1d). The distribution of NFFB among the BC 4 F 2 population in the field during the budding stage is shown in Fig. 1a. Based on the results of the phenotypic investigation, two DNA bulks with average NFFBs of 4.55 and 9.1 were constructed for QTL-seq via the selection of extreme individuals from the BC 4 F 2 population (Fig. 1b).

Quality evaluation of Illumina sequencing data
The results of Illumina sequencing and related data analyses are shown in Table 1. High-throughput sequencing of the NFFB-H-bulk and NFFB-L-bulk generated a total of 120.77 Gb of clean data, with 61.47 Gb from the NFFB-H-bulk (Clean Q30 Bases Rate, 92.58%) and 59.30 Gb from the NFFB-L-bulk (Clean Q30 Bases Rate, 92.00%). For the donor parent RIL182, the total amount of clean data was 61.50 Gb (Clean Q30 bases rate, 92.61%). For the recurrent parent G2005, the total amount of clean data was 61.53 Gb (Clean Q30 bases rate, 92.75%) ( Table 1). After trimming and filtering, over 98% of the clean reads were selected and aligned to the reference genome. The contrast ratios of all genome samples were between 98.85 and 99.48%, and the average sequencing depths for the four samples (excluding N area) ranged from 24.09 to 25.61. The the proportion of regions detected at least once in the reference genome was greater than 82.99%. The Illumina sequencing data for each sample were sufficient, and the sequencing quality was acceptable. Most of the library fragments were 350 bp (Additional file 1: Fig.  S1), which was consistent with our expectations. These results suggested that the sequencing results met the criteria for mutation detection and analysis.

Detection of mutation sites
When the clean reads were aligned to the TM-1 reference genome [29], two types of polymorphic markers, SNPs and Indels, were identified between RIL182 and G2005 and between NFFB-H-bulk and NFFB-L-bulk. Fig. 1 The NFFBs of the BC 4 F 2 population, two parents, and two bulks. a The frequency distribution of NFFBs in the BC 4 F 2 population; b Significance testing of the differences in NFFB in the two parents and two bulks. **, differences at p < 0.01; c-d The plant type and NFFBs of RIL182 and G2005 Table 1 The statistics of quality and depth for sequencing results (a)  After filtering and screening according to quality value, depth, repeatability, neighboring distance, and homozygosity of detected mutation sites, a total of 2,942,211 SNP and 391,699 Indel markers were selected, which covered all 26 chromosomes. The SNP-index of the two bulks was calculated and analyzed. A total of 449,302 polymorphic SNPs and 135,353 Indels were obtained for QTL-seq by further filtering (Additional file 2: Table S1), and the chromosomal distribution of these markers was illustrated in Fig. 2. The highest number of markers were identified on chromosome A11 (135,154 SNPs and 25,249 Indels), and the lowest numbers were detected on chromosome A03 (1689 SNPs and 3781 Indels).

QTL-seq analysis and colocalization with previous studies
We set the window size as 2 Mb and the step size as 100 K, and the mean SNP-index and mean ∆SNPindex in each window were calculated ( Fig. 3a-c). The confidence intervals of the ∆SNP-index were obtained using confidence levels of 0.1, 0.05 and 0.01 (Fig. 3c). We selected the confidence interval with a confidence level of 0.01 and the locations outside the interval were considered to be candidate regions of QTL-seq. Fifty-three regions were obtained covering 21 chromosomes, i.e., all chromosomes except A04, A05, A09, A11 and A13 (Additional file 3: Table S2). To obtain reliable QTLs, we collected most of the QTLs related to the NFFB detected in previous studies [2,6,[22][23][24]. More than 70 QTLs associated with the NFFB of cotton were reported previously using linkage mapping and association analysis. The specific physical positions of these QTLs were obtained by a BLAST search against the reference genome, and were compared with the 53 QTLs detected in the present study (Additional file 3: Table S2). Seventeen QTLs overlapped with the previous results and were located on chromosomes A06, A10, A12, D01, D02, D03, D06, D09, D10, D11 and D12 (Additional file 4: Fig. S2). Only one QTL was identified on each of chromosomes A06, A10, D01, D02, D06, D10, D11, and D12. More than 2 QTLs were detected on each of chromosomes A12 and D03, and a total of 9 QTLs were identified on chromosome D03 using different localization methods and materials.

Identification of candidate genes by expression pattern analysis
A hot spot was defined as four or more QTLs of the same trait within a 20-cM region [30]. The colocalization results showed that two QTLs (qNFFB-Dt2-1 and qNFFB-Dt3-3) were located in hotspots that contained more than 4 QTLs related to the NFFB (Fig. 4a). A total of 1238 genes were located in these two candidate QTLs. To select the genes that are predominantly expressed during flower bud differentiation, the expression profiles of genes located in major QTLs in 12 different tissues were analyzed. A total of 142 genes showed higher expression at the shoot apical meristem than in other tissues; their functional annotations are shown in  Table S3. Twelve of these genes were selected because their homologous genes in Arabidopsis were related to plant reproductive development (Fig. 4b); their FPKM values in twelve tissues are shown in Additional file 6: Table S4. To examine the expression differences of the 12 selected genes between early-maturing and late-maturing cotton materials during flower bud differentiation, qRT-PCR was used to evaluate the expression differences between the two parents (RIL182 and G2005) at four different bud growth stages (two-leaf to five-leaf stages). During these four stages, the relative expression of two genes (Gh_D02G0291 and Gh_ D03G0906) in the early-maturing material RIL182 was significantly higher than that in the late-maturing material G2005 (Fig. 4c-d). The expression profiles of residual genes are shown in Additional file 7: Fig. S3. In addition, between RIL182 and G2005, a SNP mutation (G/A) was detected at − 989 bp upstream of the start codon of Gh_ D02G0291, and an Indel deletion (TA/T) was detected at − 274 bp upstream of the start codon of Gh_D03G0906 (Additional file 8: Fig. S4). These mutations may lead to differences in gene expression between the two parents. Therefore, these two genes (Gh_D02G0291 and Gh_ D03G0906) were considered potential candidate genes for further study.

Gene structure and sequence analysis of GhAPL and GhHDA5
The genomic DNA and full-length cDNA sequences of Gh_D02G0291 and Gh_D03G0906 were extracted from CottonFGD and named GhAPL and GhHDA5, respectively. The GhAPL gene contained 6 exons and 5 introns (Fig. 5a). The open reading frame of GhAPL was 915 bp in length and encoded 304 amino acid residues. The molecular mass of the GhAPL protein was 33.27 kDa, and the isoelectric point was 6.15. Comparison of the protein sequences of GhAPL with its related proteins demonstrated that GhAPL was 68% homologous to AtAPL and the GhAPL protein included a conserved MYBlike DNA-binding domain and a MYB-CC_LHEQLE motif, which indicated that GhAPL was a MYB-CC type transcription factor (Fig. 5b). The GhHDA5 gene contained 13 exons and 12 introns (Fig. 6a). The open reading frame of GhHDA5 was 1830 bp in length and encoded 609 amino acid residues. The molecular mass of the GhHDA5 protein was 68.42 kDa, and the isoelectric point was 5.32. Multiple sequence alignment revealed that GhHDA5 was 64% homologous to AtHDA5 and the GhHDA5 protein included one Hist_deacetyl domain comprised of approximately 300 amino acids with a Zn binding site that coordinated two aspartates and one histidine, indicating that GhHDA5 belonged to the RPD3type histone deacetylases (HDACs) (Fig. 6b).

Silencing of GhAPL and GhHDA5 in cotton
A VIGS assay was performed to further verify the roles that GhAPL and GhHDA5 play in cotton early maturity. The pCLCrVA::GhPDS-silenced plants had an obvious leaves whitening phenotype, which suggested that the VIGS assay was successful (Additional file 9: Fig. S5). qRT-PCR was performed to evaluate the effects of gene silencing, and the results showed that the expression levels of GhAPL and GhHDA5 in positive plants were downregulated by approximately 41 and 36%, respectively, compared to the levels in control (pCLCrVA) plants (Fig. 7b). When flowering was observed in VA plants, the positive plants of VIGS::GhAPL and VIGS::GhHDA5 were not flowering. The positive plants flowered later than the VA plants (Fig. 7a). We also investigated the flowering time (FT), NFFB, and plant height (PH) of VA and positive plants. The VIGS::GhAPL and VIGS::GhHDA5silenced plants had later flowering and higher NFFB than VA plants (Fig. 7c-d). There was no significant difference in plant height (Fig. 7e). These results indicated that the silencing of GhAPL and GhHDA5 delayed the flowering time and increased the NFFB in cotton.

Discussion
Previous studies conducted with different populations and methods have detected hundreds of QTLs and several candidate genes involved in the early maturity of cotton. Among these QTLs associated with early maturity, more than 70 QTLs related to the NFFB of cotton have been identified via linkage mapping or association analysis. Traditional QTL mapping detected 23 QTLs related to the NFFB using F 2 or F 2:3 populations [6,22,23]. With the rapid development of high-throughput sequencing technology, SNPs have been applied for studies of cotton early-maturity traits. In combination with restriction site-associated DNA sequencing (RAD-seq) or genotyping by sequencing (GBS) technology, 52 QTLs related to the NFFB were identified by constructing a genetic map using SNP markers [2,24]. Although some QTLs related to the NFFB have been identified using traditional QTL mapping, it remains difficult to perform subsequent fine mapping to elucidate the regulatory network of the NFFB when the QTLs are located in large physical intervals. Traditional QTL mapping is generally performed by genotyping a large number of individuals in segregating populations, which is labor intensive, time consuming, and costly [31]. The major advantages of QTL-seq are that the qualified SNPs filtered between the parental lines directly serve as DNA markers, and the use of SNP-index allows accurate evaluation of the frequencies of parental alleles and the genomic contributions of the two parents to the offspring [32]. The present study used QTL-seq to identify QTLs and candidate genes related to the NFFB in cotton using a BC 4 F 2 population. Two major QTLs for the NFFB were detected by QTL-seq and colocalization, and two candidate genes related to early maturity were obtained using qRT-PCR and VIGS assays. Through the comparison of the physical positions of QTLs related to the NFFB in previous studies, seventeen QTLs were obtained, including two key QTLs (qNFFB-Dt2-1 and qNFFB-Dt3-3) that were located in hotspots related to the NFFB and were detected four and five times, respectively (Fig. 4). Among the nine QTLs, four stable QTLs (qNFFB-D2-3, qNFFB-D2-4, qNFFB-D3-1 and qNFFB-D3-2) were detected in at least two years, explaining 6.33-16.33% of the phenotypic variation (PV) with LOD scores from 3.34-9.89 [24], and one major QTL (qNFFB-D03-1) explained 32.57% of the PV with an LOD score of 17.12 [2]. Notably, among the detected QTLs related to the NFFB, 9 QTLs were mapped using different localization methods and materials on chromosome D03 (Additional file 4: Fig. S2), which contained the most candidate regions; these results are in agreement with previous reports that chromosome D03 contains potential segments and multiple sites for early-maturity traits [19,33]. These two key QTLs were found simultaneously by different populations and methods, which supports the reliability of the study results and lays a foundation for further fine mapping and candidate gene discovery.
Several candidate genes related to early maturity in cotton have been identified and verified in recent years. Five early maturity-related genes, Gh_D03G0922, Gh_ D03G0718, Gh_D03G0728, Ghir_D03G011310, and Ghir_A05G017290, were identified using genome-wide association study (GWAS) and qRT-PCR, and their molecular functions were further confirmed by VIGS or genetic transformation of Arabidopsis [19,[33][34][35]. Additionally, an important regulatory gene associated with flower bud differentiation, GhCAL, was found by analyzing the transcript dynamics of shoot apex samples collected from two early-maturing and two latematuring varieties, and its function was verified in transgenic plants [36]. The relative expression of all the genes detected in previous studies by qRT-PCR exhibited significant differences between early-maturing cotton varieties and late-maturing cotton cultivars during flower bud differentiation. Flower bud differentiation is a physiological and morphological symbol of the transition from vegetative growth to reproductive growth, and it has been reported that the stage of flower bud differentiation occurs earlier in early-maturing cotton than in late-maturing cotton [36,37]. In early-maturing cotton, the flower bud differentiation generally begins at two-leaf stage or three-leaf stage [36,38]. Therefore, qRT-PCR is a reliable method to identify candidate genes related to the early maturity of cotton during flower bud differentiation. Among the candidate genes we identified, two genes (GhAPL and GhHDA5) exhibited marked differences in expression levels between the two parents from two-leaf satge to five-leaf stage. Therefore, these genes may be associated with earlymaturity traits. GhAPL is a homolog of AtAPL, which encodes an MYB-related protein in Arabidopsis. AtAPL, a key regulator of phloem development, promotes flowering through the transcriptional activation of FLOW-ERING LOCUS T (FT) and its transport machinery component, FTIP1 [39]. Sequence analysis has revealed that the GhAPL protein contains a conserved SH[A/L] QKY[R/F] motif within the MYB-like domain, similar to the SHAQKYF-class of MYB-like proteins, and previous studies have shown that the SHAQKYF-class of MYB-like proteins in Arabidopsis play essential roles in regulating the flowering time [39,40]. GhHDA5, a member of RPD3-type HDACs, is homologous to AtHDA5, which regulates flowering time by suppressing the expression of FLOWERING LOCUS C (FLC) and MAF1 in Arabidopsis [41]. In addition, other Arabidopsis RPD3-type proteins have been reported to regulate the flowering time by different regulatory pathways [42][43][44][45]. The GhHDA5-silenced plants in our study showed obviously later flowering and higher NFFB compared to VA plants, which is consistent with the finding in a previous study that RNAi-suppressed GhHDA5 lines showed delayed flowering [46]. These two genes were preliminarily found to be related to the early maturity of cotton in this study. However, further functional verification is necessary to clarify their specific functions.
Early maturity is a complicated physiological and biochemical process involving a large number of structural, regulatory and biochemical pathway-related genes. By identifying and cloning genes related to the NFFB, the molecular improvements in cotton maturity are performed. It is of great significance for the realization of mechanized planting, the promotion of double cropping of the cotton and wheat, and the sustained and stable development of the cotton industry.

Conclusions
In this study, we used QTL-seq and colocalization to identify two key QTLs related to the NFFB on chromosomes D02 and D03 of the cotton genome. Then, two candidate genes (GhAPL and GhHDA5) related to early maturity were identified by qRT-PCR and VIGS assays. Both genes exhibited the higher expression levels in the early-maturing cotton material RIL182 than in the late-maturing material G2005 during the flower bud differentiation period, and the silencing of GhAPL and GhHDA5 delayed the flowering time and increased the NFFB compared to VA plants in cotton. These results provide a basis for the further functional verification of candidate genes related to the NFFB and contribute to the study of early maturity in cotton.

Plant materials
A population comprising 137 recombinant inbred lines (RILs) was developed in a previous study from a cross between the early-maturing variety CCRI36 and the latematuring material G2005 [24]. The NFFBs of CCRI36 and G2005 are approximately 5 and 8, respectively. Based on the phenotype and genetic composition of 137 recombinant inbred lines, RIL182 with an average NFFB of 4.5 was selected as the basic material, and the parent G2005 was used as the recurrent parent to generate BC 1 F 1 hybrids. To increase the homozygosity of the background genotype, hybrid BC 1 F 1 plants with lower NFFB were consecutively backcrossed to G2005 three times to create a BC 4 F 1 population. A BC 4 F 2 population including 278 individual plants was obtained by self-crossing for the mapping and identification of target QTLs related to the NFFB of cotton.

Identification of extreme phenotypes and construction of two DNA bulks
The BC 4 F 2 population, RIL182, and G2005 were planted in Anyang, Henan Province. The NFFB of cotton is the number of main stem nodes from the cotyledon node to the first fruiting branch node. We investigated the NFFB of the BC 4 F 2 population in the field during the budding stage. According to the results of phenotypic investigation, plants with more than 8 NFFBs were considered high-NFFB plants, and plants with fewer than 6 NFFBs were considered low-NFFB plants. Two DNA bulks for QTL-seq were constructed via the selection of extreme individuals from the BC 4 F 2 population. Twenty individuals with an average NFFB of 9.1 were sampled as the high-NFFB bulk (NFFB-H-bulk), and 20 individuals with an average NFFB of 4.55 were sampled as the low-NFFB bulk (NFFB-L-bulk).

Extraction of genomic DNA and Illumina sequencing
Genomic DNA of the female parent RIL182, male parent G2005, and extreme individuals were extracted from the young leaf tissues using the cetyltrimethyl ammonium bromide (CTAB) method [47]. RNA enzymes were used to purify the DNA pools. The quality and concentration of DNA were examined by electrophoresis on 1% agarose gels and an ultraviolet spectrophotometer, respectively. Tested and qualified DNA samples were sent to ANOROAD company (Beijing, China). Samples with sufficient purity, concentration and volume were used to build libraries using an Illumina TruSeq DNA Sample Preparation Kit (Illumina Inc., San Diego, CA, USA). The libraries were constructed by randomly interrupting the DNA, recovering the required length of DNA fragments using electrophoresis, and adding the adapter primer. Qualified libraries were sequenced using an Illumina HiSeq 2500 sequencing platform on a computer.

Filtering of the original data and sequencing quality control
Raw data, i.e., the original data of the Illumina highthroughput sequencing results, were stored in FASTQ files containing the name of each sequence, the base sequence and sequencing quality information. To ensure the quality of the information analysis, clean reads were generated by filtering the original data. Three steps were performed to filter the original data: (1) adapter-polluted reads (reads containing more than five bases polluted by the adapter) were removed; (2) low-quality reads (reads in which greater than 15% of the bases had a mass value lower than 19) were removed; and (3) reads with an N ratio greater than 5% were removed. After filtering the raw data, the proportion of clean reads was higher, indicating better sequencing quality. Sequencing quality was also evaluated by the GC content distribution and base mass distribution.

Detection of mutation sites and QTL-seq analysis
Clean short reads were mapped against the TM-1 reference genome [29] using Burrows-Wheeler Aligner (BWA) software with the BWA-MEM algorithm [48]. The mutation analysis software GATK [49] was used to detect the SNPs and Indels of the population, and further filtering was performed according to some factors, such as quality value, depth and repeatability. We detected reliable mutation sites, which were annotated accordingly using ANNOVAR software [50] and an existing genome annotation file (gff/ gtf). The SNP-index represents the proportion of mutant bases in the offspring, and we used QTL-seq approach to calculate the SNP-index of the two bulks [32]. To reduce the impacts of sequencing and comparison errors, the following criteria were used for SNP filtering: (1) parental genotypes were heterozygous; (2) the parents and progeny did not cover the corresponding site or the depth was less than 10; (3) the SNP-index was less than 0.3 or more than 0.7 in both bulks; (4) the SNPs were not located on main chromosomes. SNPs that met any of the criteria were removed. After filtering, the chromosomal distribution of the qualified SNPs and Indels was visualized using CMplot package based on R, and their ΔSNP-index values between two bulks were calculated by the following formula.
A ∆SNP-index value of 1 would suggest that the locus of the bulk DNA was entirely from G2005, a ∆SNP-index value of − 1 would indicate that the locus of the bulk DNA was completely from RIL182, and a ∆SNP-index value of 0 would indicate that the two DNA bulks had the same SNP-index at this site. Therefore, an absolute value of the ∆SNP-index closer to 1 may indicate a candidate locus for the NFFB in cotton. According to the preset sliding window (the size of the window was 2 Mb, and the step size was 100 kb), the mean of the SNP-index in each window was calculated for the qualified SNPs. The confidence intervals of the ∆SNP-index were obtained under confidence levels of 0.1, 0.05 and 0.01. We selected a genetic interval with a confidence level of 0.01 and the regions outside of the interval were considered the candidate regions of QTL localization [32].

Colocalization and candidate gene identification
To obtain reliable regions related to the NFFB, we compared our results with all of the QTLs identified in previous studies. The physical positions of markers at both ends of previously identified QTLs were found in the CottonFGD database (http:// www. cotto nfgd. org/) [51]. Depending on the specific physical location on the chromosome, regions that overlapped with other QTLs more than four times were considered major QTLs. To analyze the expression patterns of the genes located in major QTLs, the FPKM values of 12 tissues, i.e., root, stem, fiber, ovule, leaf, calycle, torus, petal, pistil, stamen, and buds of two-leaf and three-leaf stages, were obtained from publicly available transcriptomic data [29,36]. The genes showing relatively higher expression during flower bud differentiation were selected, and their homologous genes in Arabidopsis were identified. Their functions were searched in TAIR (https://www.arabidopsis.org/) and STRING (https:// string-db. org/) database to identify candidate genes associated with early maturity traits [52,53].

qRT-PCR
Total RNA was extracted from terminal buds of the RIL182 and G2005 materials at the two-leaf to five-leaf stages using a Plant RNA Purification Kit (Tiangen, China). cDNA was obtained via reverse transcription using the PrimeScript ™ RT Reagent Kit (Perfect Real Time) (Takara, Japan). The composite cDNA samples were diluted 5-fold and used as cDNA templates for the qRT-PCR experiment. The qRT-PCR experiment was performed using a 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) and UltraSYBR Mixture (Low ROX) (Cwbio, China) with three biological repeats and three technical replicates. The specific gene primer pairs used for PCR amplification are listed in Additional file 10: Table S5. G. hirsutum Actin (GhActin) was used as an endogenous control [54], and the 2 −ΔΔCt method was used to calculate the gene expression levels [55].

Gene structure and protein sequence analysis of GhAPL and GhHDA5
The CDS, protein, and genomic DNA sequences of GhAPL and GhHDA5 were extracted from the Cot-tonFGD database [51]. The intron-exon structures of the GhAPL and GhHDA5 genes were analyzed using the online software GSDS2.0 (http:// gsds. cbi. pku. edu. cn/) [56]. The conserved domains within the GhAPL and GhHDA5 proteins were predicted using the online software SMART (http:// smart. embl-heide lberg. de/) [57]. The Mw and pI were predicted using ExPASy (http:// web. expasy. org/ compu te_ pi/) [58]. DNAMAN 6.0 software was employed to perform multiple sequence alignments of the GhAPL and GhHDA5 proteins with their homologous proteins from other species.

Virus-induced gene silencing
The CDS of GhAPL and GhHDA5 were submitted to SGN-VIGS Tool (http:// vigs. solge nomics. net/) for sequence analyses and selection of two 300-bp optimal sequences, which were amplified from the RIL182 cDNA to construct the pCLCrVA::GhAPL and pCLCrVA::GhHDA5 vectors. The pCLCrVA::GhAPL and pCLCrVA::GhHDA5 plasmids were transformed into LBA4404 strains. The pCLCrVA vector served as an empty control, and the pCLCrVA::GhPDS vector served as a positive control. The LBA4404 strains carrying pCLCrVA::GhAPL, pCLCrVA::GhHDA5, pCLCrVA (negative control), or pCLCrVA::GhPDS (positive control) were mixed with the strain containing pCLCrVB (helper vector) (1:1 ratio, OD600 = 1.5) and infiltrated into two fully expanded cotyledons of 10-day-old CCRI36 seedlings via syringe infiltration to generate the VA plants, pCLCrVA::GhAPL-and pCLCrVA::GhHDA5-silenced cotton, which were cultivated in a climate-controlled greenhouse with a suitable growing environment (light/ dark cycle: 16 h at 28 °C/8 h at 22 °C). When the leaves of the pCLCrVA::GhPDS plants became white, the virusinduced silencing plants were confirmed using PCR and specific primers (Additional file 11: Table S6). The plants were transplanted into large plots and cultivated until flowering. The process of VIGS assay was the same as described previously [59].