Phenotypic variation of lint yield components under salt condition
At maturity stage, three traits including SBW, LP and BNPP were measured in 316 upland cotton accessions (Additional file 1: Table S1) with four salt conditions (Additional file 2: Fig. S1A and Additional file 3: Table S2). The average values of SBW, LP and BNPP traits ranged from 5.28g to 6.84g (Additional file 4: Table S3), 37.89% to 40.39% (Additional file 5: Table S4) and 5.01 to 7.74 (Additional file 6: Table S5) under four different salt conditions, respectively. The BNPP had the largest coefficients of variation (CV) ranging from 22.53% to 48.08%, and the LP had the smallest (8.15%~12.68%). For decreasing the environmental errors, we further evaluated the BLUP value of the three traits in the same salt condition. The BLUP value showed that the distribution of SBW was at 5~7 g (Additional file 2: Fig. S1B) and LP at 35~45% (Additional file 2: Fig. S1C). However, the BNPP showed larger difference with 3~7 in salt condition A and B, and 5~9 in salt condition C and D (Additional file 2: Fig. S1D).
Paired-samples t-test was conducted to further investigate the phenotypic variation. The SBW was slightly higher in condition A than in B, C and D conditions (Fig. 1A). However, no significant difference of LP was detected under all four salt conditions (Fig. 1B). Most significantly, the BNPP decreased with the increase of salt concentrations (Fig. 1C). In addition, the LYPP also showed the decrease with the increase of salt concentrations mainly due to change of BNPP (Fig. 1D). We compared respectively the differences of SBW, LP, BNPP and LYPP, at the highest with the lowest salt conditions, and the results showed that SBW increased by 5.29%, LP remained unchanged, BNPP decreased by 17.75%, and LYPP decreased by 16.79%. Taken together, Lint yield of cotton was mainly affected by boll number decreased under high salt concentration, although SBW slightly increase relevant to severe boll numbers decrease.
Correlation analysis was conducted among SBW, LP and BNPP under four salt conditions (Additional file 7: Fig. S2). The results showed that there was low correlation among the three traits. As for the same trait, under different salt environments, LP showed high correlation with R value from 0.75 to 0.83, followed by SBW with R value from 0.43 to 0.63. BNPP showed no or weak correlation in four salt conditions, further indicating that BNPP was largely affected by salt concentration.
GWAS of the lint yield related traits
With the genotypic data of 57,413 high-quality SNPs [17], the GWAS was conducted for the three traits with different environments and the BLUP values using six methods (“mrMLM”, “FASTmrMLM”, “FASTmrEMMA”, “pKWmEB”, “ISIS EM-BLASSO” and “pLARmEB”) of multi-loci MLM model in R package “mrMLM” [18]. In total, 854 QTNs on 26 chromosomes were identified as significantly associated with the three traits (Additional file 8: Table S6). Following the linkage disequilibrium (LD) value [17] and referring the method reported by Song et al. (2019) [19], the 200-kb upstream and downstream of QTNs were defined as one QTL. In total, 600 QTLs including 151 of SBW, 417 of LP and 112 of BNPP, were detected with 1446 times under different environments and BLUPs by six multi-loci MLM models (Table 1). For each trait, most QTLs (85 of SBW, 235 of LP and 65 of BNPP) were detected only once, implying these QTLs be apt to be affected by environmental condition (Additional file 9: Fig. S3). To improve the reliability and stability of associated QTLs, we selected QTLs detected three or more times across different methods or environments as stable QTLs for further analysis. As a result, 42, 91 and 25 QTLs were identified in SBW, LP and BNPP, respectively (Table 1).
The chromosomal distribution showed that these stable QTLs were widely distributed on 26 chromosomes, with more QTLs of SBW and BNPP located on At sub-genome than on Dt sub-genome, while QTLs of LP showed opposite (Fig. 2A). Most of SBW located on chromosomes A11 and A12, LP on A08, D06 and D13, and BNPP on A05, A12 and D07 (Fig. 2B). The vein diagram of these stable QTLs showed that no overlapped QTL was detected within three traits and most of QTLs were specific for individual trait (Fig. 3A), implying great differences in regulating pathway for the three traits. We further analyzed the QTLs of each trait under different salt conditions, only one QTL of SBW was detected under all four salt conditions, and no overlapped QTL detected in BNPP (Fig. 3B, C). However, most QTLs of LP were overlapped under different salt conditions, there were 8 QTLs of LP detected under all four salt conditions (Fig. 3D). The results suggest complex and various regulatory mechanism involved in SBW and BNPP under different salt conditions but stable and high heritability of LP against salt stress.
Identification of candidate genes in QTLs
Potential candidate genes in these stable QTL regions were extracted based on the released G. hirsutum TM-1 genome [20]. In total, 1166, 2748, and 711 candidate genes were identified in the QTL regions for SBW, LP and BNPP, respectively (Fig. 2C), with most genes distributed on chromosome A12, D13 and A12 for SBW, LP and BNPP, respectively (Fig. 2D). With GO analysis, the genes in QTL regions for SBW enriched in “embryo development” and “regulation of cell shape” (Additional file 10: Fig. S4 and Additional file 11: Table S7). The genes from LP QTLs mainly enriched in several pathways, including “regulation of organ growth”, hormone and ROS regulation such as “regulation of gibberellic acid mediated signaling pathway”, “positive regulation of reactive oxygen species metabolic process”, “brassinosteroid biosynthetic process” and “defense response by cell wall thickening”, and carbohydrate metabolism such as “glycosylation”, “glucose metabolic process”, “monosaccharide biosynthetic process” and “hexose biosynthetic process”, which was consistent with the previous reports that these Go items played the crucial roles in fiber development [11, 21]. In addition, we identified 14, 21 and 10 genes related to “Golgi vesicle transport”, “plant-type secondary cell wall biogenesis” and “glucose metabolic process”, respectively, however, none of these process-related genes were found in the QTL regions of BNPP and SBW (Additional file 12: Fig. S5 and Additional file 11: Table S7). The function of genes associated with BNPP mainly enriched in “mitotic cell cycle”, “ion transmembrane transport” and “polysaccharide catabolic process” (Additional file 13: Fig. S6 and Additional file 11: Table S7), implying that BNPP was closely related to stress response.
Genes relevant to LP
Via tissue and organ transcriptome profile [22], we identified 182 genes from LP QTLs with predominant expression during fiber development. Of them, 29, 35, 73 and 45 genes highly expressed at 10 DPA, 15 DPA, 20 DPA and 25 DPA, respectively (Additional file 14: Fig. S7 and Additional file 15: Table S8). Further, we focused on the genes in the regions of 8 overlapped QTLs under all four salt conditions and identified 10 genes predominantly expressed during fiber development (Table S9). Of them, Gh_A05G2488 encoding auxin transport facilitator family called PIN-FORMED LIKES proteins (PILS), located in a high frequency associated QTL (A05: 32377816-33100112) which was detected 21 times with multi-methods and environments. Auxin is essential for plant growth and development. Prominent auxin carriers with fundamental importance during plant development are PIN-FORMED (PIN) proteins [23, 24]. In Arabidopsis, overexpressing PILS genes could reduce hypocotyl and root growth [25]. However, auxin accumulation can promote cell initiation (-2 to 2 DPA) in the fiber cell. Transgenic assay showed that ovule-specific suppression of multiple GhPIN genes inhibited both fiber initiation and elongation in cotton [26], indicating that Gh_A05G2488 might play an important role in the fiber development. Two genes, Gh_D13G0342 encoding RAB GTPase homolog G3F (RABG3F) and Gh_A10G2138 encoding prenylated RAB acceptor protein 1 (PRA1), located on QTL (D13: 3297533-3912736) and QTL (A10: 99896949-100396471), respectively. Previous reports suggested that RAB genes played crucial roles in cell polarity growth, including elongation of pollen tubes [27] and root hairs [28]. In addition, RAB genes also involved in cotton fiber development [29]. Prenylated Rab acceptor 1 (PRA1) domain proteins are small transmembrane proteins that regulate vesicle trafficking as receptors of Rab GTPases. AtPRA1 proteins were localized in the endoplasmic reticulum, Golgi apparatus, and endosomes/prevacuolar compartments, indicating a function in both secretory and endocytic intracellular trafficking pathways [30]. These results indicates that the RAB and RAB acceptor protein have potential important functions in cotton fiber development.
Cotton seed transfer cells are enriched in callose, which regulated fiber elongation and secondary wall thickening [31]. Zhang et al. (2007) showed that the transcription factor MYB103 affects callose dissolution during the anther development in Arabidopsis [32]. We found Gh_D03G1419 encoding MYB103 transcription factor (named as GhMYB103), located in a high frequency associated QTL (D03: 42299450-43529568). Sequence analysis showed that two QTNs associated with LP located in the 3739 bp upstream (TM55216, D03: 42969311) and exon (TM55217, D03: 42973276) regions of the gene, respectively (Fig. 4A). The single nucleotide mutation (from C to G) at TM55217 locus led to the change of amino acid from leucine (L) to valine (V) (Fig. 4A). Through Student’s t test, we found that LP with A genotype in TM55216 was significantly higher than with G genotype (Fig. 4B), and with G genotype in TM55217 significantly higher than with C genotype (Fig. 4C). The two QTNs could generate 3 haplotypes including H1: AG, H2: AC and H3: GC. LP with AG and AC haplotypes were significantly higher than that with GC haplotypes. However, there was no significant difference between AG and AC haplotypes, implying that QTN TM55216 might play more important roles in LP (Fig. 4D).
In addition, we integrated the LP QTLs with GWAS signals published in previous reports [11, 12], and found three QTLs (A12: 602614-743324; D13: 58792627-59289811; A09: 4676815-5076815) overlapped with GWAS signals. The three QTLs were detected in different salt conditions. Further, 10 genes from the QTL regions were identified to be predominant expression during fiber development (Additional file 16: Table S9). Of them, ABP1 [33] and CPK17 [34] was related to hypocotyl growth or fiber development.
Taken together, the genes from these stable LP QTLs and expressed predominantly in fiber developmental stages play an important role for the LP improvement in breeding practice.
Genes relevant to BNPP
The QTLs of BNPP were easily affected in different environments. No overlapped QTL under all four salt conditions were detected and most QTLs were identified only in one or two salt environments. For example, QTL TM57617_TM57620 (D05: 24.3-24.8 Mb, detected with 16 times) and TM74225 (D10: 24.1-24.5 Mb, detected with 7 times) were identified in salt conditions C and D. However, QTL TM52041_TM52044 (D02: 49.2-49.7 Mb, detected with 6 times) was identified only in salt condition A, and QTL TM29006 (A08: 81.7-82.1 Mb, detected with 5 times) was identified only in salt condition B. These results indicated that the genes regulating the number of bolls were various in different salt conditions.
A large number of items related to ion transport were enriched with GO analysis (Additional file 13: Fig. S6 and Additional file 11: Table S7). In order to further explore the mechanism of stress tolerance on the increase of boll number, we performed RNA-seq analysis using salt stress transcriptome in Gossypium acc. TM-1. With the filter of FPKM ≥ 3, 500 genes (446 in roots and 395 in leaves) in QTL regions were obtained. With GO annotation and differential expression analysis, we further focused on 40 salt-inducible stress response genes, of them, 6 commonly identified in roots and leaves (Additional file 17: Fig. S8, Additional file 18: Table S10 and Additional file 19: Table S11). Gh_A04G1216 encoded a high-affinity K+ transporter 1 (HKT1). AtHKT1 limits the root-to-shoot sodium transportation and is believed to be essential for salt tolerance in Arabidopsis thaliana [35]. Gh_A05G3239 encoded a peroxidase superfamily protein (POD), which had been proved playing an important role in anti-oxidation under salt stress in cotton [36]. Gh_A08G1183 encoded a mitogen-activated protein kinase (MAPK), which was widely reported to be associated with salt tolerance in cotton [37]. Besides ion transport Go term, carbohydrate metabolism was active under salt stress, for example, 10 genes were identified in the enriched Go term “polysaccharide catabolic process” and 3 of them differentially expressed under salt stress. Cell cycle regulation is of pivotal importance for plant growth and development [38]. The 33 genes were identified in the enriched Go term “cell cycle” with 8 of them differentially expressed under salt stress. These candidate genes could contribute to increasing boll number under salt environment in cotton.
In order to explore the key genes and regulation mechanism of boll number under high and low salt conditions, we compared the genes located in QTL regions of BNPP between condition A and D. Totally, 204 and 265 genes were identified under salt condition A and D, respectively. GO enrichment analyses showed the significantly different function classification of genes between high and low salt conditions. Go terms related to “polysaccharide metabolic process”, “carbohydrate catabolic process” and “cell wall organization” were enriched under high salt condition A (Fig. 5A, Additional file 20: Table S12) and “cell cycle”, “ion transmembrane transport” and “regulation of signal transduction” under low salt condition D (Fig. 5B), which indicate that carbohydrate metabolism and cell structure maintenance play an crucial role under high salt condition and ion transport is more basal under low salt condition. Under high salt condition, several candidate genes associated with BNPP were detected. Gh_A11G1551 encodes a proline dehydrogenase 1 (ProDH1), also called early responsive to dehydration 5 (ERD5), which have been studied extensively, especially under abiotic stress [39]. To counteract osmotic stress caused by salt stress, some plants accumulate several kinds of compatible osmolytes, such as proline, glycine betaine, and sugar alcohols, to protect macromolecules and remains osmotic pressure equilibrium inside and outside cell membrane [40]. The expression of ProDH2, a highly homologous gene of ProDH1, can promote proline accumulation under stress conditions [41]. For energy metabolism, Gh_A05G1912 encodes an isoamylase 3 (ISA3) which contributes to starch breakdown. Atisa3 mutants have more leaf starch and a slower rate of starch breakdown than wild-type plants [42]. Under low salt condition, three candidate genes were identified to play important roles in the balance between sodium and potassium. In detail, Gh_A10G0441 encodes a potassium transporter 1 (KUP1) [43], Gh_A12G0074 encodes a high affinity K+ transporter 5 (HAK5) [44] and Gh_A12G0061 encodes a sodium hydrogen exchanger 2 (NHX2) [45].