The wilt-resistant variety Zhongzhimian 2 and the wilt-susceptible variety Junmian 1 were selected for expression analysis under wilt disease stress. The materials were provided by the Cotton Innovation Team of the College of Agriculture at Xinjiang Agricultural University and were planted in the Cotton Cultivation Laboratory of the Key Laboratory of Agriculture and Biotechnology at Xinjiang Agricultural University. Seeds of uniform fullness were selected and sown in a mixture of perlite and sterilized soil (ratio of 1:2), with three seeds per pot. The laboratory was maintained at a daytime temperature of 25.0 to 28.0°C and illuminated with incandescent lighting. Germination occurred under a cycle of 16 hours of light and 8 hours of darkness, with a relative humidity of 80%. When the cotton plants developed three true leaves, the bottom-tear root method was used to inoculate the plants with the wilt pathogen V991, with three replicates for each treatment. Root tissue samples were collected at 0, 1, 3, 6, 12, 24 and 48 h post-wilt stress. All the samples were quickly frozen in liquid nitrogen and stored at -80°C.

The genomic data of G. hirsutum, G. barbadense, G. arboreum and G. raimondii were downloaded from the CottonGen database (https://www.cottongen.org/). The hidden Markov model of the F-box domain (PF00646) and HMMER 3.2.1 software (http://hmmer.org/) were used to identify the F-box protein sequences in G. hirsutum, G. barbadense, G. arboreum and G. raimondii (Prakash et al., 2017). We verified whether the candidate sequences possessed complete F-box and LRR domains using the CDD database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi), ultimately obtaining the FBL protein sequences for G. hirsutum, G. barbadense, G. arboreum and G. raimondii. The physicochemical properties of G. hirsutum FBL proteins were predicted via the online tool ExPASy (https://web.expasy.org/protparam/).

With the default parameters of ClustalW in MEGA X software, multiple sequence alignments of FBL protein sequences from G. hirsutum, G. barbadense, G. arboreum and G. raimondii were performed (Kumar et al., 2018). The phylogenetic tree was constructed using the neighbor-joining (NJ) method with a bootstrap value of 1000. The online tool EvolView (https://www.evolgenius.info/evolview/) was employed to enhance the visualization of the phylogenetic tree. The MEME online tool (http://alternate.meme-suite.org/tools/meme) was used to analyze the motifs of the G. hirsutum FBL gene family, with the number of functional domains set to 10. The xml file generated by MEME, the nwk file of the phylogenetic tree, and the GFF3 file of gene structures were visualized using TBtools software (Chen et al., 2023).

The 2000 bp sequence upstream of the start codon of the FBL gene from G. hirsutum was extracted and submitted to the PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) online tool to predict the cis-regulatory elements in the gene promoter, and TBtools software was used for visualization (Chen et al., 2023). The transcriptome sequencing data of G. hirsutum were downloaded from the NCBI database (accession numbers: PRJNA248163 and PRJNA532694). Fastp (version 0.23.4) software was used for filtering and quality control of the raw data (Chen et al., 2018). HISAT2 software was used to align the filtered data with the G. hirsutum TM-1 reference genome (https://www.cottongen.org/species/Gossypium_hirsutum/ZJU-AD1_v2.1), and StringTie software was used for expression quantification (Kim et al., 2019). Gene expression levels were determined via the FPKM method, and expression heatmaps were created via TBtools software.

Total RNA extraction was performed using an RNA extraction kit (Tiangen, China). cDNA was synthesized from RNA using the M-MLV RTase cDNA Synthesis Kit (TaKaRa, Japan). The qRT–PCR primers used were designed using Primer Premier 5.0 software ( Supplementary Table S1 ), with the GhUBQ7 gene from G. hirsutum used as the internal control. The qRT–PCRs were conducted using the ChamQ Universal SYBR qPCR Master Mix (Novogene Bioinformatics Technology Co., Ltd., Nanjing). The reaction program consisted of 30 seconds at 95°C for initial denaturation, followed by 40 cycles of 10 seconds at 95°C for denaturation, 30 seconds at 60°C for annealing, and 20 seconds at 72°C for extension. Each sample was examined in triplicate. The quantitative data were analyzed using the 2^–ΔΔCT method (Livak and Schmittgen, 2001).

VIGS experiments involving tobacco rattle virus (TRV) were conducted on cotton (Bachan and Dinesh-Kumar, 2012). The vectors pTRV2 and pTRV1 containing specific regions of Gh_FBL43 were subsequently transformed into Agrobacterium GV3101. The restriction sites selected were BamHI and KpnI, and the primer sequences for amplifying the Gh_FBL43 fragment are shown in Supplementary Table S1 . A mixture of Agrobacterium cells containing pTRV1, pTRV2, or the derived plasmids was transformed into 7-day-old seedlings. After inoculation, the seedlings were washed with deionized water to remove excess Agrobacterium residue and were grown at 25°C under a 16 hour/8 hour light/dark cycle. After 2 weeks of cultivation, the plants were inoculated with the yellow wilt pathogen strain V991. At least 30 seedlings were treated in each trial, and the experiments were repeated three times to evaluate plant disease severity at 21 days post-inoculation.

Additionally, RNA-seq with three biological replicates was conducted on samples from the control group with empty vectors and the Gh_FBL43-silenced plants at 0 h and 24 h post-inoculation. The samples were transported to Aide Sen Biotechnology Co., Ltd., in Xinjiang (Urumqi, China) for sequencing using dry ice. RNA was extracted via a polysaccharide and polyphenol plant total RNA extraction kit, and the extraction process was carried out according to the instructions. After RNA extraction, a Nanodrop was used to determine the RNA purity (OD260/280) and concentration and to determine whether the nucleic acid absorption peak was normal. The integrity of the RNA was accurately detected with an Agilent 2100 instrument. The detection indicators included the following: the RIN value, 28S/18S ratio, whether the chromatogram baseline was increased, and the 5S peak. After the concentration and purity of the extracted RNA samples were determined, a library was constructed using an initial total amount of 1 μg for each sample. Preliminary quantification was performed using a Qubit 3.0 fluorometer, followed by detection of the inserted fragments in the library using the Qsep400 high-throughput analysis system. The effective concentration of the library was accurately quantified using the Q–PCR method to ensure library quality (effective library concentration > 2 nM). The libraries were sequenced on an Illumina HiSeq 2500 (California, USA) platform to generate 150 bp paired-end reads. Fastp (version 0.23.4) software was used to remove reads containing adapters and low-quality reads, resulting in clean data (Chen et al., 2018). The clean reads were aligned to the G. hirsutum TM-1 reference genome (https://www.cottongen.org/species/Gossypium_hirsutum/ZJU-AD1_v2.1) using HISAT2 software, and the aligned reads were quantified with StringTie (Kim et al., 2019). DEGs were screened on the basis of gene count values in various samples using DEseq2 software, applying |log2fold change| ≥ 1 and FDR < 0.01 as selection criteria (Love et al., 2014). The fragments per kilobase of exon model per million mapped fragments (FPKM) for each gene was calculated based on its length, and the expression amount mapped to that gene was calculated. R language PCAtools was used to decompose the expression data (FPKM) of all genes into n principal components to describe the characteristics of the original dataset. For all DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted via the hypergeometric test method with the R language ClusterProfiler (version 4.14.4) package (Wu et al., 2021).

To study the variation in the copy number of the FBL gene family during the evolution of cotton, a total of 171 FBL family genes were identified in the genomes of G. hirsutum, G. barbadense, G. arboreum and G. raimondii. In the two diploid cotton species (G. arboreum and G. raimondii), the number of FBL family genes is close to half of that in G. hirsutum. On the basis of the positions of 57 sequences in the chromosomes of G. hirsutum, the FBL genes in this species were named sequentially as Gh_FBL1 to Gh_FBL57 ( Supplementary Table S2 ). The 57 G. hirsutum FBL genes are distributed across 24 chromosomes, with two FBL genes (Gh_FBL28 and Gh_FBL34) located on scaffold chromosomes ( Supplementary Figure S1 ). The A subgenome contains 28 FBL genes, whereas the D subgenome contains 29 FBL genes. Except for there being one more FBL gene on chromosome D07 than on chromosome A07, the numbers of A subgenome and D subgenome genes in the remaining chromosomes are completely consistent. Overall, there remains a strong correspondence between the A subgenome and the D subgenome, which is consistent with the evolutionary relationships of cotton.

To better understand the evolutionary relationships of the cotton FBL gene family, a phylogenetic tree was constructed based on 171 FBL protein sequences from G. hirsutum, G. barbadense, G. arboreum and G. raimondii. Evolutionary tree analysis indicated that the cotton FBL gene family can be divided into five subgroups ( Figure 1 ). Group 1 contains 50 FBL genes, Group 2 contains 17 FBL genes, Group 3 contains 35 FBL genes, Group 4 contains 14 FBL genes, and Group 5 contains 55 FBL genes. The number of FBL genes in each subgroup in both G. hirsutum and G. barbadense was approximately twice that in G. arboreum and G. raimondii. This finding is consistent with the earlier analysis results and aligns with the evolutionary relationships among the cotton species, suggesting that the FBL genes are relatively conserved during the evolution of cotton. Although Group 4 has fewer members, they have been conserved throughout the evolution of cotton, indicating that they may play crucial roles in biological processes.

Analysis of the phylogenetic tree, gene structure, and motifs was conducted on the basis of the full-length sequences, coding sequences (CDSs), and protein sequences of the G. hirsutum FBL genes ( Figure 2 ). Most genes in Group 1 contain three exons and consist of two motifs (motif 3 and motif 10). In Group 2, genes contain either one or three exons and include seven motifs (motif 1, motif 3, motif 4, motif 7, motif 8, motif 9, and motif 10). Group 3 genes possess four or seventeen exons and comprise five motifs (motif 1, motif 2, motif 3, motif 7, and motif 8). Genes in Group 4 have either one or five exons and consist of four motifs (motif 1, motif 3, motif 4, and motif 8). Most genes in Group 5 contain two exons and encompass eight motifs (motif 1, motif 2, motif 3, motif 4, motif 5, motif 6, motif 7, and motif 9).

To investigate the potential functions of the FBL gene in G. hirsutum, the promoter region 2000 bp upstream of the GhFBL gene was analyzed, revealing a varying number of stress- and hormone-responsive cis-acting elements ( Figure 3 ). A total of three types of cis-acting elements related to the stress response were identified: the MYB-binding site involved in drought inducibility (MBS), the cis-acting element involved in defense and stress responsiveness (TC-rich repeats), and the cis-acting element involved in low-temperature responsiveness (LTR). Additionally, eight hormone-related cis-acting elements were identified, including the cis-acting element involved in abscisic acid responsiveness (ABRE), the cis-acting regulatory elements involved in methyl jasmonate (MeJA) responsiveness (the CGTCA motif and TGACG motif), the cis-acting element involved in salicylic acid responsiveness (the TCA element), the auxin-responsive element (the TGA element), and the gibberellin-responsive elements (the P-box, GARE motif, and TATC box). Among these, the elements related to MeJA and abscisic acid were the most abundant, with 124 and 123 such elements, respectively. Each GhFBL promoter contained a different number and type of cis-acting elements, suggesting that these elements may participate in cotton growth and development and in response to environmental stress through different signaling pathways.

Gene expression patterns are typically related to the functions of genes. The expression patterns of cotton FBL genes in different tissues and at various time points under Verticillium wilt stress were examined via RNA-seq analysis ( Figure 4 ). The analysis of tissue-specific expression patterns revealed that 78.95% of the cotton FBL genes presented tissue-specific expression, with the highest expression levels observed primarily in the calyx, petal, stamen, and torus ( Figure 4a ). Following Verticillium wilt occurrence, more than 64.91% of the cotton FBL genes presented changes in expression levels, indicating that FBL family genes may play a significant role in the response of cotton to Verticillium wilt stress ( Figure 4b ). Among them, the five genes (Gh_FBL20, Gh_FBL24, Gh_FBL36, Gh_FBL43 and Gh_FBL48) that presented the most pronounced changes were selected as candidate genes for further validation of cotton resistance to Verticillium wilt.

Through RNA-seq, five candidate genes (Gh_FBL20, Gh_FBL24, Gh_FBL36, Gh_FBL43 and Gh_FBL48) involved in the resistance of cotton to Verticillium wilt were ultimately screened. The transcription levels of these five genes were measured under Verticillium wilt stress conditions via qRT–PCR in resistant (Zhongzhimian 2) and susceptible (Junmian 1) extreme varieties ( Figure 5 ). Four genes (Gh_FBL24, Gh_FBL36, Gh_FBL43 and Gh_FBL48) presented significantly increased expression across different time points in both extreme varieties under Verticillium wilt stress, whereas Gh_FBL20 presented a significant decrease in expression in expression post-infection. Notably, Gh_FBL43 expression peaked (fold change > 9) at 24 hours post-infection, and this gene presented significant differences in expression levels between the two extreme varieties at 3 hours post-infection. Gh_FBL43 is suggested as a candidate gene for further functional analysis of cotton resistance to Verticillium wilt.

To validate the role of Gh_FBL43 in cotton resistance to Verticillium wilt, we silenced Gh_FBL43 in G. hirsutum using VIGS. Cotyledons of seven-day-old cotton seedlings were injected with Agrobacterium cells carrying pTRV1 or pTRV2-GhCLA. Fifteen days later, the new true leaves presented a distinct chlorosis phenotype, indicating that the VIGS system used in this study successfully silenced the gene ( Figure 6a ). The results of the qRT–PCR analysis of the Gh_FBL43-silenced plants compared with the control plants revealed that the expression level of Gh_FBL43 in the silenced plants was significantly lower than that in the control plants treated with an empty vector ( Figure 6b ). After confirming the silencing of the Gh_FBL43 gene, we inoculated the silenced and control plants with the V991 strain of the Verticillium wilt fungus using a root injury-based method. Twenty-one days post-inoculation, the leaves of the silenced plants presented significant yellowing and necrosis, whereas those of the control plants remained in good growth conditions with no noticeable loss of greenness ( Figure 6c ). Further phenotypic investigation revealed that the disease severity and disease index of the Gh_FBL43-silenced plants were significantly greater than those of the control plants ( Figures 6d, e ).

To clarify the pathway by which Gh_FBL43 regulates cotton resistance to V. dahliae, RNA-seq was conducted on both Gh_FBL43-silenced and control cotton plants at 0 hours and 24 hours post-inoculation with the pathogen, with three biological replicates for each sample. A total of 81.52 Gb of clean data was obtained after the raw data were filtered, with each sample yielding more than 6.05 Gb of clean data. The percentage of Q20 bases was above 98.23%, and the percentage of Q30 bases exceeded 94.58%, with the GC content exceeding 44.14% ( Supplementary Table S3 ). The Pearson correlation coefficients among the three replicates of the same sample were greater than 0.97 ( Figure 7a ). Principal component analysis (PCA) revealed that samples of the same biological replicate clustered together, indicating high reliability and reproducibility of the transcriptomic data ( Figure 7b ). Through differential expression analysis, a total of 6,255 DEGs were identified between the Gh_FBL43-silenced and control plants at 0 hours post-inoculation, of which 3,726 were upregulated and 2,529 were downregulated ( Figure 7c ). At 24 hours post-inoculation, 6,760 DEGs were identified between the two groups, consisting of 4,746 upregulated and 2,014 downregulated genes ( Figure 7d ). Overall, 10,928 DEGs were identified between Gh_FBL43-silenced and control plants pre- and post-inoculation, with 4,204 unique DEGs at 0 h and 4,709 unique DEGs at 24 h, including 2,015 shared DEGs ( Figure 7e ).

To further elucidate the functions and pathways associated with the 10,928 DEGs, a Gene Ontology (GO) enrichment analysis was first conducted on these DEGs. Notable enrichment was observed in biological process terms such as the jasmonic acid (JA)-mediated signaling pathway, response to light intensity, cellular response to fatty acid, induced systemic resistance, response to hydrogen peroxide, salicylic acid metabolic process, photosynthesis, leaf senescence, phenol-containing compound metabolic process, and plant organ senescence ( Figure 8a ). Additionally, a KEGG enrichment analysis highlighted significant enrichment in the following pathways: photosynthesis - antenna proteins, alpha-linolenic acid metabolism, flavonoid biosynthesis, phenylalanine metabolism, photosynthesis, beta-alanine metabolism, fatty acid degradation, tyrosine metabolism, glycerolipid metabolism, carotenoid biosynthesis and thiamine metabolism ( Figure 8b ).

Using k-means clustering, the 10,928 DEGs were analyzed, and the pathways involved in each cluster were annotated. Four clusters with statistical significance were identified ( Figure 9a ). Cluster 1 presented the highest expression level at 24 hours post-inoculation in the control plants, encompassing 4,355 DEGs that were significantly enriched in the flavonoid biosynthesis, phenylalanine metabolism, and alpha-linolenic acid metabolism pathways ( Figures 9a, b ). Cluster 2 presented the highest expression level in the control plants prior to inoculation, with expression levels decreasing in both the control and silenced plants post-inoculation, with 3,242 DEGs that were significantly enriched in the photosynthesis, porphyrin metabolism, and cutin, suberin and wax biosynthesis pathways. Cluster 3 had the highest expression level at 24 h post-inoculation in the silenced plants, with a slight increase in expression levels following inoculation. Cluster 3 included 1,486 DEGs that were significantly enriched in the glycerophospholipid metabolism, circadian entrainment and beta-alanine metabolism pathways. Cluster 4 presented minimal changes in expression levels before and after inoculation in the control plants, however, the expression levels in the Gh_FBL43-silenced plants decreased post-inoculation. Cluster 4 included 1,881 DEGs that were significantly enriched in the circadian entrainment, carotenoid biosynthesis and fatty acid degradation pathways.

The enrichment analysis of all the DEGs revealed that jasmonic acid (JA) metabolism, alpha-linolenic acid metabolism and flavonoid-related pathways were significantly enriched. The alpha-linolenic acid metabolism pathway is an important route for the synthesis of JA in plants. Therefore, the expression patterns of JA synthesis genes were first analyzed. The results indicated that the expression levels of most JA synthesis genes decreased after infection in the Gh_FBL43-silenced plants, whereas the control plants presented the highest expression levels 24 hours post infection ( Figure 10a ). Among the rate-limiting enzyme allene oxide synthase (AOS)-encoding genes involved in JA synthesis, all the genes except for GH_A06G0128 and GH_D06G0110 presented peak expression in the control plants at 24 hours post infection. Additionally, a gene encoding another rate-limiting enzyme, allene oxide cyclase (AOC), also presented peak expression in the control plants at this time point. With respect to the genes involved in the flavonoid synthesis pathway, all the genes except for GH_D05G2350 (flavanone 3-hydroxylase, F3H) presented higher expression levels in the control plants than in the Gh_FBL43-silenced plants, with expression levels increasing with infection and peaking at 24 hours post infection in the control plants ( Figure 10b ). These results suggest that Gh_FBL43 may regulate the resistance of G. hirsutum to wilt disease by modulating the expression of genes involved in JA and flavonoid biosynthesis.