Gossypium arboreum line Shixiya 1 originated in South China and was subsequently introduced to the Yangtze and Yellow River regions. Shixiya 1 seeds were used for EMS treatment. The CU mutants were obtained from EMS mutagenesis, which possess a stable inheritable leaf curling trait. Gossypium arboreum lines DQJ and BML were used as male parent lines to construct the F₂ population. All of the plants were sown in the experimental field of the Institute of Cotton Research, Chinese Academic Agricultural Sciences (ICR, CAAS, Anyang, Henan Province, or Sanya, Hainan Province).

In order to obtain abundant mutants with novel traits, approximately 15,000 seeds of Shixiya 1 were presoaked in phosphate buffer (100 mM, pH 7.0) for 12 h at 28°C. Then, the seeds were soaked in EMS (0.6%, ethyl methanesulfonate) phosphate buffer for 8 h in the dark at 28°C with continuous flipping and shaking during the period. Subsequently, the seeds were rinsed three times with distilled water for 30 min to remove residual EMS. Finally, they were sown in the experimental field of ICR, CAAS (Sanya). The M₂ seeds were harvested by collecting one boll per plant from the surviving 7987 M₁ plants (53.3% survival). After four generations of self-fertilization, 123 individual lines with visible phenotypes were identified out of 5,980 M₇ plants, and one of which exhibited curly leaf and was named as CU.

The CU and DQJ were crossed to construct an F₂ population (60 F₁ lines) in 2020, and an F₂ population (1,100 plants) was generated by self-fertilization of the F₁ population in 2021.

Fresh young leaves were collected from parental lines (DQJ and CU), 50 F₂ individual plants exhibiting curly leaf and 50 F₂ individual plants exhibiting flat leaf. Total DNA was extracted from the leaves following the previously reported method (Ali et al., 2019). Two DNA pools were constructed. One pool consisted of DNAs from 50 F₂ plants showing the dominant phenotype (curly leaf), and another pool consisted of DNAs from 50 F₂ plants showing the recessive phenotype (flat leaf). The experimental process is performed according to the standard protocol provided by Illumina for bulk segregation analysis (BSA)-sequencing. The DNA of each sample was randomly broken into 350-bp fragments by ultrasonic fragmentation to construct the sequencing library. The library was sequenced by Illumina HiSeq after passing the quality inspection. Each sample was sequenced at 30× coverage of the assembled genome with 150-bp paired-end reads. Clean reads were mapped to the Shixiya 1 genome ( https://www.cottongen.org/species/Gossypium_arboreum/CRI-A2_genome_v1.0 ) by BWA after filtering the raw reads. According to the positioning results of clean reads in the reference genome, Picard was used to mark duplicates and GATK for local realignment, base recalibration, and other preprocessing to ensure the accuracy of the detected SNPs. The single-nucleotide polymorphism was further detected by GATK and filtered to obtain the final SNP site set. BSA-sequencing analysis was performed using the runQTLseqAnalysis method of QTLseqr for SNP ratio calculation. A 1,000-kb window spans multiple gene-regulatory units—an optimal scale for quantitative trait locus (QTL) scans—and the 10-kb step size yields 99% inter-window overlap to ensure continuous signal capture. Accordingly, we calculated the average ΔSNP index based on a 1,000-kb sliding window with a 10-kb step size to identify regions associated with the target trait.

InDel markers are specific PCR primers designed based on nucleotide insertions or deletions (InDels) identified between parental genomes. First, the genomic positions of InDels are extracted from the Variant Call Format (VCF) file generated through BSA-seq using a Perl script. Subsequently, a local Perl script retrieves the flanking sequences surrounding each InDel locus. Finally, primers are designed using Primer 3.0 software. For InDel marker polymorphism screening, PCR products were generated using the 3G Taq Master Mix for PAGE (Red Dye) kit and subsequently resolved through 8% non-denaturing polyacrylamide gel electrophoresis.

Based on the BSA-seq results, a significant single signal was identified on Chr03, indicating that the curly leaf trait is controlled by a single locus. Combined with the Indel marker mapping results, KASP (Kompetitive Allele Specific PCR) markers were designed within the candidate region (Chr03: 4.31 Mb to 4.52 Mb) by HuaZhi Biotechnology Co., LTD. The development of KASP markers followed the workflow established in previous studies (Liu et al., 2023).

The leaves in the different stages T1 (the first true leaf at the pre-curling stage), T2 (the first true leaf in the process of curling) and T3 (after the curling degree of the first true leaf stabilized) were taken from the wild type Shixiya 1 and CU mutant, respectively. Three biological repeats were used for total RNA extraction with Tiangen RNAprep Pure Plant Plus Kit (DP441). Two cDNA libraries were constructed following Illumina standard protocol and sequenced on a HiSeq 2000 platform. All paired-end reads were aligned to the G. arboreum reference genome through the Bowtie software (Li et al., 2009; Langmead and Salzberg, 2012). Differentially expressed genes (DEGs) were analyzed by Cufflinks software for FPKM (fragments per kilo bases per million reads) calculation (Trapnell et al., 2010).The differentially expressed genes (DEGs) between Shixiya 1 and CU mutant were identified using the criteria of |log2^{Fold Change}| ≥ 1 and P-value ≤0.05. The KAAS (https://www.genome.jp/tools/kaas/) was used for pathway annotation analysis of the DEGs. The RT-qPCR assays were performed by the TB Green^® Fast qPCR Mix (Takara, RR430B). The sequences of each pair of primers used for the different RT-qPCR assays are listed in the Supplementary Table S1 , and His3 was used as the endogenous reference gene for the relative quantitation of the gene expression data (Gong et al., 2017). The RT-qPCR analysis was conducted with three biological replicates, followed by the application of the 2^−ΔΔCt method to quantify the relative gene expression level (Pfaffl, 2001; Livak and Schmittgen, 2001).

To observe the differences in cellular structure within the leaves between Shixiya 1 and CU, the third to the last leaf was taken from the wild-type Shixiya 1 and the mutant line CU, respectively. The samples were fixed in 4% (w/v) paraformaldehyde, 0.25% (w/v) glutaraldehyde, and 50 mM sodium phosphate (pH 7.2), dehydrated by passage first through an ethanol series and then through a xylene ethanol series, embedded in paraffin, and finally sectioned to 20 μm by using a Leica RM2235 rotary microtome. The paraffin sections were stained with safranin stain solution and observed under an Olympus BX53 microscope.

For the virus-induced gene silencing (VIGS) assay, we used the TRV vector (i.e., pTRV1 and pTRV2), along with VIGS positive control pTRV2–PDS and negative control (NC) empty vector pTRV2. A 294 bp GaIAA14 fragment was cloned from the leaf of Shixiya 1 cDNA with the pTRV2-GaIAA14 F/pTRV2-GaIAA14 R primer pair ( Supplementary Table S1 ). The pTRV1 helper plasmid, pTRV2, pTRV2-PDS, and pTRV2-GaIAA14 were all introduced into the Agrobacterium strain GV3101. VIGS assay and silencing efficiency check were performed with CU mutant as donor plant according to the previously published protocol (Gong et al., 2017).

The curly leaf mutant was obtained from wild-type Shixiya 1 by EMS mutagenesis, and its trait was stably inherited after multiple generations of self-purification. Here, it was designated as the CU (curly leaf) mutant. Compared with the wild-type Shixiya 1, the plant height of the CU mutant was slightly reduced ( Figures 1a, b ). Notably, the leaves of the CU mutant were cup-curled, and the leaves’ surface was bumpy and darker ( Figure 1c ). Furthermore, we sliced the third to the last leaf of CU mutant and Shixiya 1 to analyze the cellular structure difference under the microscope during development. Compared with Shixiya 1, the CU mutant exhibited more developed and tightly packed spongy tissue, along with thinner palisade tissue ( Figures 1d, e ). The abnormal development of spongy tissue may be one of the reasons for the formation of cup-curled and bumpy leaves in the CU mutant.

To identify the causal locus for the curly leaf of CU mutant, a G. arboreum line DQJ (flat leaf) and the CU (curly leaf) were used as parents to construct a F2 population through hybridization and selfing. All F₁ plants exhibited curly leaf, but less curly than those of CU mutants ( Supplementary Figure S1 ), indicating that the curly leaf allele is semi-dominant. The F₂ population produced a 1:2:1 segregation ratio of curly leaf, moderately curly leaf, and flat leaf phenotypes ( Supplementary Table S2 ), indicating that a single locus controlled the curly leaf of CU mutant. Furthermore, we identified the causal locus by bulk segregation analysis (BSA), and the results showed that the ΔSNP index greater than 95% confidence intervals was located on chromosome A03 from 3.30 to 16.10 Mb ( Figures 2a,b ) by aligning the sequence reads to the reference Shixiya 1 genome ( https://www.cottongen.org/species/Gossypium_arboreum/CRI-A2_genome_v1.0 ), demonstrating that the locus controlling the curly leaf is located on chromosome A03.

We further designed the molecular markers in the 12.8-Mb region on the basis of InDels between the two BML (BaiMiLa) and CU parents. A total of 70 markers were designed, and 24 markers showed significant polymorphisms with PCR between the two parents that were used for fine mapping in a BML × CU F₂ population ( Supplementary Table S3 ). Finally, the curly leaf locus (CU) was narrowed down to a 199.10-kb region flanked by Chr03_4317636 and Chr03_4516774 ( Figure 3a ). Based on the above-cited candidate interval, we first developed 13 KASP markers covering the genomic region from Chr03_4345697-KASP to Chr03_5016706-KASP ( Supplementary Table S4 ). Using these markers to perform genotype testing in 211 recessive lines in the DQJ × CU F₂ generation population, we found that the curly leaf trait was closely linked to the Chr03_4468542-KASP~Chr03_4516774 interval located at the 48.2-kb genomic region ( Figure 3b ).

There are six genes in the region from 4468542 to 4516774 on chromosome A03 of G. arboreum ( Supplementary Table S5 ) (Cotton Functional Genomics Database: https://cottonfgd.org/ ). We aligned the sequencing data of the six candidate genes between the two parents. The results showed that Ga03G0409 had an SNP in the second exon ( Supplementary Figure S2c ), Ga03G0407, Ga03G0410, and Ga03G0411 had SNPs at the 2000 bp region upstream of the start codon ( Supplementary Figures S2a, d, e ), and there was no difference in the nucleotide sequence of other genes between the two parents ( Supplementary Figures S2b, e ). Furthermore, the RNA-Seq analysis indicated that Ga03G0407, Ga03G0410, and Ga03G0411 did not show a clear alteration in transcription between Shixiya 1 and CU mutant. Therefore, we identified Ga03G0409 as the potential candidate gene of CU and then cloned and sequenced the putative ORF from RNA. Compared with Shixiya 1, a single-nucleotide mutation from G to A at the 301st position of the second exon occurred in the CU mutant ( Figure 3c ; Supplementary Figure S3a ), which resulted in a valine to isoleucine (V101I) mutation in the Ga03G0409 protein ( Supplementary Figure S3b ). The gene was annotated as IAA14 after the sequence alignment analysis, in which a single amino acid substitution occurs at the conserved degron motif GWPPV(I) within its domain II, resulting in the change from GWPPV to GWPPI. Further structural analysis of the protein revealed that the IAA14 protein in the CU mutant lacked a segment of α-helical structure ( Figure 3d ). This structural alteration likely impaired the normal function of IAA14, ultimately leading to leaf curling in cotton plants.

To further test the function of GaIAA14, a 294-bp fragment of GaIAA14 was cloned and inserted into pTRV2 for VIGS. At 2 weeks post-infiltration of the CU mutants, the TRV: GaIAA14 displayed flat leaves compared to that maintained curly leaves with a blank TRV: 00 vector ( Figure 4a ). The RT-qPCR analysis revealed that the expression of GaIAA14 was much lower in GaIAA14-silenced plants compared with that in blank vector control TRV: 00 plants ( Figure 4b ), indicating that GaIAA14 is the key gene for the curly leaf phenotype associated with the CU locus.

GaIAA14, a member of the AUX/IAA gene family, is an early auxin response factor (Zhao et al, 2021). Theoretically, the mutations of GaIAA14 would affect the auxin signaling pathway. Therefore, RNA-Seq was performed with the first true leaf of Shixiya 1 and CU at different stages: T1 (the first true leaf at the pre-curling stage), T2 (the first true leaf in the process of curling), and T3 (after the curling degree of the first true leaf stabilized). In total, 1,789 gene transcripts were identified. Notably, the stage-specific DEGs (1,684) in the T₃ stage were significantly more than that in the T₁ (4) and T₂ (71) stages ( Figure 5a ), indicating that the stage-specific genes involved in regulating leaf curling in the CU mutant were predominantly active during the T3 stage. A total of 1,711 differentially expressed genes (DEGs) were identified in the T₃ stage, including 760 upregulated genes and 951 downregulated genes ( Figure 5b ; Supplementary Table S6 ). To further understand the molecular functions of these DEGs in the process of cotton leaf morphogenesis, the KEGG enrichment analysis was conducted at the T₃ stage. The results showed that both upregulated and downregulated genes were associated with plant hormone signal transduction pathways, including 40 and 42 DEGs, respectively ( Figures 5c, d ; Supplementary Table S7 ). In plant hormone signal transduction pathways, the auxin signaling pathway contained the highest number of genes, totaling 15 genes (six upregulated DEGs and nine downregulated DEGs). These included six AUX/IAA family members, four auxin response factors (ARFs), three GH3 auxin-responsive promoters, and two auxin transporter-like proteins ( Supplementary Table S7 ). Collectively, these results imply that GaIAA14 regulates cotton leaf morphology, which is potentially dependent on auxin distribution and signaling.