Genome-wide association analysis reveals genetic variations and candidate genes associated with salt tolerance related traits in Gossypium hirsutum

Abstract

Background

Cotton is more resistant to salt and drought stresses as compared to other field crops, which makes itself as a pioneer industrial crop in saline-alkali lands. However, abiotic stresses still negatively affect its growth and development significantly. It is therefore important to breed salt tolerance varieties which can help accelerate the improvement of cotton production. The development of molecular markers linked to causal genes has provided an effective and efficient approach for improving salt tolerance.

Results

In this study, a genome-wide association study (GWAS) of salt tolerance related traits at seedling stage was performed based on 2 years of phenotype identification for 217 representative upland cotton cultivars by genotyping-by-sequencing (GBS) platform. A total of 51,060 single nucleotide polymorphisms (SNPs) unevenly distributed among 26 chromosomes were screened across the cotton cultivars, and 25 associations with 27 SNPs scattered over 12 chromosomes were detected significantly (−log₁₀^p > 4) associated with three salt tolerance related traits in 2016 and 2017. Among these, the associations on chromosome A13 and D08 for relative plant height (RPH), A07 for relative shoot fresh matter weight (RSFW), A08 and A13 for relative shoot dry matter weight (RSDW) were expressed in both environments, indicating that they were likely to be stable quantitative trait loci (QTLs). A total of 12 salt-induced candidate genes were identified differentially expressed by the combination of GWAS and transcriptome analysis. Three promising genes were selected for preliminary function verification of salt tolerance. The increase of GH_A13G0171-silenced plants in salt related traits under salt stress indicated its negative function in regulating the salt stress response.

Conclusions

These results provided important genetic variations and candidate genes for accelerating the improvement of salt tolerance in cotton.

Background

The competition for arable land between food crops and cotton (Gossypium spp.) has existed for a long time in China. However, the emphasis on food security has inadvertently moved cotton production to more marginal soils in saline-alkali areas and coastal beaches due to the growing population. Cotton, the largest source of textiles fiber in the world is more resistant to abiotic stresses such as high salt and drought stresses than other crops. However, excessive salt in the soil can severely affect the growth and development of cotton plants [1], resulting in a reduction in fiber yield by as much as 60% [2]. Therefore, breeding cotton varieties with improved salt tolerance could alleviate the conflict between food crops and cotton by reclaiming and utilizing saline-alkali coastal lands for production. Similarly, in the northwestern inland cotton production region of China, the availability of salt tolerant varieties would expand the area of cotton production by promoting the synchronous growth of food crops and cotton.

The genetic architecture of salt tolerance is one of the most important subjects in plant science. Salt tolerance is a complex quantitative trait, which is controlled by multiple genes and involves a variety of physiological and biochemical metabolic pathways in cotton. In addition, the expression of each gene is sensitive to external environment. In the current study, the main conventional breeding approaches for development of salt tolerant varieties are screening and collecting salt tolerant germplasm resources, then transferring elite loci by hybridization, composite hybridization and backcrossing methods. The progress of conventional breeding of salt tolerance in upland cotton is impeded due to the lack of high salt tolerant genetic resources [3], low heritability [4], genetic complexity, and the difficulties in phenotyping [5].

The development of molecular markers linked to causal genes for a trait has provided an effective and efficient approach for improving quantitative traits. Once identified, markers linked to a quantitative trait locus (QTL) such as for salt tolerance can then serve as a selection tool for rapid and efficient marker-assisted selection (MAS). QTL mapping via bi-parental populations is an important method for quantitative trait research, and has been widely employed to map a number of traits including salt tolerance in various crops [6]. In addition to QTL mapping, association mapping based on linkage disequilibrium (LD) is another approach for detecting molecular markers tightly linked with quantitative traits in a natural germplasm population. Association analysis is time- and cost-effective, and can mine for genetic variations existed in the natural population, and more importantly, it takes advantage of the recombination information that has accumulated in the natural population over the long-term evolution process, thus achieving a higher mapping resolution, possibly identifying the causative genes [7, 8]. Because of its effectiveness in QTL identification, association mapping has been widely used in crop species and plays an important role in molecular breeding [9,10,11,12,13,14].

Genotyping-by-sequencing (GBS) is a reduced representation genotyping platform, and has a broad application in crop genetics and breeding [15,16,17,18]. In general, the GBS approach begins with an enzymatic digestion to reduce genome complexity using barcoding restriction enzymes, then performs multiple sequencing of barcoding DNA fragments on the high-throughput next-generation sequencing platform. A bioinformatics analysis of indexed sequence reads is followed to identify genetic variants. Finally, genetic diversity analysis is carried out based on a sample-by-variant matrix [19]. GBS is a powerful and cost-effective tool to assess variations across populations.

Since the first molecular marker-based genetic map of cotton published in 1994 [20], cotton scientists have identified a large number of QTLs regulating important agronomic traits, including fiber quality, yield, and disease resistance. However, only limited reports have been published on salt tolerance in cotton [21,22,23,24,25]. Herein, we reported the genome-wide association study (GWAS) of salt tolerant QTLs during the seedling stage performed over 2 years of phenotyping on 217 representative upland cotton cultivars. These results provided important genetic variations and candidate genes for accelerating the improvement of salt tolerance in cotton.

Results

Phenotypic diversity analysis

In order to evaluate the phenotypic variations of salt tolerance in the GWAS population with 217 upland cotton cultivars (Additional file 1: Table S1), three traits related to salt tolerance including relative plant height (RPH), relative shoot fresh matter weight (RSFW) and relative shoot dry matter weight (RSDW) were determined. The mean values, ranges, standard deviations (SD), coefficients of variation (CV) and broad-sense heritability (H²) for these salt tolerance related traits are shown in Table 1. Great differences of the CV in both years were found for the three traits. Overall, the upland cotton cultivars in this GWAS panel clearly exhibited considerable natural variations in the three traits to salt tolerance and displayed very high genetic diversity.

Table 1 Statistics of various traits related to salt tolerance

Genetic diversity analysis

GBS produced 53,321 high-quality polymorphic single nucleotide polymorphisms (SNPs) among the 217 upland cotton cultivars, with 47,133 SNPs located in the intergenic intervals and 6188 SNPs in the coding regions. Of these, 95.8% of the loci (51,060 SNPs) were mapped onto 26 chromosomes of the cotton genome, and were selected for the GWAS analysis (Additional file 2: Table S2). These SNPs were not evenly distributed (Additional file 3: Fig.S1) with an average of 1964 SNPs per chromosome ranging from 863 to 5209. Chromosome A08 had the most SNPs (5209) with the highest SNP density (199 kb/SNP). D04 had the least SNPs (863), but A02 had the lowest SNP density with 680 kb/SNP. The genetic diversity of this population varied from 0.224 (A08) to 0.389 (A13). The polymorphism information content (PIC) varied from 0.192 (A08) to 0.305 (A13) (Table 2). The above results indicated a relatively large span in genetic diversity index and PIC in the cotton genome.

Table 2 Summary of polymorphic SNPs mapped in 26 chromosomes of Gossypium hirsutum

Population structure and LD analysis

It is important for GWAS analysis to control the effect of population structure, because population stratification could eliminate spurious associations between genotypes and phenotypes [26, 27]. STRUCTURE software was used to calculate the Bayesian clustering from K = 1 to 10 for five repetitions. LnP (D) value continued to increase from K = 1 to K = 10 without a significant inflection point (Fig. 1a). However, there was an obvious spike at the value of ΔK = 3 (Fig. 1b), suggesting that the population could be divided into 3 subgroups (Fig. 1c). Taking the corresponding Q matrix at k = 3 as the covariate could reasonably eliminate spurious association effects and improve the GWAS accuracy.

Fig. 1

Population structure analysis of 217 cotton cultivars. a Estimated LnP(K) of possible clusters (K) from 1 to 10. b Delta K based on the rate of change of LnP(K) between successive K. c Population structure of 217 upland cotton accessions based on STRUCTURE when K = 3

The LD distribution among chromosomes of the 217 upland cotton cultivars was shown in Table 2. The LD decay distance was 869 kb and 1120 kb when the r² dropped to 0.2 and 0.1, respectively. The LD decay distance was not evenly distributed among chromosomes, ranged from 435 kb (D01) to 1157 kb (A06) when r² was set at 0.2 and ranged from 712 kb (D02) to 1377 kb (A06) when r² was set at 0.1, the overall LD decay in the At subgenome was significantly higher than that in the Dt subgenome.

Association analysis

In order to explore the genetic factors underlying salt tolerance, the mixed linear models (MLMs) were performed by simultaneously accounting for population structure and relative kinship matrix to conduct a GWAS. A total of 25 significant associations (−log₁₀^p > 4) with 27 significant SNPs located on chromosomes A05, A07, A08, A09, A10, A11, A12, A13, D02, D03, D06, and D09 were detected for the three salt tolerance related traits in the 2016 and 2017 dataset. Eleven associations with 12 SNPs were detected in 2016 (Fig. 2) and nine associations with 9 SNPs were detected in 2017 (Fig. 3). In addition, five associations with 6 SNPs were detected in both 2016 and 2017 dataset. The phenotypic variance explained (PVE) by individual QTL ranged from 1.29 to 7.00% (Table 3).

Fig. 2

Manhattan and Q–Q plots for salt tolerance related traits in 2016. The horizontal dotted lines of the Manhattan plots with black color represent the genome-wide significance threshold of 0.0001

Fig. 3

Manhattan and Q–Q plots for salt tolerance related traits in 2017. The horizontal dotted lines of the Manhattan plots with black color represent the genome-wide significance threshold of 0.0001

Table 3 Summary of SNPs significantly associated with salt tolerance related traits

For RPH, nine significant associations with 10 SNPs were identified on chromosomes A07, A09, A10, A12, A13, D02, D03, and D08. The PVE ranged from 1.50 to 7.00%. Moreover, the association with one SNP on chromosome A13 and the association with two SNPs on chromosome D08 were detected in both years.

For RSDW, eight significant associations with 9 SNPs on A05, A08, A13, and D08 were detected. The PVE ranged from 1.29 to 4.85%. The association with one SNP on chromosome A08 and the association with one SNP on chromosome A13 were detected in both years.

For RSFW, eight associations with 8 significant SNPs on A07, A08, A10, A11 and D06 were detected in 2016 and 2017. The PVE ranged from 2.35 to 5.25%. The association with three SNPs on chromosome A07 was detected in both years.

Pleiotropy was also found in our GWAS results. For example, the significant SNP A07_90,682,411 was simultaneously detected to be associated with RPH in 2017 and RSFW in both environments. The significant SNP A10_84,786,908 was simultaneously detected to be associated with RPH in 2016 and RSFW in 2016. The significant SNP D08_49,014,753 and D08_49,080,865 simultaneously detected to be associated with RPH in 2016, 2017, and RSDW in 2016. In addition, the associations on A08 for RSFW and RSDW in 2016 were only 70 kb apart. Their confidence interval overlapped with each other.

Identification and preliminary function verification of candidate genes

In general, the LD decay could be used as confidence interval to identify candidate genes. Because the cotton genome has a large LD decay [28, 29], we extracted potential candidate genes within 100 kb of flanking significant markers on the basis of the published upland cotton genome sequencing database [30]. A total of 156 genes were identified in these intervals (Additional file 4: Table S3). Gene ontology (GO) enrichment analysis indicated that “dioxygenase activity”, “oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, incorporation of two atoms of oxygen” and “oxidoreductase activity, acting on single donors with incorporation of molecular oxygen” were significantly enriched using an false discovery rates (FDR) adjusted P-value of ≤0.05 as the cutoff. Of the 156 genes, 12 were differentially expressed between salt tolerant variety Miscott7913–83 and salt sensitive variety Su 12 according to the previous transcriptome sequencing results [31] (Additional file 4: Table S3). Some of these genes may be associated with salt stress, such as GH_A08G0488 and GH_A10G1620 encoding protein kinase. Protein kinase has been proved to play an important role in salt tolerance in cotton [32]. Another gene GH_A13G0171, which encodes aquaporins (AQPs), was also likely to regulate the salt stress response [33]. The confidence interval contains GH_A13G0171 was simultaneously detected in both years. We found that the salt tolerance of upland cotton cultivars with the G haplotype was significantly higher than that of cultivars with the C haplotype in both years upon a t test (Fig. 4).

Fig. 4

Box plots for the phenotypic values of the association for RPH on chromosome A13. ** indicate statistical significance at the 0.01 probability level. C: C haplotype. G: G haplotype

The three promising genes (GH_A08G0488, GH_A10G1620 and GH_A13G0171) were selected for preliminary function verification of salt tolerance in cotton. Analysis of gene expression patterns could provide important clues for gene function determination. A quantitative real-time PCR (qRT-PCR) was performed to analyze the expression levels of GH_A08G0488, GH_A10G1620 and GH_A13G0171 in roots and leaves under salt stress treatment in salt tolerant variety Miscott7913–83 and salt sensitive variety Su 12. As shown in Fig. 5, the three genes were induced by salt stress and displayed distinct expression patterns in response to salt stress in salt tolerant variety Miscott7913–83. The three genes had a much higher expression level in roots than in leaves. The gene GH_A13G0171 exhibited a significantly down-regulated expression in both root and leaf tissues after salt stress. The gene GH_A08G0488 exhibited a significantly up-regulated expression in both root and leaf tissues. The expression level of GH_A10G1620 showed an increase in leaf and no significant changes in the root tissues. As shown in Fig. 6, the three genes were also displayed distinct expression patterns in response to salt stress in salt sensitive variety Su 12. The gene GH_A13G0171 exhibited an identical expression pattern in Miscott7913–83 and Su 12. The expression levels of GH_A10G1620 and GH_A08G0488 were not significant different.

Fig. 5

Expression profiles of three promising genes under NaCl stress at three time points in Miscott7913–83. a Expression profiles of GH_A13G0171. b Expression profiles of GH_A10G1620. c Expression profiles of GH_A08G0488. M: Miscott7913–83; Y: Leaf; G: Root. The samples were collected at 0 h, 12 h and 72 h under salt stress. The cotton Actin gene was used as an internal reference. The data are mean ± SE of three biological replicates. * and ** indicate statistical significance at the 0.05 and 0.01 probability level, respectively

Fig. 6

Expression profiles of three promising genes under NaCl stress at three time points in Su 12. a Expression profiles of GH_A13G0171. b Expression profiles of GH_A10G1620. c Expression profiles of GH_A08G0488. S: Su 12; Y: Leaf; G: Root. The samples were collected at 0 h, 12 h and 72 h under salt stress. The cotton Actin gene was used as an internal reference. The data are mean ± SE of three biological replicates. * and ** indicate statistical significance at the 0.05 and 0.01 probability level, respectively

To confirm the functional roles of GH_A08G0488, GH_A10G1620 and GH_A13G0171 genes under salt stress, virus-induced gene silencing (VIGS) assay was used to repress expression of these genes in salt tolerant variety Miscott7913–83 plants. The inoculated seedlings were grown in three light incubators at 23 °C under a 16-h light and 8-h dark cycle as three biological replicates. At the developmental period when two leaves had formed, the pTRV2:: GH_A08G0488, pTRV2:: GH_A10G1620, pTRV2:: GH_A13G0171 and pTRV2:: 00 inoculated plants were treated with 350 mM NaCl. After 15 days, the plant height, fresh and dry shoot matter weight were determined and the corresponding relative values were calculated. The transcripts of the three genes in the VIGS leaves were significantly reduced compared to the negative control pTRV2:: 00 inoculated plants, indicating that they were effectively silenced (Additional file 5: Fig.S2, Fig. 7a). Compared with the control pTRV2:: 00, no significance effect on their phenotypes was observed in pTRV2:: GH_A08G0488 and pTRV2:: GH_A10G1620 inoculated plants, and the plant height, fresh and dry shoot matter weigh of GH_A13G0171-silenced plants were significantly higher than that of pTRV2:: 00 inoculated plants (Table 4, Fig. 7b), indicating that GH_A13G0171 could reduce seedling tolerance to salt stress. The corresponding markers SNP_A13_1,869,056 could be applied for improvement of cotton salt tolerance.

Fig. 7

Functional characterization of GH_A13_G0171 VIGS in Miscott7913–83. a Transcript levels of GH_A13G0171 in pTRV2:: 00 and pTRV2:: GH_A13G0171 inoculated plants. The cotton Actin gene was used as an internal reference. The data are mean ± SE of three biological replicates. ** indicate statistical significance at the 0.01 probability level, respectively. b Phenotype identification of salt tolerance of the VIGS line and its control in seedling stage

Table 4 Statistics of traits related to salt tolerance for the VIGS lines and its control at the seedling stage

Discussion

An ideal GWAS population should encompass sufficient phenotypic and genotypic diversity [34]. Thus, selection of a suitable panel with abundant genetic diversity is especially critical for GWAS analysis. In this study, the population panel consisted of 217 cultivars from different ecological regions (the Yellow River, the Yangtze River, the Northwestern inland region, the Northern China region, America and Australia) Phenotypic analysis showed that the salt tolerance related traits in two environments had wide variations and stable changing trends in different environments. Genotyping analysis revealed an average genetic diversity and PIC of 0.303 and 0.248, respectively. As a result, the 217 cotton accessions with broad geographic and genetic diversity coverage in our research provided sufficient detection power for GWAS analysis.

Recent years, association mapping as an important method for the identification of complex quantitative traits has been widely used in cotton for various traits including salt tolerance [35,36,37,38,39,40]. Jia et al. [21] identified 3 simple sequence repeats (SSR) markers significantly associated with the relative survival rate under salt stress through association mapping methods using 323 G. hirsutum germplasms. Du et al. [22] performed association analysis of 304 upland cotton cultivars and identified 95 significant associations for 10 salt tolerance related traits at the germination and seedling stages. These studies were conducted based on the SSR markers with lower mapping resolution and smaller genome coverage. Therefore, the genetic factors of salt tolerance related traits discovered with high-density SNP markers could provide more precise insight into salt tolerance in cotton than previous studies. Sun et al. [24] performed a GWAS analysis to identify marker-trait associations under salt stress using 713 upland cotton accessions through the Illumina Infinium CottonSNP63K array. A total of 23 SNPs were significantly associated with two salt-tolerance-related traits. All of the above reported associations were only detected in a single environment. QTL mapping of important agronomic traits of cotton has been carried out for nearly two decades, but few successful practicing of breeding new varieties by MAS was reported [41]. An important reason is the insufficient available information on multi-environmental detected QTLs. Salt tolerance as a complex quantitative trait controlled by multiple genes needs repeated phenotypic evaluations in multiple environments. In this study, a total of 25 significant associations with 27 SNPs scattered over 12 chromosomes were detected for salt tolerance related traits. Among these, the associations on chromosome A13 and D08 for RPH, A07 for RSFW, A08 and A13 for RSDW were expressed in both environments, indicating that they were likely to be stable QTLs. These results provided important genetic variations for accelerating the improvement of salt tolerance in cotton.

The comparison between our GWAS results and previous reports showed that most of the 25 associations were novel reported loci. The physical location of the SSR marker NAU2561 on A05 significantly associated with the relative malondialdehyde (MDA) content reported by Du et al. (2016) [22] was only 1.7 Mb from SNP_A05_ 94,406,377, which was significantly associated RSDW in 2016. The physical distance between BNL1231, which was reported by Du et al. (2016) [22] to be associated with relative superoxide dismutase (SOD) activity, and SNP_A11_117,938,337 which was associated with RSFW in 2017 was only 1 Mb. Considering larger LD decay distances in the cotton genome, it is reasonable to presume that these two loci could be the same loci.

Many studies have reported the positive effects of AQPs in improving plant tolerance to salt stresses. In this study, both gene expression and silencing analyses showed that GH_A13G0171, which encodes AQPs, was likely to regulate the salt stress response. AQPs play an important role in maintaining water homeostasis in the responses of plants to environmental stress by regulating water movement across vacuolar membranes [42]. Overexpression of MaPIP2–7 in banana improved tolerance to drought, cold, and salt stresses [43]. Overexpression of SpAQP1 promoted seed germination and root growth under salt stress in transgenic tobaccos [33]. The negative effect of AQPs also has been reported. Wang et al. [44] found that overexpression of GsTIP2;1 in Arabidopsis thaliana depressed salt tolerance and dehydration stress. Plants with overexpressed GsTIP2;1 exhibited higher dehydration rate, suggesting that this gene may mediate stress response by increasing water loss. GH_A13G0171 is homologous to Arabidopsis gene AT4G00430. The transgenic plants overexpressing PIP1;4 (AT4G00430) displayed a rapid water loss under dehydration stress, which resulted in retarded germination and seedling growth under drought stress [45]. Our data showed that GH_A13G0171-silenced plants were significantly improved in salt tolerance related traits, indicating its negative effect in regulating the salt stress response. The exact role of this gene still needs further verification by RNA interference or clustered regularly interspaced short palindromic repeats (CRISPR) technology. Nonetheless, this study adds further credence to the relationship between AQPs and plant stress response.

Preliminary functional verification by VIGS assay suggested that GH_A08G0488 and GH_A10G1620 encoding protein kinase were irrelevant to salt tolerance. However, the real relationships between the two genes and salt tolerance still require verification of transgenic overexpression technology. Protein kinase which is widespread in plants plays an important role in signal transduction and is involved in response to a variety of plant hormones and environmental signals. In order to study its role in conferring abiotic stress tolerance, Zhao et al. [32] introduced cotton GbRLK gene into Arabidopsis thaliana. The results showed that GbRLK was involved in drought and high salinity stress by activating abscisic acid signaling pathway. The overexpression of GbRLK gene may reduce water loss by regulating stress response genes, thus improving stress resistance. He et al. [46] found that the salt-induced gene TaGSK1 from wheat was a highly conserved serine/threonine protein kinase. Overexpression of TaGSK1 in Arabidopsis thaliana decreased the osmotic pressure and Na⁺ content, and improved the salt tolerance. Overall, further molecular and functional analysis of GH_A08G0488 and GH_A10G1620 would advance our understanding of the molecular mechanisms underlying salt stress.

Conclusions

It is of great significance to explore the genetic variants and the underlying candidate genes in upland cotton for salt tolerance improvement. In this study, we reported the GWAS analysis of salt tolerant QTLs during the seedling stage by GBS platform. A total of 25 associations with 27 SNPs scattered over 12 chromosomes were detected significantly associated with three salt tolerance related traits in two environments, in which five associations were simultaneously expressed in both environments. The VIGS assay revealed that GH_A13G0171, which encodes AQP PIP-type, negatively to regulates the salt stress response. All of our findings consistently suggest the identified candidate genes would be useful in salt tolerance breeding in cotton upon further validation with direct functional analyses which are underway.

Methods

Plant materials

A collection of 217 upland cotton cultivars (Gossypium hirsutum L.) were selected to assemble a GWAS panel in this study. These cultivars were collected from the Institute of Industrial Crops, Jiangsu Academy of Agricultural Sciences, kindly supplied by the Institute of Cotton Research, Chinese Academy of Agricultural Science (CRI-CAAS) and Cotton Research Center of Shandong Academy of Agricultural Sciences. Of these upland cotton cultivars, 162 were collected from China and 55 were introduced from other countries. Chinese cultivars were from four ecological areas, i.e. the Yellow River (73), the Yangtze River (62), the Northwestern inland region (23), and the early maturation area in Northern China (4).

In 2015, the 217 upland cotton cultivars were planted at the Lishui plant experiment station, Jiangsu Academy of Agricultural Sciences, China. Genomic DNA extraction referred to the method described by Paterson et al. (1993) [47]. The central natural open bolls of the 217 upland cotton cultivars were hand-harvested and ginned. The delinted seeds were then used for the subsequent determination of salt tolerance assay.

Investigation of salt tolerance related traits

The cotton seeds of each cultivar were surface sterilized, transferred into 9 cm sterile Petri dishes on filter paper, and then wetted with 35 mL distilled water for germination. The Petri dishes were incubated in a growth chamber at 32 °C with no light. The germinated seeds were then grown in cups of dimension 220 ml with equal amounts of 1:1 mix of the matrix and vermiculite, and placed in 57 cm × 41 cm plastic pallets in a greenhouse. At the developmental period when two leaves had formed, the seedlings for each cultivar with the same growth vigor were divided into three groups with ten plants each. The first ten seedlings with two replications per cultivar were selected before treatment, and their plant height, fresh and dry shoot matter weight were determined as the 0-day control; the second ten seedlings with two replications were selected for 350 mM NaCl treatment, and the third ten seedlings with two replications were treated with clean water. After 7 days, the plant height, fresh and dry shoot matter weight of NaCl and clean water treated seedlings were determined respectively. The phenotype increments under 350 mM NaCl in 7 days was calculated using the following formula: (phenotypic effects value under salt stress after 7 days - phenotypic effects value at 0 day). The phenotype increments under clean water in 7 days was calculated using the following formula: (phenotypic effects value under well-watered control after 7 days - phenotypic effects value at 0 day). The ratio of the phenotype increments under 350 mM NaCl and clean water was considered as the relative values. The phenotype of three salt tolerance related traits were calculated using the following formula: (phenotype increments under 350 mM NaCl / (phenotype increments under clean water) × 100%.

We analyzed three traits related to salt tolerance in 2 years, including RPH, RSFW and RSDW. The average trait values in each year were used for phenotypic analysis and GWAS analysis in a single environment. Statistical analysis of the phenotype of salt tolerance related traits was performed by SPSS. H² of each trait was calculated from variance components. The REML method of PROC VARCOMP in SAS/STAT software (SAS Institute Inc., Cary, NC, USA) was conducted to estimate the variance components.

Library preparation and Illumina sequencing for GBS analysis

Firstly, a GBS pre-design experiment was performed. The type of enzymes and their fragments size were evaluated based on three criteria. (1) The number of tags must be appropriate. (2) The enzymatic tags must be evenly distributed. (3) Repeated tags must be avoided. Tags about 50 bp was selected to maintain the sequence depth uniformity of different fragments. The GBS library was constructed in accordance with the pre-designed scheme and referred to the method described by Fan et al. (2018) [48].

Then, pair-end sequencing was performed on the selected tags using an Illumina high-throughput sequencing platform Illumina NovaSeq. We generated a total of 5927 Gb high quality data with an average sequencing depth of approximately elevenfold. The sequences were sorted according to the barcodes. To make sure that reads were reliable and without artificial bias (low quality paired reads, which mainly resulted from base-calling duplicates and adapter contamination) in the following analyses. C scripts were conducted to processed raw fastq format reads through a series of quality control procedures as followed: (1) removing reads with unidentified nucleotides (N) ≥10%; (2) removing reads with > 50% bases having phred quality < 5; (3) removing reads aligned to the adapter more than 10 nucleotides.

SNP detection and annotation

The clean reads were anchored to the cotton reference genome [30] using Burrows-wheeler aligner (BWA) [49]. The SAM tools software [50] was used to convert alignment files to BAM files. After alignment, the genomic variants (in GVCF format for each accession) were identified with the Haplotype caller module and the GVCF model by Genome analysis toolkit (GATK) software [51].. SNP was filtered by the perl script called ANNOVAR [52].

Population structure and LD analysis

The population structure of our mapping panel was inferred from the GBS data with STRUCTURE software [53]. Five independent runs were performed; the number of populations (K) was set from 1 to 10; the burn-in time and Markov-chain Monte Carlo replication numbers was set to 10,000. The optimal K value was determined by comparing the LnP (D) and Δk based on the rate of change in LnP (D) [54]. A Q-matrix produced by STRUCTURE listed the estimated membership coefficients in a cluster for the subsequent association analysis. The relative kinship matrix was calculated by SPAGeDi [55]. The software TASSEL 3.0 package was employed to estimate the LD parameters (r²) of adjacent loci within the same LD chromosome for each pair of SNP loci [56].

Association mapping

The TASSEL 3.0 software package [56] was employed to construct association tests of salt tolerance related traits. The MLMs was performed by simultaneously accounting for multiple levels of Q-matrix and K-matrix according to the methods described by Yu et al. (2006) [57].

Expression analysis of salt-induced genes

Total RNA was extracted using an improved CTAB method [58]. For qRT-PCR, PrimeScript RT Master Mix (Perfect Real Time) from TaKaRa (Dalian, China) was used to synthesize first-strand cDNA. The qRT-PCR program was performed on an ABI QuantStudio 5 RealTime PCR System. The qRT-PCR reaction system referred to the method described by Xu et al. (2017) [59]. qRT-PCR was carried out with three biological replicates. The cotton Actin gene was employed as a reference gene. The specific primers for qRT-PCR are listed in Additional file 6: Table S4.

Virus-induced gene silencing assay

About 300–500 bp fragments overlapped the C-terminal and 3′ untranslated region (UTR) from 3 differentially expressed genes (GH_A08G0488, GH_A10G1620 and GH_A13G0171) were isolated in Miscott7913–83. The specific primers were listed in Additional file 6: Table S4. These fragments were sequenced and cloned into double enzyme (EcoRI and BamHI) digested pTRV2, generating a new vector named pTRV2:: gene. The help plasmid pTRV1, pTRV2:: gene, VIGS positive control pTRV2:: CLA1 and negative control pTRV2:: 00 (empty vector) were introduced into Agrobacterium tumefaciens strain GV3101 by heat shock method. The VIGS assay was performed according to the protocol described by Gao et al. (2011) [60]. A ratio of 1:1 mixed pTRV1 and pTRV2:: gene were injected into the seedlings with mature cotyledons but without a visible true leaf via a syringe infiltration. The inoculated seedlings were grown in a light incubator at 23 °C under a 16-h light and 8-h dark cycle. Two weeks later, the inoculated plants RNA was isolated to detect whether the target gene was still expressed. For each gene, a total of 100 seedlings were knocked out for subsequent salt tolerance identification.