High-Density Linkage Map Construction and Mapping of Salt-Tolerant QTLs at Seedling Stage in Upland Cotton Using Genotyping by Sequencing (GBS)

Over 6% of agricultural land is affected by salinity. It is becoming obligatory to use saline soils, so growing salt-tolerant plants is a priority. To gain an understanding of the genetic basis of upland cotton tolerance to salinity at seedling stage, an intra-specific cross was developed from CCRI35, tolerant to salinity, as female with Nan Dan (NH), sensitive to salinity, as the male. A genetic map of 5178 SNP markers was developed from 277 F2:3 populations. The map spanned 4768.098 cM, with an average distance of 0.92 cM. A total of 66 QTLs for 10 traits related to salinity were detected in three environments (0, 110, and 150 mM salt treatment). Only 14 QTLs were consistent, accounting for 2.72% to 9.87% of phenotypic variation. Parental contributions were found to be in the ratio of 3:1, 10 QTLs from the sensitive and four QTLs from the resistant parent. Five QTLs were located in At and nine QTLs in the Dt sub-genome. Moreover, eight clusters were identified, in which 12 putative key genes were found to be related to salinity. The GBS-SNPs-based genetic map developed is the first high-density genetic map that has the potential to provide deeper insights into upland cotton salinity tolerance. The 12 key genes found in this study could be used for QTL fine mapping and cloning for further studies.


Introduction
Currently, it is estimated that over 6% of agricultural land is affected by salinity [1]. To satisfy the increasing population demand, it is becoming obligatory to use saline soils, either by reclamation or by growing salt-tolerant plants [2]. To mitigate the challenge of salt stress, selective breeding programs are adopted in order to improve the performance of the genotypes developed [3]. However, the majority of traits with economic significance, such as salinity tolerance, are regulated by quantitative trait loci (QTL) and thus are environment-dependent [4]. Using conventional breeding techniques to enhance genetic improvement has limitations due to low efficiency, slow speed, and being expensive for some traits like salt tolerance [4]. Molecular breeding, such as marker-assisted breeding (MAS) and genomic selection (GS), have the potential to overcome the inefficiencies of conventional breeding [5]. This is equivalent to approximately 83.13-foldhaploid genome coverage of raw paired-end Illumina reads by sequencing whole-genome shotgun (WGS) libraries of homozygous cv. "TM-1" compared to Li Fuguang et al. (2015) [29] in their study, which generated a total of 445.7 Gb of clean reads, or 181-fold haploid genome coverage, of raw paired-end Illumina reads by sequencing whole-genome shotgun (WGS) libraries of homozygous cv. "TM-1" had fragment lengths ranging from 250 bp to 40 kb. The average GC content of the sequences is 38.25%, with a Q20 score of 94.66%.

SNP Detection and Annotation
The parents, CCRI35 and NH, were heterozygous and homozygous lines with AC and AA genotypes, respectively. The genotype AC × AA, consisting of 93,384 markers, was used for further analysis. Among the 93,384 markers, the low-coverage sequences of the F 2:3 populations (coverage less than 75%) were filtered out, leaving 24,549 markers. Markers with significant distortion (p < 0.001) were filtered out and 6405 markers were retained with the purpose of determining bin markers.

Phenotypic Variation between Parents
There was a wide range of phenotypic variation among the F 2:3 populations in all the traits measured and in both environments: MDA, EC; GR; FW; SL; LFW; SLW; DLW; RWC and CHL Figure 1. All traits exhibited a normal segregation pattern, with equal distribution ( Figure S1). In the control, there was no significant difference observed between the parents, except for chlorophyll content (CHL). However, under salt treatment in all environments, whether 110 mM or 150 mM, all traits were significantly reduced in the salt-susceptible genotype compared to the resistant parent; MDA, LFW, SLW, GR, and SL showed a significant difference. However, in MDA, higher concentration denoted sensitivity; its concentration in a resistant accession was significantly lower than the sensitive accession ( Figure 1).

ANOVA and Heritability Analysis of Salt Tolerance Traits for the Two Parents and the Progeny
All the data collected were analyzed for ANOVA using mixed-model analysis. The ANOVA results revealed significant differences between the genotypes, environment, and their interactions for all the traits ( Table 1).
The measured trait's heritability percentage was much higher on the morphological traits compared to the physiological traits. The highest percentage of heritability was achieved in fresh weight (FW), at 87.4%, while the lowest level of heritability was noted in relative water content (RWC), at 42.3%.

Correlation Analysis
Correlation analysis was carried out using the mean values of the five seedlings in each replicate for the two salt treatments, 110 mM and 150 mM.

Salt Treatment (110 mM)
Under 110 mM, a greater percentage of the traits were positively correlated. However, MDA had no significant correlation to any of the measured traits. Negative correlation was noted on EC/DLW, EC/CHL, SL/CHL, and RWC/CHL, while no significant correlation was observed between EC/FW, EC/LFW, EC/SLW, EC/RWC, GR/RWC, GR/CHL, FW/CHL, SLW/RWC, and DLW/RWC (Figure 2a).

Salt Treatment (150 mM)
In an increased salt concentration of 150 mM, MDA was negatively correlated to all traits except for EC, but EC was negatively correlated to GR and DLW. Moreover, GR, FW, SL, LFW, SLW, DLW, RWC, and CHL were positively correlated. However, no correlation was noted between RWC/GR, EC with FW, LFW, RWC and CHL; DLW/SL (Figure 2b).

Construction of the Linkage Maps
A total of 6405 markers and phenotypic data of the F 2:3 , an intra-specific cross of two upland cottons, were utilized for generating the intra-specific linkage map. The genotyping data of 5178 GBS markers were used for mapping after removing the distorted and the unlinked markers in Join Map 4.0 (Table S1). Linkage maps of the 5178 GBS markers generated 26 linkage groups (Figure 3 and Figure S2 and Table 2). The 26 LGs were designated as A01 to A13 for the A t sub-genome and D01 to D13 for the D t sub-genome. The total map distance was 4768.098 cM, higher than the most current linkage map with a map distance of 4450 cM of cotton genome [30]. The map developed in this study is the most highly saturated intra-specific map of upland cotton ever developed. The average distance between adjacent markers was 0.92 cM. The A t sub-genome spanned 2611.43 cM and consisted of 3313 markers with 13 linkage groups. The average distance in A t sub-genome was 0.79 cM, with a maximum gap of 26.598 cM between adjacent markers. In the D t sub-genome, 13 linkage groups were assigned, comprised of 1865 markers spanning 2156.67 cM, with an average of 1.156 cM. The maximum gap between adjacent loci was 30.082 cM ( Table 2). Chromosomes; A02, D02, A01, A05, A03, D01 and A10 had more markers compared to other chromosomes such as D06 and D13. For instance, Chr A02 had 705 loci with map distance of346.314 cM with average distance of 0.49 cM. Chr D06 was the smallest with only 16markers, and a total length of 79.084 cM.

Salt Treatment (150 mM)
In an increased salt concentration of 150 mM, MDA was negatively correlated to all traits except for EC, but EC was negatively correlated to GR and DLW. Moreover, GR, FW, SL, LFW, SLW, DLW, RWC, and CHL were positively correlated. However, no correlation was noted between RWC/GR, EC with FW, LFW, RWC and CHL; DLW/SL (Figure 2b).

Construction of the Linkage Maps
A total of 6405 markers and phenotypic data of the F2:3, an intra-specific cross of two upland cottons, were utilized for generating the intra-specific linkage map. The genotyping data of 5178 GBS markers were used for mapping after removing the distorted and the unlinked markers in Join Map 4.0 (Table S1). Linkage maps of the 5178 GBS markers generated 26 linkage groups (Figures 3 and S2 and Table 2). The 26 LGs were designated as A01 to A13 for the At sub-genome and D01 to D13 for the Dt sub-genome. The total map distance was 4768.098 cM, higher than the most current linkage map with a map distance of 4450 cM of cotton genome [30]. The map developed in this study is the most highly saturated intra-specific map of upland cotton ever developed. The average distance between adjacent markers was 0.92 cM. The At sub-genome spanned 2611.43 cM and consisted of 3313 markers with 13 linkage groups. The average distance in At sub-genome was 0.79 cM, with a maximum gap of 26.598 cM between adjacent markers. In the Dt sub-genome, 13 linkage groups were assigned, comprised of 1865 markers spanning 2156.67 cM, with an average of 1.156 cM. The maximum gap between adjacent loci was 30.082 cM ( Table 2). Chromosomes; A02, D02, A01, A05, A03, D01 and A10 had more markers compared to other chromosomes such as D06 and D13. For instance, Chr A02 had 705 loci with map distance of346.314 cM with average distance of 0.49 cM. Chr D06 was the smallest with only 16markers, and a total length of 79.084 cM.    Ratio: number of markers less than (<) 5 cM divided by total number of markers within chromosome.

Clustered QTLs Region
QTL clusters are regions of the genome in which large quantities of QTLs are co-localized [31]. A total of eight clusters for six traits were detected. The highest number of consistent QTLs mapped was three and all were identified in the marker interval of mk12058_D03-mk12186_D03on D03 (c17), as in Figure 4. This region was designated as Cluster 5, from 3,459,294 bp to 42,120,184 bp. Cluster 5 harbored QTLs for DLW, FW, and SLW with the following proportions: 4, 5, and 5, respectively, which explained the phenotypic variance ranging from 2.72% to 9.87% (Table 3 and Table S2). The lowest number of major QTLs was identified in Clusters 1, 7, and 8, which harbored QTLs for EC, FW, and MDA with the following proportions: 5, 3, and 2, respectively. Furthermore, we focused our analysis on the 14 stable QTLs across multiple environments based on their broad-sense heritability and phenotypic variation explained by a single QTL. These QTLs were mapped in eight chromosomes; A02, A06, A12, D01, D03, D06, D08, and D13, with five, six, five, eight, 14, four, three, and two QTLs, respectively (see Table 3 and Figure S2).
A total of 66 QTLs were identified for 10 traits, but only 14 QTLs were consistent in at least two environments (Table S2 and Figure S2). The distribution of the QTLs within the identified chromosomes is illustrated in Table 3. Of the 14 consistent detected QTLs, five were located on the A t sub-genome while the remaining nine were located on the D t sub-genome (Table 3). As to the contributions of the parents toward the QTLs, 10 QTLs were contributed by the sensitive parent (NH), while only four QTLs were contributed by the resistant parent (CCRI35). Moreover, only eight chromosomes out of 26 were found to harbor consistent QTLs for six traits (EC, FW, SLW, SL, DLW, and MDA) related to salt tolerance. Four types of gene actions were revealed by the genetic effects, of which four genes exhibited dominant effects (D), 13 partial dominance (PD), 28 overdominance (OD), while only two had an additive effect (A). OD was observed for most of the traits in response to salt tolerance, as shown in Table 3.

Identification of Candidate Genes
The putative candidate genes identification was done based on the QTLs consistency. The overall number of consistent QTLs was grouped into eight clusters with varying marker flanking intervals: Cluster 1 (mk19866 to mk1778_A02); Cluster 2 (mk4886_A06 to mk5000_A06); Cluster 3 (mk18878 to mk9189_A12); Cluster 4 (mk10837_D01 to mk10918_D01); Cluster 5 (mk12058_D03 to mk12186_D03); Cluster 6 (mk17780 to mk13562_D06); Cluster 7 (MulMa512_D08 to mk16015_D08), and finally Cluster 8, which had a marker interval from mk11123 to mk18547_D13. Clusters 1 to 8 were located on chromosomes A02, A06, A12, D01, D03, D06, D08, and D13, respectively. The highest number of stable QTLs was detected in Cluster 5 (D03), with three QTLs linked to DLW, FW, and SLW (see Table S3 and Figures 5 and 6).      In all the clusters, 5596 genes were mined, of which 176 were found to be related to biotic and abiotic stress. The stress-related genes were found to belong to the NAC, WRKY, HAD, C2H2, and RING/U-box superfamily proteins (Table S3). In order to identify the most robust candidate genes for salt tolerance, all 176 genes were analyzed using "TM-1" RNA-seq expression profiles data at different time points of salt treatment 1, 3, 6, and 12 hin leaves part area (http://mascotton.njau.edu.cn). Out of 176 genes, 12 genes were highly expressed in salt treatment (see Table S3 and Figures 5 and 6). Of the 14 major QTLs identified, 12 genes were predicted to be the putative candidate genes associated with salt-tolerant traits.
In order to understand which QTL is preferable out of the 14 detected, we based our research on the additive effect and gene action (GA). Therefore, the parental contribution was determined by the use of additive effect. For positive and negative additive effects, alleles came from CCRI35 (the tolerant parent) and NH (the sensitive parent), respectively. GA was obtained from the ratio dominant effect by additive effect, |d/a| (see Table 3). To get good QTLs for salt tolerance, we focused our analysis on the QTLs with positive additive effect, coming from the tolerant parent, CCRI35. Therefore, seven QTLs were involved with the four following traits: FW, MDA, SL, and SLW (Figure 7). The genes' actions were partial dominance: PD (with FW trait), dominance: D (with FW trait), and overdominance: OD (with MDA, SLW, FW, and SL traits). The highest number of QTLs had an OD effect. This result suggests that these QTLs can be used for further studies in salt tolerance. Moreover, the 12 genes identified were classified into three groups, with two, two, and eight genes in Group1, Group2, and Group3, respectively. The best expression pattern was noted in Group1 (see Figure 6 and Table S3). In all the clusters, 5596 genes were mined, of which 176 were found to be related to biotic and abiotic stress. The stress-related genes were found to belong to the NAC, WRKY, HAD, C2H2, and RING/U-box superfamily proteins (Table S3). In order to identify the most robust candidate genes for salt tolerance, all 176 genes were analyzed using "TM-1" RNA-seq expression profiles data at different time points of salt treatment 1, 3, 6, and 12 hin leaves part area (http://mascotton.njau.edu.cn). Out of 176 genes, 12 genes were highly expressed in salt treatment (see Table S3 and Figures 5 and 6). Of the 14 major QTLs identified, 12 genes were predicted to be the putative candidate genes associated with salt-tolerant traits.
In order to understand which QTL is preferable out of the 14 detected, we based our research on the additive effect and gene action (GA). Therefore, the parental contribution was determined by the use of additive effect. For positive and negative additive effects, alleles came from CCRI35 (the tolerant parent) and NH (the sensitive parent), respectively. GA was obtained from the ratio dominant effect by additive effect, |d/a| (see Table 3). To get good QTLs for salt tolerance, we focused our analysis on the QTLs with positive additive effect, coming from the tolerant parent, CCRI35. Therefore, seven QTLs were involved with the four following traits: FW, MDA, SL, and SLW ( Figure 7). The genes' actions were partial dominance: PD (with FW trait), dominance: D (with FW trait), and overdominance: OD (with MDA, SLW, FW, and SL traits). The highest number of QTLs had an OD effect. This result suggests that these QTLs can be used for further studies in salt tolerance. Moreover, the 12 genes identified were classified into three groups, with two, two, and eight genes in Group1, Group2, and Group3, respectively. The best expression pattern was noted in Group1 (see Figure 6 and Table S3).

Collinearity Analysis
The syntenic blocks were obtained by comparing the 5178 marker position in the genetic map and physical map using AD cotton genome. The results showed that most of the markers had good collinearity ( Figure 8). However, some chromosomes, such as A10 and A12, showed poor syntenic blocks. This could be attributed to the mutation effect.

Collinearity Analysis
The syntenic blocks were obtained by comparing the 5178 marker position in the genetic map and physical map using AD cotton genome. The results showed that most of the markers had good collinearity ( Figure 8). However, some chromosomes, such as A10 and A12, showed poor syntenic blocks. This could be attributed to the mutation effect.

Discussion
This study was carried out in a controlled environment and the F 2:3 population and their parents were exposed to salt at the seedling stage since much evidence points to the fact that salt tolerance is a stage-specific phenomenon and seedlings are the most vulnerable to salinity stress [32]. The main aim of this research was to determine the response of cotton seedlings to salt stress when grown in a salt environment. The seedlings were highly assimilated to the natural environment by using soil as the rooting media. Soil substrate type should be considered when characterizing plant growth and physiological responses to salinity [33].
The plant materials used in developing the mapping population had varying salt stress tolerance. Significant differences between the two parents for all the measured traits were observed, both under stress and control conditions. Under stress, all the traits exhibited a significant reduction in LFW, SLW, GR, SL, and MDA. In the sensitive cultivar (NH), MDA concentration was higher than in the tolerant cultivar (CCRI35). MDA is a biochemical component of the cell; its release is triggered by exposure to stress. The lower MDA level in the tolerant cultivar is an indication of the ability of the plant to tolerate salt stress; this result is consistent with previous publications [34,35]. Moreover, the low concentration of MDA in the tolerant cultivar can be attributed to either an internal mechanism to convert the released MDA into other non-toxic compounds or minimizing the release of MDA to a minimal threshold with no major effect on the plant cell. Previous reports have shown that MDA in salt-sensitive plants is more pronounced than in salt-tolerant ones [36,37]. Furthermore, similar results were found by [38], stipulating that with regard to MDA content reduction in tolerant cultivars, it has been shown that low MDA content inhibits membrane damage by ROS and hence confers tolerance [39,40].
Under controlled conditions, in all the measured traits, no significant statistical difference was noted between the tolerant parent (CCRI35) and the sensitive parent (NH), except for CHL concentration. The sensitive cultivar produced more chlorophyll than the tolerant one. The increased level of chlorophyll content in the sensitive cultivar can be explained by the ability of the plant to maximize its photosynthetic apparatus, preparing for any sudden alteration that may occur within its environment, thus enhancing its survival. This result was consistent with the findings of [38], which reported that the chlorophyll (a), chlorophyll (b), and total chlorophyll contents were significantly lower in tolerant cultivars in comparison to sensitive cultivars. Moreover, plants under salt stress conditions are often observed to have reduced chlorophyll content [40,41].
Under stress, five traits showed significant difference (SL, GR, SLW, MDA, and LFW), while the rest showed no significant difference. However, as for chlorophyll concentration (CHL), in all the genotypes, the tolerant cultivar showed relatively higher concentrations in chlorophyll but statistically, no significant difference. It has been reported that salt-tolerant plants exhibit increased or unchanged chlorophyll content under halophytic conditions, whereas chlorophyll concentrations decrease in salt-sensitive species, indicating that this parameter can be considered a biochemical marker of salt tolerance in plants [42,43]. The higher number of chloroplasts per leaf area unit may be the main reason for the increased chlorophyll concentration in tolerant plants under salt stress. However, in the sensitive cultivar, salt stress may lead to stomatal closure, which in turn reduces the availability of carbon (IV) oxide in the leaves and therefore inhibits carbon fixation, exposing chloroplasts to excessive excitation energy, which in turn could increase the generation of reactive oxygen species (for example, MDA or other toxic elements) and induce oxidative stress [44]. High accumulation of sodium in plant tissues has been reported as one of the factors contributing to the reduction of photosynthetic pigments and rate of photosynthesis [45,46]. Chlorophyll is an important factor in plant photosynthesis; long exposure to salt stress may lead to the impairing of the photosynthetic apparatus of the plant, and in turn lead to a reduction in the plant's performance in terms of yield and the quality of the salt-sensitive cultivars.
Most of the traits were positively correlated. MDA was negatively correlated to all traits in 150 mM of salt treatment but no significant correlation was noted in 110 mM of salt treatment. Plant growth and development are impaired with an increase in salt concentration [47]. Previous studies have shown that 150 mM NaCl results in significant growth reduction between the salt-tolerant and less salt-tolerant cotton cultivars [48,49]. A negative correlation was observed between EC and the following traits; GR, DLW, and DLW with RWC at 150 mM salt concentration. The results obtained are in agreement with previous findings in Brassica napus under salt stress: EC was significantly negatively correlated with RL, RDW, and SFW [50].
The extent of transmission of traits from the parents to the offspring is determined by levels of heritability and hence traits with higher heritability percentage are easier to manipulate [51].  RWC (42.3). These results were in agreement with previous reports that indicated that heritability estimates in cotton are moderate to high for SH, RL, RDW, SFW, RFW, and SDW under salt conditions [52]. Selection, done based on low heritability, may be biased towards environmental influence on the genetic make-up of the plant and thus is not an effective way to improve plants [53]. In this study, the number of QTLs identified for high heritability traits was consistent, except for LFW.
Transgression was also a phenomenon noted in plant species. Ten QTLs were contributed by the sensitive parent (NH), while only four were contributed by the resistant parent (CCRI35). The positive and negative signs of the additive effect at the different loci showed the parental lines' contribution and confirmed the transgressive segregation pattern observed at the phenotypic level. Transgressive segregation for morphological and agronomical traits always follows the Mendelian pattern [54]. This observation was also reported by [55,56] in RILs population of Arabidopsis thaliana cultivated in a nitrogen environment.
In this study, four types of gene actions (GA) were revealed by the genetic effects: four genes exhibited dominant effects (D), 13 partial dominance (PD), 28 overdominance (OD) and two had an additive effect (A). Overdominance (OD) was observed for most of the traits in response to salt tolerance. This result was in disagreement with Oluoch Q, (2016) [25], who stated that partial dominance (PD) was observed for most traits in response to salt tolerance. In Oluoch et al. (2006), the marker and map size were relatively smaller than in our map; their final linkage map consisted of 1295 markers that amplified 1342 loci [25], and this could hamper the effective determination of the gene action identification. In this study, we developed a high-density linkage map with ≤1 cM of marker intervals, which enabled us to identify QTLs with more accuracy and higher resolution than earlier reports.
A total of 14 consistent QTLs for six traits were detected, which explained the phenotypic variance ranging from 2.72% to 9.87%. Five QTLs were located in the A t sub-genome, while the remaining nine were located in the D t sub-genome. This finding was consistent with earlier reports [25,57] that suggested that the D t sub-genome harbored genes or QTLs for stress biotic, abiotic, and fiber quality. Oluoch et al. (2016) stipulated that, of the 11 significant QTLs detected, only one (qSdw-Chr9-1) on chromosome 9 was located in the A t sub-genome, while the remaining 10 were located in the D t sub-genome. Meanwhile, approximately 58 QTLs were identified on the A t sub-genome chromosomes, whereas 107 QTLs were localized on the D t sub-genome chromosomes [57]. Saeed et al. (2014) highlighted the contribution of the D t sub-genome of tetraploid cotton to abiotic stress tolerance [14]. QTLs were also present at unique loci with defined underlying markers on different chromosomes, and none were shared by two or more traits. However, some were present in close proximity to each other, such as QTLs on chromosome D03 with the following traits: MDA and FW.
Of the 12 key genes found in this study, 11 genes were located in the D t sub-genome. Only Gh_A12G0454 gene on chromosome A12, which belongs to the zinc finger (C2H2 type) family protein, was located in the A t sub-genome. This result could explain the high number of QTLs in the D t sub-genome compared to the A t sub-genome in upland cotton. This finding is consistent with reports by Oluoch et al. (2016) and Jamshed et al. (2016) [25,57]. Moreover, the 12 key genes were localized in six families: RING/U-box superfamily protein, NAC, WRKY, Haloacid dehalogenase-like hydrolase (HAD) superfamily protein, zinc finger (C2H2 type) family protein, and MYB. Gh_D03G0935, Gh_D06G0163, and Gh_D06G0462 genes belonged to RING/U-box superfamily protein and were located on Clusters 5 and 6, respectively. The U-box genes in Arabidopsis such as AtPUB23 have been reported to show strong upregulation in the roots under salt and drought conditions [58]. The highest expression was observed with Gh_D03G0654 and Gh_D06G0417 genes in D03 and D06, respectively ( Figure 6 and Table S3). These two genes belonged to Group1 and could be helpful in further salt-tolerant gene studies. Furthermore, C 2 H 2 (94 genes), MYB (225 genes), NAC (170 genes), and WRKY (222 genes) unigenes were mostly upregulated under salt stress. MYB, NAC, and WRKY were also highly enriched at 4 h and 24 h of salt stress [15].

Greenhouse Experiment
Seeds of two upland cotton (G. hirsutum) lines, CCRI35 (salt-tolerant) and Nan Dan Ba Di Da Hua or NH (salt-sensitive) were used. The salt-tolerant cultivar was used as the female plant while the salt-sensitive one was the male plant. All plant materials were obtained from the National Mid-Term Gene Bank of the Institute of Cotton Research, Chinese Academy of Agricultural Sciences [59]. The seeds were then sterilized in 70% ethanol for 15 s, and put in 4% sodium hypochlorite for 15 min. The seeds were later submerged in sterile water for 12 h to rinse off the treatment chemicals. The seeds were sown in sterile silica sand with 0, 110, and 150 mM salt concentrations in planting boxes measuring 60 cm by 35 cm and depths of 12 cm in a greenhouse, where the conditions were optimized to 28/22 • C day/night, 60-80% relative humidity, and a 14-h photoperiod under 450 µmol·m −2 ·s light intensity. The experiment was done by laying out a complete randomized block design (CRD) with three replicates, with 10 and four seeds sown in 60 cm long and 35 cm wide rows, respectively. Two experimental periods were chosen, the first one from April to June 2015, and the second one from April to June 2016, corresponding to the cotton cultivation seasons in the Yellow River region of China with different salt concentrations of 110 mM or 150 mM. In both experiments, the plant population size was 277, with three replications. The two salt concentrations, 110 mM and 150 mM, were chosen based on the proportion of Hoagland [60]. Seven days after sowing (D.A.S), five seedlings with uniform growth for each treatment were sampled for all traits. Seed germination (GR) was determined by counting the seedlings in each replication. Five fresh seedlings were weighed for fresh weight (FW) in each replicate; leaf fresh weight (LFW) and stem length (SL) were measured immediately after collection. The sampled leaves were submerged in distilled water overnight and their saturated leaves were weighed (SLW). The leaf samples were then oven dried at 60 • C for 48 h; after attaining constant mass, all were weighed to obtain the dry leaf weight (DLW). Relative water content (RWC) was calculated using the formula described by [61]. The electric conductivity (EC) of the sample leaves was determined by using a sliced leaf portion of 0.5 g dipped in 100 mL of ddH 2 O; initial measurements were taken at the start, while the final measurements were taken after 48 h, as explained by [62]. For chlorophyll content (CHL) estimation, we used a chlorophyll meter (SPAD 502 m, Minolta, Osaka, Japan), and the average of five measurements on the same leaf was taken.

DNA Quantification and Qualification
The DNA of the entire population of 277 individuals of the F 2:3 generation, together with 10 samples for each parent, was extracted by the CTAB method [63]. Fresh and young leaves were collected and immediately frozen in liquid nitrogen. Each sample was then crushed in liquid nitrogen into a fine powder, then immediately added to a CTAB solution. For each 100 mg of homogenized tissue, we used 500 µL of CTAB extraction buffer. The contents were mixed and thoroughly vortexed.
The homogenized mixture was then incubated in a 60 • C water bath for 30 min. Following the incubation period, we centrifuged the homogenate for 5 min at 12,000× g. After centrifuging, the supernatant was transferred to a new tube. Then 5 µL of RNase solution was added to remove RNA and the result was incubated at 32 • C for 20 min. An equal volume of chloroform/isoamyl alcohol (24:1) was added and the sample was vortexed for 5 sand centrifuged for 1 min at 12,000× g to separate the phases. We transferred the aqueous upper phase to a new tube; the method was then repeated until the upper phase was clear. The upper aqueous phase was then transferred to a new tube. DNA was precipitated by adding 70% volume pre-refrigerated isopropanol and incubating at −20 • C for 15 min. The precipitated DNA samples were then centrifuged at 12,000× g for 10 min. The supernatant was then decanted without disturbing the pellet and subsequently washed with 500 µL ice cold 70% ethanol twice, then absolute alcohol. Then DNA pellets were later dissolved in 20 µL TE buffer (10 mM Tris, pH 8, 1 mM EDTA) [64]. DNA degradation and contamination were monitored on 1% agarose gels. DNA purity was checked using the Nano Photometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA). The ratio of absorbance at 260 nm and 280 nm was used to assess the purity of DNA. The DNA samples with a ratio of~1.8 were selected as pure [65]. The DNA concentration was measured using Qubit ® DNA Assay Kit in Qubit ® 2.0 Fluorimeter (Life Technologies, Carlsbad, CA, USA). The Qubit ® dsDNA HS (High Sensitivity) Assay Kits make DNA quantitation easy and accurate. The kits include concentrated assay reagent, dilution buffer, and pre-diluted DNA standards. The reagents were diluted by the buffer solution, then we added 1-20 µL of each DNA sample. The concentrations were read using the Qubit ® Fluorometer; only DNA samples with a concentration range of 10 pg/µL to 100 ng/µL were used (https://tools.thermofisher.com/content/sfs/manuals/Qubit_dsDNA_HS_Assay_UG.pdf).

GBS Library Preparation, Sequencing, and SNP Genotyping
Genotyping-by-sequencing (GBS) is an efficient method of large-scale genotyping, which is based on reduced-representation libraries (RRL) and high-throughput sequencing. First, we performed a GBS pre-design experiment to confirm the effectiveness of the GBS protocol and the quality of the output data. The enzymes and sizes of restriction fragments were evaluated using training data. Three criteria were considered: (i) The number of tags must be suitable for the specific needs of the research project; (ii) the enzymatic tags must be evenly distributed through the sequences to be examined; (iii) repeated tags must be avoided. These considerations improved the efficiency of GBS. To maintain the sequence depth uniformity of different fragments, a narrow length range was selected (about 50 bp).
Next, we constructed the GBS library in accordance with the pre-designed scheme. For the F 2:3 populations, genomic DNA was incubated at 37 • C with MseI (New England Biolabs, NEB, Ipswich, MA, USA), T 4 DNA ligase [66], ATP [66], and MseI Y adapter N containing barcode. Restriction-ligation reactions were heat-inactivated at 65 • C and then digested for additional restriction enzyme NlaIII at 37 • C. The restriction ligation samples were purified with Agencourt AMPure XP (Beckman, Brea, CA, USA). Then a PCR reaction was performed using purified samples, Phusion Master Mix [66] universal primer, and index primer to add index, complete i5 and i7 sequence. The PCR productions were purified using Agencourt AMPure XP (Beckman) and pooled, then run out on a 2% agarose gel. Fragments with 375-400 bp (with indexes and adaptors) in size were isolated using a Gel Extraction Kit (Qiagen, Hilden, Germany). These fragment products were then purified using Agencourt AMPure XP (Beckman), which was diluted for sequencing.
Genotyping-by-sequencing (GBS) was carried out as outlined by Elshire et al. (2011) [67], integrating three 96-well plates across 288 barcodes for library preparation and sequencing. For SNP calling, the raw sequence data for the 277 F2 progeny plus the F1 progenitor were processed through the TASSEL 3.0 GBS pipeline [68] using Gossypium_hirsutum_v1.1.fa as the reference genome [69], obtained from Nanjing Agriculture University's Cotton Research Institute (http://mascotton.njau. edu.cn/info/1054/1118.htm) for alignment and the Burrows-Wheeler Aligner (BWA) mem [70] with default parameters. The output consisted of variant call format (VCF) file version 4.1 [71], including SNPs present in at least 40% of the progeny and with a minor allele frequency (MAF) 0.1. Subsequently, the VCF was filtered using vcf tools version1.12a [71] and TASSEL [72] versions 3.0 and 4.0. A total of 93,384 SNPs were identified in 277 F 2:3 progeny by TASSEL 3.0, then a custom filtering process was applied for alignment. The filtering was based on keeping sites with a minimum read depth of 6% and 75% completeness by site across progeny and by progeny across sites. The results were produced as a TASSEL hapmap file.
Finally, using custom perl script markers heterozygous in the F1 progenitor and witha co-dominant 1:2:1 segregation among the F 2:3 , progenies were identified by a chi-squared (χ 2 ) goodness-of-fit test at α ≤ 0.01. These were reformatted to be imported in JoinMap ® 4.1 [73] for linkage group determination. A total of 26 linkage groups were formed; each linkage group was matched to its corresponding chromosome using BLASTN for the marker sequence (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Construction of the Linkage Maps and Data Analysis
ANOVA was performed using the greenhouse experiment data of the two seasons, 2015 and 2016. A mixed procedure was used; the genotypes and the environments were fixed as factors in order to detect the heritability. A post hoc test (Turkey's) was used to compare means. The broad-sense heritability (H) was calculated using the formula described by [31]: H = σ 2 G/σ 2 G + (σ 2 e/r) where σ 2 G is the genotypic variance; σ 2 e: phenotypic variance and r: replication. All the data were analyzed using R software version 3.4.1 [74].
Linkage map analysis was conducted using Join Map 4.0 [73] with a recombination frequency of 0.40 and a LOD score of 2.5 for the F 2:3 population. The Kosambi mapping function was used to convert the recombination frequencies to map distances. For the two experiments, each data point represented the mean of three replications, with each treatment consisting of 40 plants. Salt-tolerance-related traits EC, GR, SL, FW, LFW, SLW, DLW, RWC, and CHL were used to conduct a QTL analysis. WinQTL Cartographer 2.5 was used to detect QTLs using composite interval mapping by [75]. In the composite interval mapping (CIM) method, model 6, forward-backward regression method with 1 cM walking speed, a probability into and out of the model of 0.01, and window size set at 10 cM were used. A stringent logarithm of odds (LOD) threshold value was estimated by 1000 permutation test for all traits and used to determine the significant QTLs with a significance level of p = 0.05. However, QTLs in varying environments with an LOD threshold of at least 2.5 were considered common QTLs based on the explanation by Lander and Kruglyak [76].
QTL nomenclature was done based on previous criteria [77]. The phenotypic variance observed proportions as illustrated by each QTL was estimated by coefficient of determination R 2 (%) as a percentage. The additive and dominance effects from QTL cartographer results were used to calculate the genetic effects (|d/a|). The results were used to classify the QTL as additive (A) (0-0.20), partially dominant (PD) (0.21-0.80), dominant (D) (0.81-1.20), or over-dominant (OD) >1.20, according to [78]. The graphic presentation of the linkage group and QTLs marked was created by Map Chart 2.2 [79] and the R/qtl package was used to generate the genetic map in R software version 3.4.1 [74].
In this study, only consistent QTLs were used to identify the crucial candidate genes for salinity tolerance. The genes were identified through the available resources [80] (http://mascotton.njau.edu.cn). The physical position of the SNP markers flanking major QTLs for salinity tolerance was used to find the genes located in each cluster. The function of the identified genes was determined through gene annotation. Furthermore, the expression profiles of the candidate genes were analyzed by mapping the genes in the "TM-1" RNA-seq transcriptome data of cotton leaf tissues at different time points of salt treatment (1, 2, 3, 6, and 12 h) (http://mascotton.njau.edu.cn). The putative candidate genes were mapped using VLOOK UP Excel formula and the output FPKM for each gene was expressed to log10 to get the final gene expression. These values were used to generate the heat map by using R script.

Collinearity Analysis
Based on the high-density genetic linkage map and sequences corresponding to markers, collinearity analysis between genetic and physical maps within chromosome was carried out using a BLASTN (https://blast.ncbi.nlm.nih.gov/Blast.cgi) search with E ≤ 1 × 10 5 , identity ≥80%, and matched length ≥200 bp. Next, the best hit for each marker was chosen and all the best hits were illustrated intuitively using online drawing tools (http://circos.ca/).

Conclusions
We identified 66 QTLs, of which 14 QTLs were consistent in eight chromosomes in at least two environments. A total of 5596 genes were identified; 12 genes showed preferential expression in leaf tissues at different time points (1, 3, 6, and 12 h) of salt treatment in the "TM-1" transcriptome data. The genes found in this research work were in three groups ( Figure 6); 11 genes were from the D t sub-genome, and only one (Gh_A12G0454 in chromosome A12) was identified from the A t sub-genome. Cloning and gene saturation need to be done to verify their specific functional role in cotton plants under salt stress.