Comprehensive identification and expression analysis of CRY gene family in Gossypium

Background The cryptochromes (CRY) are specific blue light receptors of plants and animals, which play crucial roles in physiological processes of plant growth, development, and stress tolerance. Results In the present work, a systematic analysis of the CRY gene family was performed on twelve cotton species, resulting in 18, 17, 17, 17, and 17 CRYs identified in five alloteraploid cottons (Gossypium hirsutum, G. barbadense, G. tomentosum, G. mustelinum and G. darwinii), respectively, and five to nine CRY genes in the seven diploid species. Phylogenetic analysis of protein-coding sequences revealed that CRY genes from cottons and Arabidopsis thaliana could be classified into seven clades. Synteny analysis suggested that the homoeolog of G. hirsutum Gh_A02G0384 has undergone an evolutionary loss event in the other four allotetraploid cotton species. Cis-element analysis predicated the possible functions of CRY genes in G. hirsutum. RNA-seq data revealed that Gh_D09G2225, Gh_A09G2012 and Gh_A11G1040 had high expressions in fiber cells of different developmental states. In addition, the expression levels of one (Gh_A03G0120), 15 and nine GhCRY genes were down-regulated following the PEG, NaCl and high-temperature treatments, respectively. For the low-temperature treatment, five GhCRY genes were induced, and five were repressed. These results indicated that most GhCRY genes negatively regulate the abiotic stress treatments. Conclusion We report the structures, domains, divergence, synteny, and cis-elements analyses systematically of G. hirsutum CRY genes. Possible biological functions of GhCRY genes in differential tissues as well as in response to abiotic stress during the cotton plant life cycle were predicted. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08440-9.

The numbers of cryptochromes vary among plant species, ranging from three in Arabidopsis to seven in soybean [8,9]. All higher plants studied to date have two phylogenetically diverged clades of cryptochromes, CRY1 and CRY2 [10,11]. Most cryptochromes in plants are involved in regulation of gene expressions of the plant life cycle [12][13][14]. Among them, CRY1 could inhibit the hypocotyl elongation in Arabidopsis [4], and the grain dormancy and germination in barley [15]. Arabidopsis CRY1 also controls photomorphogenesis through the regulation of H2A.Z deposition [16]. CRY2 could interact with CIB1 or SPA1 to regulate the floral initiation in Arabidopsis [17,18], and suppress the leaf senescence in soybean [19]. Furthermore, CRY1 and CRY2 act together to stimulate the stomata opening and development in Arabidopsis [20,21].
Among the many processes regulated by cryptochromes, responses to biotic and abiotic stresses, such as drought, salinity, heat, and so on, are one of the most active research topics in plant biology [22]. It has been demonstrated that cryptochromes in Arabidopsis enhance plant resistance to Pseudomonas syringae [23] and drought [21]. In tomato, CRY1a could modulate water deficit and osmotic stress responses [24] as well as mediate long-distance signaling of soil water deficit [25]. In rice, suppression of CRY1b improved salt tolerance as a result of down-regulation of the melatonin and brassinosteroid biosynthetic genes [26]. Overexpressing the wheat cryptochromes TaCRY1a and TaCRY2 in Arabidopsis led to higher sensitivity to salt stress in the transgenic plants [27]. In addition, the cry1 mutant of A. thaliana had a greater germination and seedling survival rate than the WT in salt-stressed conditions, and the mutant plants exhibited enhanced tolerance to salinity [28]. Obviously, these results support a role of cryptochromes acting as a negative regulator in plant response to salinity.
In the current study, we performed a genome-wide screening of CRY genes in cottons, based on data gathered from recent whole-genome sequencing results. We used in silico approach to identify CRY genes in Gossypium species, focusing gene structures, conserved domains, synteny, cis-elements, and the phylogenetic relationships. Moreover, the tissue-specific expression patterns and the transcriptional responses of GhCRYs to abiotic stresses were examined. Our data provide inspirations for further research of cotton CRYs, as well as for molecular design of cotton cultivars with desired traits.

Identification and chromosomal location of CRY family genes in G. hirsutum
Cryptochromes are a class of photolytic flavin proteins, which act as UV-A/blue light receptors and play an important role in plant growth and development [50]. These proteins are defined by the presence of a FADbinding domain of DNA photolyase [51,52]. Hmmersearch against the G. hirsutum genome database with the conserved domains (PF00875 for DNA photolyase domain) identified 18 CRY genes (Table 1), which are dispersed over 14 of the 26 G. hirsutum chromosomes, with most, but not all, homoeologs conserved in the two (A and D) subgenomes (Fig. 1).

Structural organization of GhCRY genes
Less than two-fold variation in length was detected in the predicted coding sequences (CDS) for the recovered GhCRY s, from 1380 bp for Gh_D06G1145 to 2,268 bp for Gh_A02G0384/Gh_D02G0436 (Table 1), which translate to proteins ranging from 459 amino acids (aa) (52.29 kDa) to 755 aa (85.81 kDa). Predicted isoelectric points (pI) for members of this family vary widely, from 5.55 to 9.43. All of the putative GhCRY proteins have DNA photolyase domain in the N-terminal region (Table 1). Twelve GhCRY proteins have FAD binding 7 domain, four proteins have cryptochrome C domain and two have hydrolase 4 domain in the C-terminal region, respectively (Table 1).
While all putative GhCRY genes contain introns (Fig. 2), they also exhibit considerable variations, in both length and number. In general, homoeologous GhCRY genes show highly similar intron patterns, however, among different homoeologous pairs the genes vary in both intron numbers (3 to 13) and lengths. One of the homoeologous gene pairs does exhibit divergence in structure, namely Gh_A06G0969 vs Gh_D06G1145, which contain 12 and 10 introns, respectively. Characterization of parental genes (both containing 12 introns in the diploids) for the homoeologs suggests that this structural variation was descendant divergence rather than inherited. In addition, phylogenetic relationship of the GhCRY gene family is not consistent with the intron/exon structures characterized (Fig. 2).

Phylogenetic analysis of CRY family genes in Gossypium
The general but incomplete conservation of CRY genes between the two subgenomes of G, hirsutum prompted us to ask whether the minimal loss and/or gain occurred before or after the marriage of the two diploid progenitors. We specifically assessed this using the protein-coding sequences of 62 cotton CRY genes (G. hirsutum, 18; G. barbadense, 17; G. raimondii, 9; G. arboreum, 9; and G. herbaceum, 9) with 3 Arabidopsis thaliana CRY genes for phylogenetic analysis (Fig. 3). Seven clades (I-VII) were robustly supported with each of the A. thaliana genes associated with clades I, II and VI, respectively, and the reamining four clades were composed of Gossypium CRY genes only.
Overall, the expected diploid-polyploid topology is reflected in the tree for each set of orthologous/homoeologous genes, indicating general preservation during diploid divergence and through polyploid evolution. That is, the number of CRY genes in tetraploids was generally additive with respect to the model diploid progenitors, with each homoeolog (A t or D t ) sister to their respective counterparts in the diploid species. Clades I and II had the most CRY genes, and the other five clades contained equal numbers (Fig. 3). In clades I and II, genes related to AtCRY1 and AtCRY2 exhibit duplications in Gossypium species, which indicate a duplication event in Gossypium compared to A. thaliana. In addition, the Gossypium CRY genes of clade I have a sister relationship with AtCRY1, and clade II have the closest relationship with AtCRY2, and clade VI was classified with AtCRY3. Therefore, it is speculated that the function of CRY genes in these clades of cotton is similar to their homologs in Arabidopsis.
Although the CRY family exhibits general conservation, a few deviations were noted. For example, Clade II exhibits evidence of homoeolog loss; that is, the A t copy of GB_D02G0441 is missing from G. barbadense genome, whereas both copies (Gh_A02G0384/Gh_ D02G0436) exist in G. hirsutum. This gene loss might specific to G. barbadense after divergence of the two allotetraploid species.

Divergence of CRY genes in allotetraploid G. hirsutum and its diploid progenitors
The CRY genes in the two diploid species were then compared with G. hirsutum A t -and D t -subgenome homoeologs ( Fig. 4, Additional file 1: Table S1). To explore the evolutionary relationship and possible functional divergence of CRY genes between the allotetraploid cotton and their extent putative diploid progenitors, the nonsynonymous substitution (Ka) and synonymous substitution values (Ks) and the Ka/Ks ratios for each pair of the genes were calculated (Additional file 1: Table S1). By comparing the Ka and Ks values of 18 orthologous gene sets between the allotetraploid and the respective diploid genomes, we found that the Ka and Ks values are higher in the D t subgenome than in the A t subgenome (Fig. 4a, b). These results indicate that GhCRY genes in the D t subgenome tend to have experienced faster divergence than their A t counterparts. However, the Ka/Ks ratios of D t subgenome was lower than that of A t subgenome (Fig. 4c), indicating that GhCRY genes in A t subgenome were subjected to positive selection during the course of

Dynamic evolution of CRY family genes in Gossypium
We further evaluated the general preservation of CRY genes in 12 Gossypium species, Gossypioides kirkii and Arabidopsis thaliana (Fig. 5). In A. thaliana only three CRY genes were identified, whereas and the relative of Gossypium, Gossypioides kirkii, has eight. All Gossypium species surveyed recovered a minimum of five putative CRY genes in G. australe (G 2 ) to 18 in G. hirsutum (AD 1 ). Among the D genome species, CRY gene copy number varied from a minimum of 7 in G. thurberi (D 1 ), to 9 in both G. raimondii (D 5 ) and G. turneri (D 10 ). The two cultivated A-genome species of G. herbaceum (A 1 ) and G. arboreum (A 2 ) also have 9 copies. The sister-species of A-genome, G. longicalyx (F 1 ), contains 6 CRY genes. The . Notably, this high copy number in tetraploid is slightly more than double the copy number in diploid, likely reflective of the duplicated history of cotton. Comparatively, in G. hirsutum (AD 1 ) the CRY copy number is stable after polyploidization, whereas the other four allotetraploid cottons included in this analysis all appear to have undergone a homoeolog loss (17 versus18), while G. hirsutum has retained it on the A02 chromosome (Gh_A02G0384) (as stated above, Fig. 3).

Chromosomal distribution and synteny analysis of Gossypium CRY genes
Based on these Gossypium genomes, the location of CRY genes and the length of chromosomes from two diploid species and five allotetraploid species were used to analyze the chromosomal distribution and synteny (Fig. 6). High similarity was found in the chromosomal distribution patterns among these seven cotton species. The CRY genes were unevenly distributed on chromosomes with divergence detected between the diploid and allotetraploid species. For instance, no CRYs were found on Chr 02 of two A-genome species nor on Chr A02 of the allotetraploid specie of G. barbadense (AD 2 ), G. tomentosum (AD 3 ), G. mustelinum (AD 4 ) and G. darwinii (AD 5 ), except G. hirsutum (AD 1 ) which has one (Gh_A02G0384) on this chromosome (Fig. 6A). There were totally 104 CRYs distributed throughout the 80 chromosomes, comprising 38 located on the A or A t (sub)genomes and 42 located on the D or Dt (sub)genomes. The majority of CRYs were located on the proximate or the distal ends of the chromosomes. In addition, there were nine collinear gene pairs between G. raimondii and G. arboreum, nine between A t and D t subgenomes of G. hirsutum, and eight for the other four allotetraploid subgenomes.

Cis-element analysis of CRY genes in G. hirsutum
Cis-elements are involved in the responsive to corresponding stimulations to regulate the expression of genes [53]. In this study, a 1.5-kb fragment upstream to the start codon of each CRY gene of G. hirsutum was extracted to investigate putative cis-elements in the mediation of gene expression using the PlantCARE server [54]. In total 581 cis-elements among 18 GhCRY genes were identified, ranging from 22 in Gh_A02G0384 to 47 in Gh_D03G1520 (Fig. S1). Some cis-elements were predicted to mediate the phytohormone (ABRE) and stress (TC-rich repeats) responses (Fig. S1). 13 GhCRY gene promoters possess at least one abscisic acid responsiveness element (ABRE), and 11 GhCRY genes have at least one antioxidant response element (ARE). There were 9 GhCRYs which harbor the elements involved in the MeJA-responsiveness (CGTCA-motif and TGACG-motif ), four have the cisacting element involved in salicylic acid responsiveness (TCA element), two have gibberellin-responsive element (GARE-motif ), and three have auxin-responsive element (TGA-element and AuxRR-core). Of the 18 G. hirsutum CRY genes (GhCRY s), all have light responsiveness element (Box 4 and G box) and at least two MYB binding site elements (MYB), three have meristem expression element (CAT-box) and four (Gh_D02G0436, Gh_D03G1520, Gh_ A05G2282 and Gh_D05G2543) have low-temperature responsiveness element (LTR). The remaining elements related to stress, such as thoese of defense and stress responsiveness (TC-rich repeats), W-box recognized by WRKY transcription factors, wound-responsive element (WUN-motif ), and MYB-binding site involved in drought inducibility (MBS), were also detected.

Expression patterns of GhCRY genes in different G. hirsutum tissues
The expression profile of a gene family can provide valuable clues to the possible functions of the gene. Analysis of 18 GhCRY s showed that most genes differ in spatial expression patterns. For instance, the expression levels of Gh_A09G2012, Gh_D09G2225, Gh_A11G1040 and Gh_D11G1195 in root, stem, leaf, torus, stamen, pistil and calycle were significantly higher than those of other GhCRY genes (Fig. 7a). Gh_A06G1059 presented the highest expression level in petal (Fig. 7a). In addition, four genes (Gh_A09G2012, Gh_D09G2225, Gh_A11G1040 and Gh_D11G1195) also showed high expression in seed, root and cotyledon samples at different time points post seed germination (Fig. 7b). In ovule samples of different developmental stages, Gh_A09G2012 and Gh_D09G2225 had the highest expression levels followed by Gh_A11G1040 and Gh_D11G1195 (Fig. 7c). In fiber samples of different developmental stages, Gh_A09G2012 had the highest expression at 20-and 25-dpa (Fig. 7d), suggesting that this CRY gene might play a role at the cell wall thickening stage. Gh_D09G2225 was preferentially expressed in 0-dpa ovule. These two GhCRYs are homologous to Arabidopsis AtCRY2 (Fig. 3). Gh_A11G1040 showed the highest expression in 5-and 10-dpa fiber cells, which suggests that it may play a role at the elongation stage of fiber development (Fig. 7d). These three genes could be taken as candidate genes for subsequent transformation experiments to dissect their functions in cotton fiber development.

Expression changes of GhCRY genes in G. hirsutum under different stresses
Cotton is often subjected to a variety of abiotic stresses during its growth and development. Therefore, we analyzed the expression changes of CRY genes under simulated drought (PEG 6000), salt (NaCl), heat and cold abiotic stresses from RNA-seq data (Fig. 8). At different time points of PEG6000 simulation drought condition, expressions of most GhCRY genes were not changed (|lg 2 (Fold change)|≥ 1| as the threshold of differentially expressed genes). For instance, the expression of GhCRY genes did not change after PEG treatment for 3 and 6 h (Fig. 8a). However, Gh_A12G2401 was down-regulated after 1 h of PEG treatment (Fig. 8a). After 12 h of PEG treatment, the expression of one GhCRY gene (Gh_A03G0120) was repressed (Fig. 8a). These results indicated that only a limited nimber of GhCRY genes responded to drought stress. Under salt stress, the expressions of GhCRY gene did not change after 1, 3 and 6 h of NaCl treatment (Fig. 8b). However, a prolonged (12 h) NaCl treatment down-regulated 15 GhCRY genes (Fig. 8b), indicating that they were negatively related to salt stress.
At the four time points of high temperature stress, the expression of GhCRY gene did not change after 1 and 3 h of high temperature treatment (Fig. 8c). After 6 h under high temperature, the expression levels of four GhCRY genes (Gh_A05G2282, Gh_D06G1145, Gh_A12G2401 and Gh_D12G2528) decreased (Fig. 8c), and after 12 h of high temperature the down-regulaged GhCRY genes increased to nine (Fig. 8c). Notably, Gh_D06G1145 gene was inhibited at both 6 and 12 h, suggesting that this gene might be a key factor in adjusting the response to high temperature stress in G. hirsutum.
As for the four time points of low temperature stress, there were four GhCRY genes which were up-regulated after 1 h of low temperature treatment (Fig. 8d). And 3 h in low temperature the expressions of five GhCRY genes were induced, and two were reduced (Fig. 8d). After 6 h in low temperature the expression levels of eight GhCRY genes were elevated and two were declined (Fig. 8d). Finally, extending the low temperature treatment to 12 h resulted in five GhCRY s were up-regulated and five down-regulated (Fig. 8d). Among them Gh_A06G0969 and Gh_D06G1145, both homologous to AtCRY3, were induced at all the four-time points under the low temperature condition, suggesting their involvement in plant toleranbce to low temperature stress.
To confirm the accuracy of these GhCRY candidates in response to these stresses, the expressions of four GhCRY genes (Gh_A05G1941, Gh_A05G2282, Gh_A06G1059 and Gh_A12G2401) under hot, cold, salt and PEG stresses were examined by qPCR (Additional file 3: Figure  S2). These four genes were mostly repressed under different stress conditions, which verified their expression patterns under stress detected by RNA-seq and qPCR data were generally congruent.

Discussion
In plants the function of CRYs vary not only among individual CRY members but also among species. Arabidopsis thaliana has three CRYs, CRY1 and CRY2 are located in the nucleus and play a multifaceted role in various aspects of plant growth and development [1,55]. For instance, CRY1 primarily regulates photomorphogenic responses related to the inhibition of hypocotyl elongation, anthocyanin accumulation and cotyledon expansion, while CRY2 plays a role in the hypocotyl inhibition, circadian clock and photoperiod-dependent flowering [56]. However, CRY3 is a DASH protein located in chloroplasts and mitochondria [51], which works to repair UV-damaged DNA in a light-dependent manner [57]. Overall, the cryptochrome-mediated photoresponses remain unclear with the existing differences in plant species as well as their physiological responses [58,59].
It has been reported that plant cryptochromes were involved in the adversity stress response [22,24,60,61]. In our study, all the GhCRY genes identified were negatively related to the PEG, NaCl and high temperature treatments. The negative regulation pattern of CRY factors in response to the drought, salt and osmotic stresses is paobably common in plant species. For instance, in Arabidopsis CRYs play an important role in drought stress tolerance [21]. Overexpressing the CRY1 protein in Arabidopsis resulted in excessive water loss whereas cry-1cry2 double mutant plants were clearly more droughttolerant than the wild type. In addition, introducing Triticum aestivum CRYs (TaCRY1a and TaCRY2) into Arabidopsis plants reduced osmotic stress tolerance, including drought and salt stresses [27]. Meanwhile, overexpression of Sorghum bicolor SbCRY1a in Arabidopsis rendered the transgenic plants oversensitive to salt stress [28]. In tomato (Solanum lycopersicum L.), CRY1a modulated the water deficit response under osmotic stress conditions, further increasing tomato growth by reducing malondialdehyde (MDA) and proline accumulation [24]. In addition, the tomato cry1a mutant plants showed enhance drought tolerance wing to the increased leaf relative water content [25]. In Brassica napus overexpressing CRY1 resulted in plants that were very sensitive to osmotic stress, whereas the antisense silencing plants were more tolerant [60]. However, the relationship between CRYs and abscisic acid (ABA), and the role of other blue light photoreceptors in modulating water loss under drought or salt stresses is still unclear, thus further research is rewarding.
In the current study, nine GhCRY genes were found negatively related to the high temperature treatments, at least based on expression patterns. However, there were five GhCRY genes which were induced after a 12-h low temperature treatment, along with the repression of an equal number of the GhCRY s. These results demonstrated the functional divergence among the CRYs in G. hirsutum. Plant cryptochromes have been implicated in adaptations to the changing environmental factors, for instance, a report showed that low temperatures would increase the biological activity of CRY [62]. Consistently, as reported herein in cotton each CRYs may behave differently in response to low temperature stresses. Although there have been many studies on plant CRYs in regulating plant growth and development, as well as the response abiotic stresses, our understanding of cotton CRYs is still preliminary. This genome-wide survey paves the way for in-depth research of the function of cotton cryptochromes, which in turn will add valuable data to cotton breeding.

Conclusions
We systematically analyzed cotton CRY family genes and their expressions using bioinformatic approaches. We analyzed gene structures, chromosomal locations, intron-exon organizations, phylogenetic relationships and expression patterns in different cotton tissues and under different stress conditions to predict their possible biological functions. In particular, the GhCRY highly expressed in cotton fiber cells were identified. The decreased expressions of several GhCRY genes in response to multiple abiotic stress implies their involvement in the regulation of growth and development under the abiotic stress conditions. Together, our results provide candidate genes to facilitate the functional identification of the CRY genes in cotton that are important in modulating plant growth, development and stress tolerance.
The presence of conserved domains in each Arabidopsis and Gossypium gene was verified using the SMART conserved domain search tool [67] and Pfam databases [68].

Chromosomal location and gene structure analyses
Chromosomal locations for each of the above identified GhCRYs were extracted from the genome annotation gff3 file [41]. Chromosomal locations of the predicted GhCRYs were visualized using TBtools [69], and the exon-intron structure of each gene was displayed using the online tool GSDS 2.0 [70]. The number of amino acids, molecular weight (MW), and theoretical isoelectric point (pI) of putative GhCRYs proteins were determined using the ProtParam tool [71].

Sequence alignment, Ka, Ks and phylogenetic analyses
Complete protein-coding sequences for CRY genes from Gossypium and AtCRY were aligned using MAFFT with the G-INS-i algorithm [72]. The nonsynonymous substitutions rate (Ka) and synonymous substitution rate (Ks) were calculated using the DnaSP 6.0 [73]. The NJ phylogenetic tree was constructed using MEGA version 6.0 [74] by sampling 1000 bootstrap replicates.

Analysis of Cis-acting element in promoter regions of GhCRYs
The upstream sequences (1.5 kb) [75] of the GhCRYs genes were retrieved from G. hirsutum genome sequence based on the gene locations [41]. Then, the retrieved promoter sequences were submitted to PlantCARE [54] to identify the potential Cis-acting element.

Chromosomal mapping and synteny analysis of CRY genes in diploid and allotetraploid Gossypium species
CRY genes were mapped on chromosomes using TBtools [69] software. Blastn was used to determine CRY gene synteny. Then, TBtools [69] software was applied to express the syntenic relationship of the homologous gene pairs.

Expression patterns of GhCRYs in different tissues and stress conditions
Raw RNA-Seq data for G. hirsutum seed, root, stem, leaf, torus, petal, stamen, ovary, calyx, ovule (-3 dpa, -1 dpa, 0 dpa, 1 dpa, 3 dpa, 5 dpa, 10 dpa, 20 dpa, 25 dpa, 35 dpa) and fiber (5 dpa, 10 dpa, 20 dpa, 25 dpa) were downloaded from the NCBI Sequence Read Archive (PRJNA 248,163) [41], represented by one library each. Reads were mapped to the G. hirsutum genome [41] via HISAT2 software with default parameters, and read abundance calculated via StringTie [76,77]. Read counts were normalized in R3.2 using RUVSeq [78] and the internal control reference gene GhUBQ7, which is detected at relatively constant levels across different cotton samples [79]. Potential batch-effects were corrected by an improved version of ComBat, ComBat-seq [80]. Gene expression was estimated by Ballgown [81], using fragments per kilobase million (FPKM) values to calculate the gene expression levels across libraries. Expression levels of G. hirsutum leaf RNA-Seq data (in FPKM) for each GhCRY gene under drought, salt, heat and cold stress (time points: 0, 1, 3, 6, 12 h) were retrieved from the ccNET database [82]. Genes were considered differentially expressed if expression varied more than twofold change with a p-value of less than 0.05. TBtools [69] was used to display the gene expression patterns from the calculated FPKM values.

Plant cultivation and stresses treatment
G. hirsutum cv. R15 [83], were planted in a controlled environment at 28 °C day/20 °C night, with a 16-h light/8-h dark photoperiod [84]. For stress treatments, 28-day old plants were treated with 200 mM NaCl, PEG6000, 42 °C (hot) and 4 °C (cold), respectively. Leaves from stress-treated plants were collected at 12 h post-stresses treatments for further expression analyses. All plant tissues were frozen in liquid nitrogen immediately after collection and stored at -80 °C until RNA extraction. These treatments were sampled with three biological repeats.

RNA extraction, cDNA synthesis and qRT-PCR expression analyses
Total RNAs from cotton tissues were extracted using the RNAprep pure plant kit (TIANGEN, Shanghai, China) as described [85]. The resulting RNAs were treated with DNase I prior to synthesizing cDNA with oligo (dT) primers and M-MLV Reverse Transcriptase (Invitrogen); these products were diluted fivefold before use. For quantitative real-time PCR (qRT-PCR), Primer5 software was used to design gene-specific forward and reverse primers (Additional file 4: Table S2). Analyses were performed with SYBR-Green PCR Mastermix (TaKaRa) on a cycler (Mastercycler RealPlex; Eppendorf Ltd, Shanghai, China). The G. hirsutum histone-3 (GhHIS3) genes were used as internal references [81,86], and the relative amount of amplified product was calculated following the 2 −∆∆Ct method [87]. For the G. hirsutum samples, relative expression levels among different stresses were normalized by calibrating with the leaves sample from the wild plants. The leaves sample was washed with DEPC sterile water three times before extracting the RNA.