The cotton material used in this study was the G. hirsutum variety Xinlu Zhong 61, which was provided by the Economic Crop Research Institute of the Xinjiang Academy of Agricultural Sciences. Seeds of appropriate maturity, plumpness, and uniformity were selected and sown into pots filled with a mixture of perlite and sterilized soil (1:2 ratio), with four seeds per pot. The daytime temperature in the laboratory was maintained at 25–28°C, with incandescent lighting provided on a 16-h light/8-h dark cycle for germination. After germination, seedlings with consistent growth were selected, and the substrate attached to the roots was gently removed before seedlings were transferred to a hydroponic box for cultivation. The hydroponic box was a 30 L plastic container, and the culture mixture was Hoagland’s nutrient mixture, which was changed weekly. A foam board, cut to fit the culture box, was placed on top of the box, with uniformly spaced holes (2 cm in diameter, 15 holes per board) for plant insertion. The stems were wrapped in sponge strips and inserted into the holes in the foam board, with the roots submerged in the nutrient mixture. Continuous and uniform aeration was provided by an air pump to ensure proper root respiration. When the seedlings reached the three-leaf stage, they were treated with 100 mM NaHCO₃ to induce alkaline stress. The plants were grown in a growth chamber with a daytime temperature of 28°C, a nighttime temperature of 25°C, and a 16-h light/8-h dark cycle. Samples were collected at 0 h, 2 h, 4 h, 8 h, 12 h, and 24 h of stress treatment, rapidly frozen in liquid nitrogen, and stored at -80°C for subsequent RNA-seq, metabolomic sequencing, and qRT–PCR analyses.

The samples were sent on dry ice to Mavimetabolism (Wuhan, China) for RNA-seq. RNA extraction was performed using TRIzol reagent (Invitrogen). A certain amount of total RNA was extracted and fragmented into smaller pieces. The fragmented mRNA was then mixed with primers, and first-strand cDNA was synthesized using PCR. Next, a second-strand synthesis reaction was performed, and the second-strand products were recovered. The resulting cDNA was subjected to end repair, the addition of an A base, and adapter ligation. PCR amplification of the ligated products was performed, followed by purification and recovery. The final library was tagged, completing cDNA library construction. Library quality was assessed using an Agilent 2100 instrument and Q–PCR. The constructed library was sequenced using the Illumina HiSeq 2500 platform. After the raw sequencing data were obtained, Fastp software (Chen et al., 2018) was used to remove adapter sequences, filter low-quality reads, and eliminate sequences with a greater than 5% N content, resulting in clean reads suitable for analysis. 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 (Pertea et al., 2016). The alignment results were quantified using featureCounts. Setting the FDR standard to less than 0.01 ensures that fewer than 1% of the DEGs screened are caused by random differences. FDR<0.01 and |log₂fold change|>1 were used as the criteria for screening DEGs (Love et al., 2014). GO and KEGG enrichment analyses were performed on all DEGs using the clusterProfiler (version 4.14.4) package (Wu et al., 2021), with statistical significance determined by hypergeometric tests. The protein sequences of all DEGs were submitted to PlantTFDB (https://planttfdb.gao-lab.org) for analysis and prediction to obtain differentially expressed transcription factors (TFs).

After vacuum freeze-drying, 50 mg of the sample was weighed, mixed with 1000 μL of extraction solution (methanol/acetonitrile/water, 2:2:1 v/v) and vortexed for 30 seconds. Steel beads were added, and the sample was processed with a 45 Hz grinder for 10 minutes, followed by ultrasonication for 10 minutes (in an ice–water bath). The sample was then left to stand at -20°C for 1 h, followed by centrifugation at 12,000 rpm for 15 minutes. The supernatant (500 μL) was carefully collected and transferred to an EP tube. The extract was dried in a vacuum concentrator and then redissolved in 160 μL of extraction solution (acetonitrile/water, 1:1, v/v). The mixture was vortexed for 30 seconds, ultrasonicated in an ice–water bath for 10 minutes, and centrifuged again at 12,000 rpm for 15 minutes. The supernatant was carefully collected for subsequent analysis. Metabolites were analyzed using a Waters Acquity I-Class PLUS ultrahigh-performance liquid chromatography system coupled with an AB Sciex Qtrap 6500+ high-sensitivity mass spectrometer. Chromatographic separation was performed on a Waters Acquity UPLC HSS T3 column (1.8 μm, 2.1*100 mm) with an injection volume of 2 μL. The mass spectrometry conditions were as follows: electrospray ionization (ESI) temperature of 550°C; ion spray voltage (IS) of 5,500 V in positive ion mode and -4,500 V in negative ion mode; and ion source gases I (GSI), II (GSII), and curtain gas (CUR) set to 50 psi, 55 psi, and 35 psi, respectively, with collision-induced dissociation (CID) parameters set to moderate intensity.

Metabolite identification was performed with the in-house database GB-PLANT using secondary mass spectrometry data. Isotope signals, redundant signals from K⁺, Na⁺, and NH₄ ⁺ ions, and fragment ions corresponding to higher-molecular-weight compounds were removed. The metabolites were quantified in multiple reaction monitoring (MRM) mode with a triple quadrupole mass spectrometer. After the mass spectrometry data for different samples were obtained, the peak areas of all the metabolites were corrected. Principal component analysis (PCA) was conducted on the metabolite content data matrix using R software. Metabolite classification and pathway functional annotations were performed using the KEGG database (http://www.genome.ad.jp/kegg/) to identify the major biochemical metabolic and signal transduction pathways involved. Partial least squares regression (PLSR) was applied to establish a model of the relationships between metabolite levels and sample categories for predictive modeling. DAMs were identified using the criteria of a fold change >2 or a fold change < 1/2 and a P value <0.05.

The gene expression profiles of the DEGs were subjected to a coexpression analysis using the dynamic branch cutting method in the R package WGCNA (Langfelder and Horvath, 2008). The weighting coefficient β was chosen to yield a correlation coefficient of approximately 0.8 with a certain level of gene connectivity to ensure a scale-free network. In this study, β=8 was selected as the weighting coefficient. The network was constructed via an automatic network construction function, blockwise modules, resulting in multiple valid modules, each containing a different number of genes. Modules with a similarity greater than 0.75 were merged using minModuleSize = 30 and Merge Cut Height = 0.25 as the criteria. The module eigengene (ME) was calculated and correlated with hormone levels and different treatment durations. Specific modules were selected using the criteria of r>0.80 and P<0.05. The coexpression network was visualized using Cytoscape (Shannon et al., 2003) (version 3.10.0) software.

Homologs of the candidate genes were identified using the BLASTn (Altschul et al., 1990) function on the CottonGen (Yu et al., 2021) website (https://www.cottongen.org/), and specific primers were designed using DNAMan (Lynnon Corporation, Canada). Total RNA was extracted from samples collected after various durations of stress with the RNAprep Pure Polysaccharide Polyphenol Plant Total RNA Extraction Kit (Tiangen, Beijing, China) according to the manufacturer’s instructions, and first-strand cDNA was synthesized using a reverse transcription kit (abm). Cotton Ubiquitin7 (GhUBQ7) was selected as the internal reference gene, and qRT–PCR amplification of the relevant genes was performed on an Applied Biosystems™ 7500 Fast Real-Time PCR System (three biological replicates). Relative gene expression was analyzed using the 2^−ΔΔCt method (Livak and Schmittgen, 2001), and the data were visualized using GraphPad Prism version 8.0.1 for Windows. All primers used in this study are listed in Supplementary Table S2 .

RNA-seq was performed on 18 samples of G. hirsutum subjected to alkaline stress at six time points (0 h, 2 h, 4 h, 8 h, 12 h, and 24 h). After filtering, we obtained a total of 153.89 Gb of clean data, with each sample yielding more than 6.65 Gb of clean data. The Q30 base percentage exceeded 95.90%, and the GC content was above 43.55% ( Supplementary Table S1 ). The correlation between biological replicates not only indicates the reproducibility of experimental procedures but also reflects the reliability of the DEG identification and helps identify potential outliers. The Pearson correlation coefficients for the three biological replicates of the same sample were greater than 0.96 ( Figure 1a ). Principal component analysis (PCA) revealed that samples from the same biological replicate clustered together, confirming the reliability and reproducibility of the RNA-seq data ( Figure 1b ). PCA revealed a significant separation trend among the samples from the alkaline stress control group and groups exposed to alkaline stress for different times. As the duration of stress increased, the distance between the samples from the treatment groups and the control samples gradually increased, suggesting that the response of G. hirsutum to alkaline stress became more pronounced with prolonged treatment time.

The differential expression analysis revealed a total of 751 DEGs at 2 h of alkaline stress compared with 0 h, including 344 upregulated and 407 downregulated genes, with 148 unique DEGs ( Figures 2a, b ). At 4 h of alkaline stress, 1,367 DEGs were identified, with 652 upregulated and 715 downregulated genes and 624 unique DEGs, compared with those at 0 h. At 8 h of alkaline stress, 1,317 DEGs were identified, with 813 upregulated and 504 downregulated genes and 486 unique DEGs, compared with those at 0 h. At 12 h of alkaline stress, 1,770 DEGs were identified, with 797 upregulated and 973 downregulated genes and 353 unique DEGs, compared with those at 0 h. At 24 h of alkaline stress, 4,598 DEGs were identified, with 1,802 upregulated and 2,796 downregulated genes and 2,840 unique DEGs, compared with those at 0 h. In total, 6,610 DEGs were identified in plants under alkaline stress, with 71 common DEGs detected among the groups. KEGG enrichment analyses were conducted for the DEGs at different time points of alkaline stress to reveal the dynamic changes in the response of G. hirsutum to alkaline stress ( Figure 2c ). The DEGs at 2 h were significantly annotated to the MAPK signaling pathway, phenylpropanoid biosynthesis pathway and flavonoid biosynthesis pathway. The DEGs at 4 h were significantly annotated to glucosinolate biosynthesis, linoleic acid metabolism, glycerophospholipid metabolism and zeatin biosynthesis pathways. The DEGs at 8 h were significantly annotated to the phenylpropanoid biosynthesis, starch and sucrose metabolism, glycerolipid metabolism and fatty acid elongation pathways. The DEGs at 12 h were significantly annotated to the flavonoid biosynthesis, alpha-linolenic acid metabolism, linoleic acid metabolism and phenylpropanoid biosynthesis pathways. The DEGs at 24 h were significantly annotated to the flavonoid biosynthesis, glutathione metabolism, phenylalanine metabolism and alpha-linolenic acid metabolism pathways. Flavonoid biosynthesis and alpha-linolenic acid metabolism pathways were significantly annotated at 12 h and 24 h, during which the greatest number of DEGs was observed, indicating that these time points are the most intense periods of the response of G. hirsutum to alkaline stress and that flavonoid biosynthesis and alpha-linolenic acid metabolism play important roles during these periods. GO enrichment analyses were performed on the 6,610 DEGs. The GO analysis revealed significantly enriched biological processes, including the cellular response to phosphate starvation, the salicylic acid catabolic process, the hormone catabolic process, the flavonoid biosynthetic process, the regulation of the cell cycle process, the organic acid catabolic process, the negative regulation of leaf senescence, the jasmonic acid-mediated signaling pathway, inorganic anion transmembrane transport, and the phenylpropanoid metabolic process ( Figure 2d ).

K-means clustering was applied to the 6,610 DEGs, identifying 6 statistically significant clusters. A KEGG pathway enrichment analysis was performed for each cluster ( Figures 3a, b ). Cluster 1 presented a decrease in expression at 2 h of stress, an increase at 4 h, and a gradual decrease thereafter, with the lowest expression observed at 24 h. This cluster contained 2,657 DEGs and 189 TFs. Significant enrichment was detected in the following pathways: the cell cycle, cysteine and methionine metabolism, and flavonoid biosynthesis. Cluster 2 showed a gradual decrease in expression after 4 h of stress, with the lowest expression observed at 24 h. This cluster contained 504 DEGs and 24 TFs. Significant enrichment was detected in the flavonoid biosynthesis and betalain biosynthesis pathways. Cluster 3 presented an increase in expression at 2 h of stress, followed by a decrease at 4 h and then a gradual increase, with expression peaking at 24 h. This cluster contained 1,017 DEGs and 165 TFs. Significant enrichment was detected in the phototransduction, carotenoid biosynthesis, and circadian entrainment pathways. Cluster 4 presented an increase in expression at 2 h of stress and a slight decrease from 4 h to 8 h, followed by a gradual increase, with the highest expression occurring at 24 h. This cluster contained 1,088 DEGs and 76 TFs. Significant enrichment was detected in the pathways of alpha-linolenic acid metabolism, fatty acid degradation, and glutathione metabolism. Cluster 5 showed a rapid increase in expression under stress, with expression peaking at 4 h, followed by a gradual decrease. This cluster contained 631 DEGs and 77 TFs. Significant enrichment was detected in the zeatin biosynthesis, thiamine metabolism, and phenylalanine metabolism pathways. Cluster 6 presented an increase in expression under stress, with expression peaking at 8 h, followed by a gradual decrease. This cluster contained 713 DEGs and 48 TFs. Significant enrichment was detected in the phenylalanine metabolism, cytochrome P450, and circadian rhythm pathways.

A total of 4,225 metabolites were identified in 18 samples of G. hirsutum subjected to alkaline stress at six time points (0 h, 2 h, 4 h, 8 h, 12 h, and 24 h) via UPLC–MS. PCA revealed that samples of the same biological replicates clustered together, indicating that the metabolomic data were reliable and reproducible ( Supplementary Figure S1 ). The 4,225 identified metabolites were classified into 10 categories. Amino acids and their derivatives accounted for 35.10% of the total, organic acids accounted for 16.77%, and benzene and substituted benzene derivatives accounted for 11.31% ( Figure 4a ). Through a differential abundance analysis, a total of 489 differentially accumulated metabolites (DAMs) were identified at 2 h compared with 0 h, with 206 upregulated and 283 downregulated DAMs, including 135 unique DAMs ( Figures 4b, c ). At 4 h, 469 DAMs were identified, with 192 upregulated and 277 downregulated DAMs, including 105 unique DAMs. At 8 h, 766 DAMs were identified, with 382 upregulated and 384 downregulated DAMs, including 278 unique DAMs. At 12 h, 647 DAMs were identified, with 220 upregulated and 427 downregulated DAMs, including 180 unique DAMs. At 24 h, 698 DAMs were identified, with 387 upregulated and 311 downregulated DAMs, including 200 unique DAMs. A total of 1,684 DAMs were identified, including 57 shared DAMs. The KEGG enrichment analysis of the 1,684 DAMs revealed significant enrichment in pathways related to flavonoid biosynthesis, starch and sucrose metabolism, photosynthesis, alpha-linolenic acid metabolism, anthocyanin biosynthesis, folate biosynthesis, phenylalanine metabolism, phenylpropanoid biosynthesis, carbon fixation via the Calvin cycle, the pentose phosphate pathway, and ascorbate and aldarate metabolism ( Figure 4d ).

Using k-means clustering, 6 statistically significant clusters were identified from 1,684 DAMs, and metabolic classification was performed for each cluster ( Figures 5a, b ). Cluster 1 showed a gradual increase in metabolite levels under stress, with metabolite levels peaking at 12 h, followed by a decrease. This cluster contained 209 DAMs, classified mainly as amino acids and their derivatives and benzene and substituted benzene derivatives. Cluster 2 maintained a relatively stable metabolite level before 12 h, and the metabolite level increased sharply at 24 h. This cluster contained 320 DAMs, classified mainly as amino acids and their derivatives and benzene and substituted benzene derivatives. Cluster 3 showed a decrease in metabolite levels at 2 h and 8 h, with increases at 4 h and 12 h. This cluster contained 238 DAMs, classified mainly as amino acids and their derivatives, and organic acids. Cluster 4 showed a gradual increase in metabolite levels under stress, with the levels peaking at 8 h, followed by a decrease. This cluster contained 344 DAMs, classified mainly as amino acids and their derivatives and benzene and substituted benzene derivatives. Cluster 5 presented a slight increase in metabolite levels after 2 h of stress, followed by a gradual decrease. This cluster contained 252 DAMs, classified mainly as amino acids and their derivatives, benzene and substituted benzene derivatives, and organic acids. Cluster 6 showed a gradual decrease in metabolite levels under stress, with the minimum levels observed at 24 h. This cluster contained 321 DAMs, classified mainly as amino acids and their derivatives and organic acids.

The KEGG enrichment analysis of all the DAMs and DEGs revealed that both groups were significantly enriched in the flavonoid biosynthesis, phenylalanine metabolism, alpha-linolenic acid metabolism, and ascorbate and aldarate metabolism pathways ( Figure 6a ). First, we analyzed the expression patterns of genes in the flavonoid biosynthesis pathway. In addition to the flavonol synthase (FLS) and hydroxycinnamoyl transferase (HCT) genes, other flavonoid biosynthesis genes, such as chalcone isomerase (CHI), chalcone synthase (CHS) and dihydroflavonol 4-reductase (DFR), presented the highest expression at 8 h ( Figure 6b ). Moreover, the levels of flavonoid metabolites, such as delphinidin 3-O-beta-D-sambubioside, malvidin 3-(6’’-p-coumarylglucoside), peonidin 3-(6’’-p-coumarylglucoside), petunidin 3-(6’’-p-coumarylglucoside), pelargonidin 3-sambubioside 5-glucoside, and glycyphyllin, peaked at 24 h ( Figure 6c ). On the other hand, the levels of naringenin chalcone, peonidin-3-O-alpha-arabinopyranoside, petunidin 3-(6’’-acetylglucoside), 5H-benzylohepten-5-one, and 1,8-bis((2R,3R)-3,5,7-trihydroxy-2H-1-benzopyran-2-yl)-3,4,6-trihydroxy, and 2’,3,4,4’,6’-pentahydroxychalcone peaked at 0 h, with their levels decreasing after exposure to stress. We calculated the correlation between genes and metabolites to further explore the relationships between the expression of flavonoid biosynthesis genes and metabolite levels, visualizing those with an absolute correlation coefficient greater than 0.8 and a p value less than 0.05 ( Figure 6d ). A total of 25 flavonoid biosynthesis genes were significantly correlated with 9 flavonoid pathway metabolites, with the expression of 14 genes exhibiting negative correlations with the levels of 6 metabolites and the expression of 16 genes showing positive correlations with the levels of 5 metabolites. Specifically, GH_D11G1872 (CHI) expression was significantly negatively correlated with the peonidin 3-(6’’-p-coumarylglucoside), delphinidin 3-O-beta-D-sambubioside, and malvidin 3-(6’’-p-coumarylglucoside) levels.

The alpha-linolenic acid metabolism pathway is an important route for the synthesis of jasmonic acid (JA) in plants (Choudhary et al., 2025). Both the alpha-linolenic acid metabolism and JA-mediated signaling pathways were significantly enriched in differentially expressed genes (DEGs) and differentially accumulated metabolites (DAMs) at multiple stages. To this end, an analysis was conducted on the changes in the levels of genes and metabolites associated with other JA biosynthesis pathways, with an initial focus on the variations in the expression levels of JA biosynthesis genes ( Figure 7a ). The rate-limiting enzymes of JA synthesis, allene oxide cyclase (AOC, GH_D08G0423) and allene oxide synthase (AOS, GH_A04G1388, GH_A12G0282, and GH_D06G0112), presented the highest expression after 4 h of stress. The JA level peaked at 12 h of stress, whereas the level of methyl dihydrojasmonate peaked at 24 h of stress ( Figure 7b ). The correlations between genes and metabolites were calculated to further explore the relationships between the expression of JA biosynthetic genes and metabolite levels, and those with absolute correlation coefficients greater than 0.8 and p values less than 0.05 were visualized ( Figure 7c ). A total of 14 genes were significantly correlated with 4 metabolites in the JA biosynthesis pathway, with the expression of 7 genes exhibiting significant negative correlations with JA and methyl dihydrojasmonate levels and the expression of 9 genes exhibiting significant positive correlations with methyl dihydrojasmonate, (+)-7-iso-JA, and colfosceril palmitate levels. GH_D10G0585 (LOX2S) expression was significantly positively correlated with methyl dihydrojasmonate and (+)-7-iso-jasmonic acid levels and significantly negatively correlated with JA levels.

TFs are key regulators of gene expression and play important roles in plant growth, development, and stress response mechanisms (Yan et al., 2025; Yang et al., 2025). In this study, we identified a total of 579 differentially expressed TFs (DE-TFs) from the 6,610 DEGs. The major TF families identified included the AP2/ERF, MYB, bHLH, WRKY, NAC, bZIP, C2H2, HD-ZIP, and GRAS families ( Figure 8a ). Using k-means clustering, we identified 5 statistically significant clusters of these DE-TFs, and we further examined the TFs with the highest fold changes in each cluster ( Figure 8b ). Cluster 1 showed a gradual decrease in TF expression under stress, with the lowest expression occurring at 24 h. The TFs with the greatest fold changes in expression in this cluster were GH_D11G1947 (MYB), GH_A13G0333 (MYB) and GH_A11G0465 (MYB). Cluster 2 presented the highest TF expression at 0 h and 4 h. The TFs with the greatest fold changes in expression in this cluster were GH_D11G1947 (MYB), GH_A09G2571 (AP2/ERF) and GH_A08G1918 (AP2/ERF). Cluster 3 showed an increase in TF expression, with the expression peaking at 24 h. The TFs with the greatest fold changes in expression in this cluster were GH_D05G2112 (MADS), GH_D08G0585 (bHLH) and GH_A05G2081 (MADS). Cluster 4 presented the highest TF expression at 4 h. The TFs with the greatest fold changes in expression in this cluster were GH_D08G1271 (bHLH), GH_A08G2170 (AP2/ERF) and GH_A08G1275 (bHLH). Cluster 5 presented the highest TF expression at 8 h. The TFs with the greatest fold changes in expression in this cluster were GH_D09G2504 (AP2/ERF) and GH_D11G0488 (MYB).

WGCNA was performed on the expression matrix of the 6,610 DEGs to identify the core genes related to alkali tolerance in cotton, which resulted in 13 distinct coexpression modules ( Figure 9a ). The correlation between each module and the duration of stress was calculated. The results indicated that the pink module was highly correlated with 0 h, the green module was correlated with 4 h, the brown module was correlated with 8 h, the magenta module was correlated with 12 h, and the blue module was correlated with 24 h ( Figure 9b ). By calculating the kME (eigengene connectivity) value of each module gene, the gene with the highest absolute kME value was used as the hub gene of each module. Three hub genes were identified for each module, resulting in a total of 15 candidate genes associated with salt tolerance in cotton ( Figure 9c ).

The functions of these 15 candidate genes were annotated based on the homologous genes in A. thaliana ( Table 1 ). GH_A03G1732 (NAC), GH_A05G2081 (MADS), GH_A13G2571 (MYB), and GH_D02G1891 (NAC) encode four TFs that are involved primarily in the response to stress. GH_A10G1801 encodes a late embryogenesis abundant (LEA) protein that is crucial for maintaining osmotic pressure and protecting cell membrane structures. GH_A01G0048 encodes a 3-ketoacyl-CoA synthase (KCS) that is involved in the biosynthesis of very-long-chain fatty acids. GH_A01G2352 encodes a glutathione S-transferase (GST) that is involved in scavenging ROS. GH_A05G4235 and GH_D05G2419 encode a chalcone isomerase (CHI) and a chalcone synthase (CHS), respectively, which are involved in flavonoid biosynthesis. GH_A07G2044 encodes the UDP-glycosyltransferase 82A1 (UGT82A1), which catalyzes the glycosylation of flavonoids. GH_A09G0143 encodes a calmodulin-like (CML) protein, a calcium signal effector. GH_A09G0432 encodes a lysine histidine transporter (LHT) that is involved in amino acid absorption and transport via roots. GH_A12G0029 encodes a sulfate transporter (SULTR) that is involved in sulfate transport. GH_A12G0282 encodes an allene oxide synthase (AOS) that is involved in JA biosynthesis. GH_D06G2326 encodes an anthocyanidin-3-O-glucosyltransferase (UFGT) that is involved in the biosynthesis of anthocyanins.

The relationships between these 15 candidate genes and salt tolerance in cotton was further explored using qRT–PCR to analyze their expression patterns in cotton under salt stress at different time points ( Figure 9d ). Among these genes, GH_A01G0048, GH_A03G1732, GH_A05G2081, GH_A05G4235, GH_A07G2044, GH_A09G0143, GH_A09G0432, GH_A10G1801, GH_A12G0282, GH_A13G2571, GH_D05G2419 and GH_D06G2326 presented significant increases in expression under salt stress. Three genes (GH_A01G2352, GH_A12G0029 and GH_D02G1891) presented significant decreases in expression. In summary, 15 candidate genes associated with cotton salt tolerance, including 4 TFs, were identified through WGCNA and qRT–PCR.