Phenotypic performance and correlation among parameters
We obtained data on 14 traits from 11 environments for 200 materials under two different low-temperature stresses, including germination, emergence, and osmoregulation substance-related traits of upland cotton (Table S2). The cold tolerance coefficient was used to perform descriptive statistics and correlation analysis. HL, TL, and SSI under CC stress, and DW and GI under DVC stress were partially normally distributed, whereas the others were normally distributed (Fig. 1). In the CC environment, stress delayed the ER of seeds after 5 d as 0.00–1.00. The indices with higher coefficients of variation were ER, RL, and SAAI. Under DVC environment, GR was reduced (0.02–0.81) with a mean of 0.12. GR and GI had higher coefficients of variation than that of the other parameters.
A certain correlation existed between the corresponding traits, with a positive correlation between seedling emergence traits (ER, HL, RL, and TL) under the two conditions (Fig. 2). Varieties with higher emergence rates also had longer hypocotyls and roots, and higher number of whole seedlings after some days of recovery. Among the varieties with higher ER after the 4°C stress, the SS content increased, and the SAA content decreased. Under DVC stress, GR and DW were positively correlated with ER and total length, respectively.
For the 14 traits obtained in two chilling stress, SS and SAA contents were obtained only in one environment, therefore, the interaction between environment and cultivar could not be calculated. The variances of the remaining 12 traits all showed environmental variation > germplasm × environmental (G × E) interaction variation > germplasm variation, and there were extremely significant differences (Table 1). After CC treatment, ER-CC, RL-CC and TL-CC had lower generalized heritability (H2) ranging from 14.50 to 24.33, indicating that these traits were more affected by constant low temperature. The H2 of ER-DVC, HL-DVC and DW-DVC was smaller (8.49–24.87) after DVC treatment, indicating that DVC treatment had a greater effect on the emergence rate, hypocotyl length and dry matter weight. It showed that the measured traits were mainly affected by the environment, and secondly by G × E, and different traits were affected to different degrees.
Table 1 Multi-factor ANOVA for phenotypic parameters under CC and DVC stress
* p < 0.05, and ** p < 0.01 statistical significance. The same as below
SNP identification and population structure analysis
Based on re-sequencing, a total of 2,606,756 highly consistent SNPs were produced. The statistical results for the distribution of each chromosome are shown in Fig. 3. The average physical distance between SNPs was 2.75 Kb. According to the phylogenetic tree, the cotton accessions were divided into four groups, Groups I, II, III, and IV, which contained 77, 18, 20, and 85 accessions, respectively (Fig. 4a). Based on the first three axes of PCA, Groups II and IV were distinguished from the other two accessions, which was consistent with the results of the phylogenetic tree (Fig. 4b).
Next, the geographical distribution of the four subgroups was analyzed. In Group I, YR, FR, and YZR accounted for the largest proportion of the total number of sources (42.37, 52.38, and 62.50%, respectively). In Group III, all varieties were NW (10.00%), YR (13.56%), FR (4.76%), and YZR (12.50%), and no varieties from NEM. Group II had the lowest number of varieties, including five geographical sources with a relatively small proportion. The accessions from NW and NEM made up a large part of Group IV, and the proportion of each region was 61.25% and 46.67%, respectively (Fig. 4c). The four subgroups showed different phenotypes under CC and DVC stress. Broadly, Group IV had better performance in the seedling stage (ER, HL, RL, TL, and DW) and germination processes (GR and GI), and the cold resistance coefficients were relatively higher under both conditions. Under CC stress, Group I had the lowest emergence rate and dry weight, and HL, RL, and TL showed little difference among Groups I and III. In addition, the SS content showed an increasing trend after 24 h of 4°C treatment, and the varieties in Group III showed a maximum increase, whereas the SAA content showed a decreasing trend. No significant difference was observed among the four subgroups. Under DVC stress, Groups I and III exhibited moderation of the seeding and seedling stage. Group II had the lowest emergence and germination rates (Fig. 4d). Overall, Group IV performed the best under both treatments and had the strongest cold resistance.
Identification of candidate genes by GWAS
Subsequently, a GWAS was performed using 14 quantitative traits. A total of 577 significant association SNPs were identified at -log10(p) ≥ 6.41 (p = 1/n; Bonferroni correction) and FDR < 0.05. These significant associations included 229 and 348 significant associations for CC-and DVC-related traits, respectively (Table S4). The reoccurring and very closely linked SNPs in multiple tolerance indices of all significant markers detected in a range of 200 kb (distances when dropped to half of the maximum value) were seen as common loci with strong signals (Ge et al. 2022). A threshold value (− log10 (p) ≥ 6.41) equivalent to p ≤ 0.0000004 was used to indicate a significant association between SNPs. This threshold value was more stringent than those reported in other low-temperature GWAS studies (Zeng et al. 2021; Thapa et al. 2020). A threshold value of p ≤ 0.0001 (− log10 (p) ≥ 4) was used to identify the most significant SNPs present in more than one environment; however, SNPs with a higher association threshold (− log10 (p) ≥ 6.41) in at least two environments were considered common. Next, SNPs that recurred between thresholds 4–6.41 were screened by four algorithms. A total of 35 quantitative trait nucleotides (QTNs) were significantly associated with cold tolerance indices (Table 2). Of these 35 QTNs, 14 harbored at least two SNPs from various chilling indices. qGhDW-A02-1 and qGhRL-D11-1 were localized using indicators under CC conditions. The QTNs qGhHL-D03-2, qGhHL-D06-1, qGhRL-D05-2, and qGhRL-D10-2were all localized as indicators of cold resistance under DVC conditions. The remaining 29 QTNs were obtained by the joint association of phenotypic indicators under the two low-temperature stresses. It is assumed that the regulatory mechanisms of cotton under the two different low-temperature stresses are mostly consistent; however, we observed a few differences.
Table 2
SNPs (− log10(p) ≥ 4) associated with traits, including QTN name, trait, co-located QTNs, chromosome (Chr), physical location, number of significant SNPs (− log10(p) ≥ 6.41), and gene region
QTN name
|
Trait
|
Co-located QTNs in trait
|
Chr
|
Start (bp)
|
End (bp)
|
Number of significant SNPs
|
Gene region
|
qGhRL-A01-1
|
RL-CC5
|
TL-CC5; RL-DVC5; TL-DVC5
|
A01
|
9510779
|
9766034
|
1
|
Gh_A01G0554-Gh_A01G0560
|
qGhDW-A02-1
|
DW-CC3
|
SAA-CC1
|
A02
|
44385186
|
44478779
|
1
|
Gh_A02G1006
|
qGhRL-A02-2
|
RL-CC4
|
DW-CC4; TL-CC4; RL-DVC5; TL-DVC5
|
A02
|
2896600
|
3078996
|
1
|
Gh_A02G0240-Gh_A02G0264
|
qGhER-A03-1
|
ER-DVC3
|
RL-CC6; TL-CC6; DW-DVC2
|
A03
|
87267372
|
87749071
|
1
|
Gh_A03G1236-Gh_A03G1238
|
qGhGI-A03-1
|
GI-DVC1
|
GR-DVC1; HL-DVC1
|
A03
|
79414427
|
79603749
|
1
|
Gh_A03G1094-Gh_A03G1095
|
qGhRL-A03-3
|
RL-CC2
|
RL-DVC5; DW-DVC5; DW-CC6;
HL-DVC2; TL-CC2; TL-DVC1; RL-DVC1
|
A03
|
96024997
|
96273909
|
2
|
Gh_A03G1526-Gh_A03G1543
|
qGhRL-A04-1
|
RL-CC2
|
DW-DVC5; TL-CC2; HL-DVC4;
HL-DVC1
|
A04
|
17321042
|
17703996
|
2
|
Gh_A04G0452-Gh_A04G0456
|
qGhRL-A04-2
|
RL-CC6
|
TL-CC6; HL-CC4; DW-CC4
|
A04
|
49951029
|
50096820
|
1
|
Gh_A04G0730
|
qGhHL-A05-1
|
HL-DVC2
|
GI-DVC1; TL-DVC2; RL-DVC2;
HL-CC4; DW-CC3
|
A05
|
8053310
|
8228229
|
2
|
Gh_A05G0782-Gh_A05G0809
|
qGhRL-A05-2
|
RL-CC6
|
DW-CC5; DW-DVC5; DW-CC6
|
A05
|
16032367
|
16791751
|
1
|
Gh_A05G1573-Gh_A05G1609
|
qGhRL-A06-1
|
RL-CC6
|
TL-CC2; RL-CC2
|
A06
|
68371212
|
68711833
|
1
|
Gh_A06G1102-Gh_A06G1107
|
qGhHL-A07-1
|
HL-DVC5
|
RL-CC3; HL-DVC4; DW-CC2
|
A07
|
14032784
|
14379401
|
1
|
Gh_A07G0817-Gh_A07G0833
|
qGhDW-A10-1
|
DW-DVC5
|
DW-CC6
|
A10
|
5430537
|
5499846
|
5
|
Gh_A10G0504-Gh_A10G0509
|
qGhGI-A10-2
|
GI-DVC2
|
RL-DVC5; TL-CC5; GR-DVC2; ER-CC5; HL-CC5
|
A10
|
4906273
|
4997697
|
4
|
Gh_A10G0479-Gh_A10G0480
|
qGhER-A12-1
|
ER-DVC4
|
TL-DVC2; HL-CC3; HL-CC2; TL-CC3;
|
A12
|
538629
|
735546
|
1
|
Gh_A12G0025-Gh_A12G0047
|
qGhDW-A13-1
|
DW-DVC5
|
DW-DVC4; DW-CC4; ER-DVC4
|
A13
|
3044291
|
3074621
|
1
|
Gh_A13G0257-Gh_A13G0258
|
qGhHL-A13-2
|
HL-DVC2
|
ER-DVC4; RL-DVC4; TL-DVC3;
TL-DVC2; GI-DVC1; RL-CC6; RL-DVC3
|
A13
|
5298800
|
5555034
|
2
|
Gh_A13G0396-Gh_A13G0408
|
qGhDW-D02-1
|
DW-DVC4
|
DW-DVC3; DW-DVC5;ER-DVC5;
DW-CC6; DW-CC4; RL-CC5
|
D02
|
53946149
|
54310473
|
2
|
Gh_D02G1563-Gh_D02G1575
|
qGhTL-D02-2
|
TL-CC2
|
SS-CC1
|
D02
|
52484587
|
52875858
|
1
|
Gh_D02G1511-Gh_D02G1519
|
qGhDW-D03-1
|
DW-DVC5
|
DW-CC4; HL-CC1; HL-CC5; TL-CC5;
DW-CC4; DW-CC6; GI-DVC1; HL-DVC4
|
D03
|
41956639
|
42386456
|
3
|
Gh_D03G1338-Gh_D03G1370
|
qGhHL-D03-2
|
HL-DVC2
|
TL-DVC2; DW-DVC2; RL-DVC5;
TL-DVC5; GI-DVC1
|
D03
|
36113046
|
36261250
|
1
|
Gh_D03G1082-Gh_D03G1088
|
qGhHL-D05-1
|
HL-DVC2
|
TL-DVC2; RL-DVC2; TL-DVC3; RL-DVC3; RL-DVC5; TL-CC5; HL-CC5;
RL-CC5; TL-CC5
|
D05
|
6399248
|
6809890
|
7
|
Gh_D05G0783-Gh_D05G0823
|
qGhRL-D05-2
|
RL-DVC5
|
ER-CC6; TL-DVC5; RL-DVC4; TL-DVC4
|
D05
|
52787137
|
52806498
|
1
|
Gh_D05G3304
|
qGhHL-D06-1
|
HL-DVC2
|
GI-DVC1; GR-DVC1; TL-DVC2
|
D06
|
9313582
|
9706669
|
6
|
Gh_D06G0578-Gh_D06G0593
|
qGhHL-D05-2
|
HL-DVC2
|
GI-DVC1; TL-DVC2
|
D06
|
9708291
|
11280704
|
18
|
Gh_D06G0594-Gh_D06G0652
|
qGhRL-D05-1
|
RL-DVC5
|
GI-DVC1; ER-CC2; TL-DVC5
|
D07
|
26421009
|
26722023
|
2
|
Gh_D07G1497Gh_D07G1501
|
qGhDW-D09-2
|
DW-CC2
|
TL-CC2; HL-CC2; ER-CC4; ER-CC1;
RL-CC1; RL-CC3; RL-CC4; TL-CC4;
ER-DVC4; HL-DVC5; DW-DVC2
|
D09
|
5962380
|
6453208
|
3
|
Gh_D09G0185-Gh_D09G0212
|
qGhER-D09-2
|
ER-CC4
|
RL-DVC5; ER-DVC4; GI-DVC1;
DW-DVC2
|
D09
|
6813191
|
7122473
|
2
|
Gh_D09G0223-Gh_D09G0233
|
qGhGI-D10-1
|
GI-DVC1
|
DW-DVC5; DW-DVC3; TL-CC3;
SAA-CC1
|
D10
|
23804915
|
24048758
|
1
|
Gh_D10G1289-Gh_D10G1295
|
qGhRL-D10-2
|
RL-DVC2
|
TL-DVC2; GI-DVC1
|
D10
|
24400297
|
24616040
|
1
|
Gh_D10G1311-Gh_D10G1313
|
qGhRL-D11-1
|
RL-CC2
|
RL-CC3; TL-CC2; RL-CC6; RL-CC3
|
D11
|
50773560
|
51019195
|
1
|
Gh_D11G2494-Gh_D11G2499
|
qGhRL-D11-2
|
RL-CC2
|
TL-CC2; DW-DVC3; RL-CC3
|
D11
|
60175737
|
60210242
|
1
|
Gh_D11G2947-Gh_D11G2948
|
qGhRL-D12-1
|
RL-CC1
|
TL-CC5; DW-DVC2; RL-CC4; TL-CC4
|
D12
|
4436632
|
4592816
|
1
|
Gh_D12G0320-Gh_D12G0329
|
qGhRL-D12-2
|
RL-CC6
|
DW-DVC5; DW-CC6
|
D12
|
52426228
|
52640575
|
1
|
Gh_D12G1938-Gh_D12G1951
|
qGhGI-D13-1
|
GI-DVC1
|
HL-DVC4; TL-CC5; TL-CC3
|
D13
|
4999879
|
5182668
|
1
|
Gh_D13G0434-Gh_D13G0443
|
The numbers after CC and DVC represent different environments. CC has 6 ambient repetitions marked 1–6; DVC has 5 ambient repetitions marked 1–5. The same as below
To further screen important candidate genes from colocalized QTNs, SNPs with stable occurrence and strong signals in multiple environments were examined in the EMMAX model. Two regions were detected, including the significant SNPs provided in Table S5 and S6. The notable targeting areas were on Chr A10, around 5.44 Mb, which was detected in two traits, and on Chr D09, between 6.06 and 6.27 Mb, which was detected in three traits (Table 3). LD decay (r2) was analyzed to identify associated SNPs located in one LD (r2 > 0.6) region. LD regions were observed on Chr A10 and D09.
Table 3
Significant SNPs (− log10(p) ≥ 4) associated with the tolerance trait presenting more than one marker, including physical location and gene range
Trait
|
Chr
|
Start (bp)
|
End (bp)
|
Gene region
|
DW-DVC5
|
A10
|
5430537
|
5442083
|
Gh_A10G0493-Gh_A10G0509
|
DW-CC6
|
A10
|
5430537
|
5499846
|
Gh_A10G0493-Gh_A10G0505
|
ER-CC4
|
D09
|
6245032
|
6269214
|
Gh_D09G0194-Gh_D09G0201
|
TL-CC4
|
D09
|
6390331
|
6455466
|
Gh_D09G0200-Gh_D09G0212
|
DW-CC2
|
D09
|
6062602
|
6273073
|
Gh_D09G0186-Gh_D09G0201
|
Two LD blocks were identified on Chr A10 and D09 for haplotype analysis to characterize the genetic basis of chilling tolerance in our population. The LD block between the 5.43 and 5.44 Mb genome regions on Chr A10 was significantly linked with DW in the two environments (Fig. 5a-b). A total of 13 genes were identified in the LD block. Gh_A10G0500 was selected as the candidate gene because its transcript levels were highly divergent (Table S7). Gh_A10G0500 was 98.43% similar to AT5G07990 in Arabidopsis, involved in flavonoid biosynthetic process, and responsed to UV and auxin. Before the cold treatment (0 h), the gene Gh_A10G0500 was highly expressed in the cold-sensitive variety, with 5.2 times more expression than that in the cold-resistant variety. After 6 and 12 h of chilling stress, the expression levels decreased; however, they were still approximately four times higher than those of the cold-resistant variety. The physical positions of the four SNPs in this candidate gene were 5439916 bp (C/G), 5440164 bp (A/G), 5440494 bp (T/A), and 5440711 bp (G/A). 5440711 bp was in the intron region and the other three were in the exon region, two of which were non-synonymous SNPs. The position 5440164 bp resulted in amino acid substitutions from Lys (K) to Arg (R), and the position 5440494 bp resulted in amino acid substitutions from Met (M) to Lys (K) (Fig. 5d). The cultivars in this study carried three haplotypes for these three SNPs: A (n = 108), B (n = 33), and C (n = 32) (Fig. 5c). Gh_A10G0500 may be the key gene that regulates the dry weight under chilling conditions.
HapC had a higher dry weight than HapA and HapB groups, indicating that HapC was a favorable rare haploid type for better growth under stress. Under DVC condition, HapC could increase the DW by 6.69–8.13% compared with the other two haplotypes. Simultaneously, it could increase the DW by 6.32–9.22% under the CC condition. (Fig. 5e). The number of accessions carrying HapC was 16.00% of the total accessions, and Groups I and IV in HapC were 25.00% and 75.00%, respectively. For haplotypes A, B, and C, the percentage of Group IV increased gradually, while in the other three groups it gradually decreased. Accessions from Groups II and III were not included in the HapC group. From the perspective of the evolutionary tree classification, it was further verified that HapC had stronger cold resistance than the other groups.
In addition, the pleiotropic regions associated with DW-CC, ER-CC, and TL-CC were located on Chr D09, indicating that the region had an important regulatory role under the CC condition. Two SNPs were significantly associated (-log10(p) ≥ 6.41), one located in the linkage region LD09-1 (6.07–6.15 Mb), and the other in the region LD09-2 (6.20–6.27 Mb) (Fig. 6a-b). In the first LD block (LD09-1), there were two haplotypes: HapA-1 had more accessions (n = 163) than the group HapB-1 (n = 16) (Fig. 6c). In the second LD block (LD09-2), HapA-2 had fewer accessions (n = 75) than HapB-2 (n = 108). There were 15 genes in the two regions, 11 genes were located in the LD09-1, and four genes were located in the LD09-2. Combining the low-temperature transcriptome data between HapA (XLZ84) and HapB (XLZ16), Gh_D09G0189 was significantly differentially expressed in the two materials, and the expression of the remaining genes was almost zero or not significantly different (Table S8, Fig. S1). Gh_D09G0189 had the highest homology with AT5G63980 in Arabidopsis which is negatively regulated by signal transduction, is related to the circadian rhythm and responsive to cold temperatures. Based on the name of the AT5G63980, Gh_D09G0189 was named GhSAL1. Gh_D09G0189 was found to have base alterations in the promoter (+ 2000 bp) and coding sequences of the two haplotypes (Table S9). Haplotype analysis showed that the promoter and coding regions divided the 200 materials into two haplotypes, and the material types in haplotype were similar.
Fifteen SNPs in the gene were significantly associated with traits, seven of which were in exons and eight in introns. Three non-synonymous SNPs were located at positions 6069983 bp (A/G), 6070038 bp (T/A), and 6070257 bp (T/A), which resulted in a change in amino acids from Glu (E) to Lys (K), Glu (E) to Val (V), and Thr (T) to Val (V), respectively. Four synonymous polymorphisms were identified: 6068075 bp (T/G), 6069997 bp (G/A), 6070384 bp (T/A), and 6070692 bp (T/C) (Fig. 6e). The cultivars carried two haplotypes for these 15 SNPs, the common (n = 163) HapA and the rarer (n = 16) HapB (Fig. 6c). DW, ER, and TL in the two haplotypes considerably differed under CC stress. GhSAL1HapB was an excellent haplotype, which could increase seedling emergence by 19.04% compared to GhSAL1HapA. After seedling emergence, the dry matter weight and total length increased by 11.26% and 7.69%, respectively (Fig. 6e), suggesting that Gh_D09G0189 may have different regulation methods at different developmental stages and under different low-temperature conditions.
GhSAL1 negatively regulates cold tolerance in cotton
These results suggest that GhSAL1 was the most likely candidate gene associated with cold response in cold-tolerant accessions. To further evaluate the role of GhSAL1, we constructed a gene silencing line using a VIGS system. The survival rate of GhSAL1 silenced lines was significantly higher than that of the negative control by 46.33% after the CC treatment (Fig. 7b–e). Gh_D09G0189 encoded a bifunctional nucleotide phosphatase with 3'(2'),5'-bisphosphate nucleotidase and inositol polyphosphate 1-phosphatase activities, which was highly specific for its substrates PAP and IP3, where it removed the 1-phosphate group from the IP3 second messenger. PAP showed no significant difference between the gene silencing line and negative control before and after the chilling treatment (Fig. 7f). IP3 showed significant differences even before treatment. After 6 h of CC treatment, the difference was 2.12 times greater, thus highly significant. Therefore, we hypothesized that GhSAL1 mainly functions as an inositol polyphosphate 1-phosphatase that mediates the catabolism of IP3 under CC treatment.