Global gene expression in cotton (Gossypium hirsutum L.) leaves to waterlogging stress

Cotton is sensitive to waterlogging stress, which usually results in stunted growth and yield loss. To date, the molecular mechanisms underlying the responses to waterlogging in cotton remain elusive. Cotton was grown in a rain-shelter and subjected to 0 (control)-, 10-, 15- and 20-d waterlogging at flowering stage. The fourth-leaves on the main-stem from the top were sampled and immediately frozen in liquid nitrogen for physiological measurement. Global gene transcription in the leaves of 15-d waterlogged plants was analyzed by RNA-Seq. Seven hundred and ninety four genes were up-regulated and 1018 genes were down-regulated in waterlogged cotton leaves compared with non-waterlogged control. The differentially expressed genes were mainly related to photosynthesis, nitrogen metabolism, starch and sucrose metabolism, glycolysis and plant hormone signal transduction. KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis indicated that most genes related to flavonoid biosynthesis, oxidative phosphorylation, amino acid metabolism and biosynthesis as well as circadian rhythm pathways were differently expressed. Waterlogging increased the expression of anaerobic fermentation related genes, such as alcohol dehydrogenase (ADH), but decreased the leaf chlorophyll concentration and photosynthesis by down-regulating the expression of photosynthesis related genes. Many genes related to plant hormones and transcription factors were differently expressed under waterlogging stress. Most of the ethylene related genes and ethylene-responsive factor-type transcription factors were up-regulated under water-logging stress, suggesting that ethylene may play key roles in the survival of cotton under waterlogging stress.


Introduction
Waterlogging is a well-known abiotic stress that can severely damage crop production worldwide. It usually occurs in areas with high rainfall and/or poor drainage [1]. Low oxygen is known to be a serious factor in the adverse effects of waterlogging [2]. The low oxygen condition leads to less adenosine triphosphate (ATP) production and has negative impacts on PLOS ONE | https://doi.org/10.1371/journal.pone.0185075 September 27, 2017 1 / 24 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 cellular energy status. Therefore, a key feature for the acclimation to low oxygen environment is to activate genes encoding proteins and enzymes for anaerobic fermentation, glycolysis, transcription factors, and signaling pathways in order to allow biological and physiological adjustments to the low oxygen conditions [3].
Waterlogging usually results in a rapid decrease in photosynthetic rate, ranging from 10 to 90% in different plant species [19]. Thus, waterlogging response of plants has been strongly associated with regulation of genes related to photosynthesis and photosystem II (PSII). The expression of Chlorophyll a-b binding protein 4 (LHCB4), a gene involved in the light-harvesting complex of PSII, was specifically reduced in waterlogged sesame [20]. Our previous study also showed that the expression of GhLHCB was significantly down-regulated by waterlogging stress [21].
Plant hormones play important roles in regulating developmental processes and signaling networks involved in plant responses to abiotic stresses. The induction of hormones under waterlogging stress is involved in signaling cascades, including increases in ethylene [22,23], abscisic acid (ABA) [24] and gibberellic acid (GA) [25], and a reduction in cytokinin (CK) [26] and auxin (IAA) [10]. ABA is an important signal which can be induced by some abiotic stress in regulating stomatal conductance and photosynthesis and transpiration [27]. An increase in ABA concentration was noted in leaves of flooded alfalfa [28]. Bai et al. [29] also found an increased ABA content in Malus leaves under hypoxia stress, indicating that ABA is a key signal in mediating responses to waterlogging.
Transcriptional factors (TFs) also play an important role in waterlogging response. Two ethylene response factors (ERFs) of Jatropha, JcERFVII-2 and JcERFVII-3, were noticeably induced in response to waterlogging [30]. The Snorkel (SK) and Submergence-1A (Sub-1A) are two pivotal genes for rice to adapt to different types of waterlogging; although both genes encode ethylene-responsive factor-type transcription factors, they function in opposite ways [31,32]. Studies of Arabidopsis have shown that group VII ERF participated in oxygen sensing through the N-end rule pathway [33]. HRE1, HRE2, RAP2.2, and RAP2.12 are four members of group-VII ERFs; overexpression of these four genes significantly improved low oxygen survival by promoting the expression of genes involved in low oxygen adaptation [34][35][36]. Like ERF, the expression of MYB and WRKY has also been regulated in plant under waterlogging stress, which activates many vital processes to cope with stress [37].
Cotton is a crop sensitive to waterlogging. In recent years, the growth response and yield loss of cotton to waterlogging stress have been extensively investigated [38][39][40][41][42]. Waterlogging stress decreased photosynthetic rate, stomatal conductance, and plant height, as well as biomass and lint yield [40][41][42][43]. Bange et al. [42] reported that intermittent waterlogging caused an approximately 38% reduction in cotton yield. Hodgson [44] showed that waterlogging stress of 4, 8, 16 and 32 h significantly decreased cotton yield. When subjected to waterlogging for 3, 6, 9 and 12 d, cotton suffered a lint yield reduction of 16, 24.1, 39.5 and 50.2%, respectively [45]. Our previous work also showed considerable yield loss due to waterlogging for 10 d, 15 d and 20 d at flowering stage [21]. Although the morphological and physiological changes including plant growth reduction and yield loss under anaerobic conditions have been well documented, limited work has been done to elucidate changes and adaptations of waterlogged cotton at the molecular level. In this study, we performed a transcriptome analysis using RNAseq to investigate gene expression patterns in the waterlogged cotton leaves to provide new insights into the molecular responses of cotton to waterlogging under field conditions.

Plant material and waterlogging treatment
Experiments were conducted in an electrical rain-shelter at the Experimental Station of Shandong Cotton Research Center, Linqing (115˚42'E, 36˚61'N), Shandong, China, in 2015. Twelve bottomless pools (3 m × 4 m) separated by concrete walls (13 cm thick and 1.5 m deep) to prevent soil water movement in different pools were established under the rain shelter. A commercial cotton (Gossypium hirsutum L.) cultivar, K638 developed by the Cotton Research Center, Shandong Academy of Agricultural Sciences, Jinan, was used in the experiment. Aciddelinted seeds of K638 were sown and allowed to grow in bottomless pools.
Waterlogging stress was established by over irrigation at peak flowering. The over irrigation of waterlogged plots was sustained till the water layer reach 20 cm above the soil surface for 10, 15 or 20 days. Plots with soil moisture of 60-70% under normal irrigation acted as non-waterlogging controls. The experiment was arranged into a randomized complete block design with 3 replications. The 10-, 15-and 20-d waterlogging treatment started on 15, 10 and 5 July, 2015. At the end of each waterlogging treatment, the surface water in each plot was manually removed and the soil moisture recovered to normal in a few days through evaporation.
After 10-, 15-and 20-d treatments, the fourth leaves on the main-stem from the top were harvested, frozen in liquid nitrogen and stored at -80˚C for physiological measurement and RT-PCR analysis. The samples of each biological replicate were pooled from 6 plants to avoid any potential effects of position within the field. Leaves of 15-d waterlogged cotton plants were analyzed by RNA-Seq for global gene expression study.
The procedures of field management such as sowing, seedling, fertilizer application, plant pruning, pest control and chemical control were implemented according to local agronomic practices as described in Zhang et al. [21].

RNA extraction and Solexa sequencing
RNA isolation Kit (Hua yue yang Biotechnology, China) was used for total RNA extraction. For Solexa sequencing, total RNA from 6 representative individual plants of each treatment was mixed into one biological replicate. Approximately 20 mg of total RNA was used. Tag libraries were prepared using the Illumina Gene Expression Sample Prep Kit, following the manufacturer's protocol, as described in Luan et al. [46]. The libraries were then sequenced using an Illumina HiSeq 2500 with 50-bp single-end (SE) reads each.

Gene annotation
The database of G. hirsutum L. genome (http://mascotton.njau.edu.cn/html/Data/Genome fhsequence/2015/04/20/16ab0945-19e9-49f7-a09e-8e956ec866bf.html) was used as a reference sequence to align and identify the sequencing reads. To map the reads to the reference, the alignments and the candidate gene identification procedure were conducted using the mapping and assembly with qualities software package [47]. The reference sequences were converted to binary FASTA format, and each Solexa read data subset (corresponding to one lane on the instrument) was transformed from Solexa FASTQ to Sanger FASTQ format. As recommended, each subset was separately mapped to the reference. These maps were then merged to form general maps for assembling the consensus sequences [48].

Identification of differentially expressed genes (DEG) and functional analysis
The RPKM (Reads Per kb per Million reads) method was used to normalize the influence of different gene length and sequencing level on the calculation. Thus the result can be used directly for comparing the gene expression of different samples. The gene expression level was calculated by the formula [49]: Where C represents the number of reads that uniquely aligned to the gene, N is the total number of reads that uniquely aligned to all genes in the specific sample, and L is number of bases of the gene. The p-value corresponding to differential transcript expression in two samples was determined from Audic's algorithm [50], and False Discovery Rate (FDR) method was applied to determine the threshold of P-values in multiple tests. The DEGs were obtained after filtering using a FDR of 0.001 and an absolute value of log 2 Ratio!1.
GO enrichment analysis was performed for functional categorization of differentially expressed transcripts using agriGO software [51] and the P-values corrected by applying the FDR correction to control falsely rejected hypotheses during GO analysis. The pathway analysis was conducted using KEGG (www.genome.jp/kegg/).

Real-time PCR (RT-PCR) analysis
To validate the results of the DGE-based analyses, the expression of 30 genes was determined using quantitative RT-PCR. The leaves from 9 representative individual plants of waterlogged and control plant were harvested and every 3 leaves were combined into one biological replicate and then extracted for total RNA. The samples used for qRT-PCR experiments and RNA-Seq analysis were not the same. Quantitative RT-PCR was performed according to the instructions provided for the Bio-Red iCycler iQ system. Each sample was run in triplicate on Bio-red IQ2 Sequence Detection System and Applied Biosystems software using first-strand cDNAs and SYBR Green PCR Master Mix. The first-strand cDNA was synthesized by Superscript II reverse transcriptase (Invitrogen) following procedures described in the manufacturer's guidelines. Gene-specific primers were designed according to the gene sequences using the Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA) and then synthesized commercially (Shanghai Sangon Biological Engineering Technology & Services Co., Ltd., Shanghai, China). The specific primers for the selected genes and internal control gene (actin) used for qRT-PCR are listed in Table 1. The amplification of β-actin was used as an internal control to normalize all data. Thermal cycling was performed at an initial denaturation step at 95˚C for 3 min followed by 40 cycles at 95˚C for 10 s, at annealing temperatures of 60˚C for 10 s, and at 72˚C for 10 s.

Physiological measurements
Net photosynthetic rates of the fully expanded young leaf (4th leaf) on the main-stem from terminal were taken between 09:00 and 11:00 h on cloudless days, using a LI-6400 portable photosynthesis system (Li-Cor, Lincoln, NE, USA) with an integrated fluorescence chamber head (LI-COR 6400-40 Leaf Chamber Fluorometer). The photosynthesis-related parameters were as follows: PAR 1500 μmol m −2 s −1 at the leaf level; the leaf temperature of the instrument was set to 25˚C; the CO 2 concentration was 400 μmol mol −1 ; and the relative humidity inside the LI-COR leaf chamber varied between 45 and 60%. The 4th main-stem leaves of 9 randomly selected plants were sampled and washed with distilled water, and immediately frozen in liquid nitrogen and stored at −80˚C for physiological measurement.
Nitric oxide was detected with the nitric oxide (NO) assay kit (Nan-jing Jiancheng Bioengineering Institute, China) according to the manufacturer's guidelines. The content of soluble protein and soluble sugar was measured following procedures described in the manufacturer's guidelines of their assay kits (Nan-jing Jiancheng Bioengineering Institute, China).
The 4th main-stem leaves from other 9 randomly selected plants per waterlogged plot were oven-dried for 72 h at 80˚C and milled to analyze for N by Kjeldahl method.

Statistics
Analysis of variance was performed with the function of completely randomized design using Data Processing System (DPS) [52]. Means were separated using Duncan's multiple range tests at p = 5%.

RNA-Seq analysis and identification of differentially expressed genes
RNA-Seq analysis was performed to determine the transcriptome response to waterlogging stress in cotton. We generated 11.7 and 11.4 million reads from libraries of non-waterlogged and waterlogged cotton plants, and obtained 10.9 and 10.8 million clean tags from non-waterlogged and waterlogged cotton libraries (SRA submission number: SRP095435) after the low quality tags were filtered out. The gene sequences of G. hirsutum genome were used as reference to align and identify the sequencing reads. This allowed for the mapping of approximately 90% of the distinct clean tags that passed our filters, representing more than 10.8 million reads per library ( Table 2). Putative differentially expressed genes were finally selected depending on the expression profiles and whether: a) the average fold change between two treatment genes was more than or equal to two folds; and b) the FDR was less than 0.001. Transcriptome analysis identified 1812 DEGs with 794 up-regulated genes and 1018 down-regulated genes.

Functional classification of differentially expressed genes
Gene ontology (GO) analysis was performed by mapping each differentially expressed gene into the records of the GO database (http://www.geneontology.org/). The GO annotation of these genes is presented in Fig 1. Waterlogging stress affected a wide spectrum of physiological processes as described in the relevant GO analysis. Most of the genes related to hormone response, auxin-activated signaling pathway, regulation of hormone levels, defense response, gibberellin biosynthetic process, programmed cell death, abscisic acid glucosyltransferase activity and xyloglucan endotransglucosylase activity were up-regulated under waterlogging stress (Fig 1). On the contrary, most of the genes related to photosynthesis, light harvesting, light reaction, electron transport chain, cysteine biosynthetic process, carbohydrate metabolic process, xyloglucan metabolic process and hydrogen peroxide-mediated programmed cell death were down-regulated under waterlogging stress (Fig 1).
The GO terms were in accord with KEGG pathway analysis; e.g., enrichment of nitrogen metabolism, starch and sucrose metabolism, glycolysis and plant hormone signal transduction. In addition, significant changes in photosynthesis, flavone and flavonol biosynthesis, cysteine and methionine metabolism, steroid biosynthesis, glutathione metabolism, phenylpropanoid biosynthesis, spliceosome, as well as circadian rhythm pathways were found ( Table 3). The 44 differently expressed photosynthesis related genes were all down-regulated, and 19 of 20 differently expressed antenna proteins related genes were down-regulated under waterlogging treatment. Similarly, most of the nitrogen metabolism and amino acid metabolism related DEGs were down-regulated under waterlogging treatment. However, 24 of 28 differently expressed flavonoid biosynthesis related genes were up-regulated under waterlogging (Table 3).
After 10-, 15-and 20-d waterlogging, the expression levels of the 4 LHCB genes were all down-regulated (Fig 2C-2F). Consistent with the decreased expression of LHCB; the leaf Pn rate in waterlogged cotton was reduced by 22.3, 27.2 and 38.4% after 10-, 15-and 20-d waterlogging (Fig 3). After 10-, 15-and 20-d waterlogging, the concentrations of Chl a were reduced Waterlogging decreased soluble sugar contents by down-regulating the expressions of carbon metabolism related genes Waterlogging inhibited the expression of starch and sucrose metabolism related genes. However, the expressions of glycolysis related genes were up-regulated ( Table 4). The expression patterns of 2 ADH genes were analyzed by real-time PCR at 10-, 15-and 20-d waterlogging stress and their expression was increased by waterlogging stress (Fig 2L and 2M). The accumulation of soluble sugar in the main stem leaf was decreased by waterlogging stress, being consistent with the decreased expression of sucrose metabolism related genes. The soluble sugar content in the main stem leaves of waterlogged cotton was reduced by 8.7, 16.0 and 14.6% at 10-, 15-and 20-d waterlogging treatment, respectively (Fig 4). Similarly, the soluble protein content was decreased by 15.8, 17.5 and 26.6%, respectively (Fig 4).
Waterlogging decreased the NO contents by down-regulating the expression of nitric oxide biosynthetic related genes Many genes involved in nitric oxide biosynthetic process were down-regulated under waterlogging (Table 4). Nitrate reductase is a key enzyme in catalyzing the conversion of nitrate (NO −3 ) to nitrite (NO -2 ) and nitric oxide (NO). Nitrite reductase is responsible for the reduction of NO −2 to ammonium and conversion of NO −2 to NO. These two genes were down-regulated under waterlogging. The NO concentration in the leaf was reduced by 76.2, 77.3 and 74.4% at 10-, 15-and 20-d waterlogging treatment, respectively (Fig 5). Hormone related genes and transcription factor genes There were 27 differently expressed hormone-related genes, of which 21 were up-regulated in waterlogged cotton compared with non-waterlogged control (Table 5). Compared with nonwaterlogged control, the 3 JA and 2 SA related DEGs were all up-regulated, but the 3 IAA related DEGs were down-regulated in waterlogged cotton ( Table 5). The 9 differently expressed GA-related genes, which include 4 GA biosynthesis genes, were all up-regulated in leaf of waterlogged cotton compared with non-waterlogged control ( Table 5). Five of the 8 differently expressed ABA-related genes were up-regulated in waterlogged cotton compared with nonwaterlogged control (Table 5). Interestingly, the 2 ABA biosynthesis genes, 9-cis-epoxycarotenoid dioxygenase (NCED), were all down-regulated, but the 2 ABA degradation genes, abscisic acid 8'-hydroxylase (CYP707A), were all up-regulated under waterlogging stress (Table 5).
To determine if GA-related genes in leaves of waterlogged cotton were differentially regulated as described in RNA-Seq data and check their expression patterns under different duration of waterlogging stress, the expression patterns of 2 GA biosynthesis genes (GA3ox1, GA3ox2) and GA receptor genes (GID1B, GID1-3) were analyzed by real-time PCR after 10-, 15-and 20-d waterlogging stress. The expression level of GA3ox1, GA3ox2, GID1B and GID1-3 was up-regulated by waterlogging stress. The expression of these genes in cotton leaves increased greatly after waterlogging stress and most of them peaked after 20 d-waterlogging (Fig 2G-2K). Thirty-four waterlogging-regulated TFs were identified, with 26 up-regulated and 8 down-regulated (Table 6). These TFs including 13 ERFs, 15 MYBs and 6 WRKYs with 9 ERF, 13 MYB and 4 differentially expressed WRKY TFs were up-regulated by waterlogging  stress. The expression level of three AP2/ERFs was up-regulated after 10, 15 and 20 days stress and two of them peaked after 15-d stress, the other one peaked after 10-d waterlogging stress (Fig 2N-2Q).

Confirmation of Solexa expression patterns by RT-PCR analysis
To test the reliability of Solexa sequencing, RT-PCR analysis was performed with specific primers for 30 genes, which have been identified by Solexa sequencing in which 13 genes were up-regulated and 17 genes were down-regulated. The results showed that 28 of the 30 genes had the same expression profiles as the original Solexa sequencing, and the original Solexa pattern was validated in 93.3% of the cases. This was not the case for other genes presumably because of the mutations within the primer sites or because the RNA used for Solexa sequencing and qRT-PCR were extracted from different plants. The expression patterns of the 30 genes were highly consistent with the Solexa sequencing ratios, with a relative R 2 of 0.8973 (Fig 6).

Discussion
Waterlogging has become one of the serious problems limiting cotton yield in recent years [43,44]. Understanding its effects and the mechanisms would help minimize yield loss. Our previous study showed that waterlogging at flowering stage significantly inhibited plant growth and reduced biological and economical yield of cotton due to reductions in net photosynthetic rate, leaf area, boll density and boll weight [53]. In this study, sequencing method was employed to compare differential gene expression profiles of cotton leaves subjected to 15-d waterlogging stress.

Global gene transcription changes in waterlogged cotton leaf
Waterlogging stress increased the expression of genes in glycolysis and some catabolism pathways, but reduced the expression of synthesis pathways and oxidative phosphorylation genes in poplar and maize [5,8]. The down-regulation of these energy-related and O 2 -consuming metabolic pathways suggested that plants initiate several responses to alleviate the impact of low O 2 during waterlogging periods. Our data showed that many genes with potential roles in carbohydrate metabolism, nitrogen metabolism and photosynthesis were down-regulated; in ethylene synthesis and perception, regulation of transcription and glycolysis were significantly up-regulated in cotton under waterlogging stress. Waterlogging stress decreases leaf photosynthesis, and many genes involved in the photosynthesis pathway were down-regulated [54,55]. Thus plant growth decreased with waterlogging stress due to decreasing photosynthetic rate following induced damage to cellular and photosynthetic machinery [56]. Our data showed that many genes involved in photosynthesis were down-regulated in leaves under waterlogging. Most of the genes related to chlorophyll and light-harvesting complex were down-regulated in leaves under waterlogging. Remarkably, 4 chlorophyll a/b-binding (LHCBs) genes which involved in the light-harvesting complex of photosystem II (PSII) were significantly down-regulated by waterlogging stress (Fig 2C-2F). The decreased plant growth under waterlogging stress may be due to the suppressed expression of photosynthesis and metabolism related genes.   Carbon and energy metabolism as affected by waterlogging Plants subjected to waterlogging stress shift their metabolism from oxidative phosphorylation to anaerobic fermentation to maintain ATP production by down-regulating storage metabolism, changing Suc synthase to Suc hydrolysis, and inhibiting mitochondrial respiration [3,[57][58]. Under waterlogging conditions, O 2 limits oxidative phosphorylation and plant cells must alternate metabolic pathways to produce ATP. Plant ethanolic fermentation is activated under low-oxygen stress: PDC firstly converts pyruvate to acetaldehyde, and then ADH converts acetaldehyde to ethanol. The expression of PDC1, PDC2 and PDC4 is strongly up-regulated under waterlogging in rice, which improves the tolerance under long-term anoxia [59]. Hypoxia and anoxia induced significantly the expression of Arabidopsis PDC1and PDC2, and  these two genes play key roles in the tolerance of submergence by mutant and transgenic experiments [60]. ADH is an important enzyme involved in ethanolic fermentation and induction of ADH during anaerobiosis has been observed in many plant species such as mungbean [61], pigeonpea [62] and rice [63]. ADH activity has been reported to increase under anoxia in all parts of plants such as roots [64], shoots [65,66], seedlings [67] and coleoptile [68]. The activity of ADH in flood-tolerant sorghum (Sorghum Bicolor L.) cultivar SSG-59-3 was significantly higher than the sensitive variety S-308 [69]. The maize mutant which is deficient in one of its ADH genes, was more sensitive to flooding injury and died earlier than wild type plants [70]. In addition, transgenic lines over-expressing ADH produced significantly more ethanol than wild type plants, indicating an increase in ethanol fermentation [71]. In the present study, the expression of ADH was specifically enhanced in waterlogged cotton, which was parallel to our previous observation that greater activity of ADH occurred in waterlogged cotton, suggesting that ADH may have important roles in maintaining cotton growth under waterlogging stress. In addition, many genes involved in the synthesis of sucrose and starch were down-regulated (Table 4) in waterlogged cotton, which may be the reason for the reduction in total soluble sugar content (Fig 4).

Nitrogen metabolism as affected by waterlogging
Waterlogging results in root anoxia and increased denitrification, leading to significant N loss from soil and potential nitrous oxide (N( 2 )O) emissions [72]. Plant repressed nutrient uptake and inhibited biosynthetic processes in order to avoid the occurrence of complete anoxia under waterlogging stress [5]. In this study, the soluble protein and N contents (Fig 7) in the leaves of waterlogged cotton were reduced, which might be due to the increased expression of nitrogen metabolism related genes and decreased expression of amino acid biosynthesis related genes.
Nitric oxide (NO) is an essential endogenous signal molecule involved in multiple physiological processes in plants [73]. It can also act as a secondary messenger in environmental stress signal transduction [74,75]. It was reported that NO could increase the photosynthetic pigments and net photosynthesis rate of leaves [76]. It could facilitate anaerobic survival of plants [77][78][79]. Chen [80] and Song [81] showed that providing NO by spraying SNP (sodium nitroprusside), a NO donor, alleviated the damage to plants caused by waterlogging stress. It has been observed that long periods of low oxygen survival can be achieved when NO −3 is provided because NO can be synthesized from NO −3 via NR and NiR [77,[82][83]. In this study, the NO content in the waterlogged cotton leaves decreased (Fig 5), possibly due to the decreased expression of NR and NiR under waterlogging stress.

Hormone responses to waterlogging
Ethylene is considered to be the first warning signal in plants under hypoxia stress. Ethylene concentration increased under waterlogging stress to stimulate shoot elongation in wetland plants and it is likely to be a key player in the interactions among hormones like ABA and GA [84][85][86][87][88]. In this study, the ACC oxidase 1 (ACO1), a gene involved in ethylene synthesis, was up-regulated in cotton leaves (Table 5). Furthermore, 9 of the 13 differently expressed ERF TFs were up-regulated under waterlogging stress, suggesting that ethylene may play an essential role in waterlogging response of cotton. Gibberellin has been reported to enhance waterlogging tolerance in rice and Rumex palustris under submerged condition [89,90]. GA could induce the elongation of internode to bring out rice leaves from water surface for aerobic respiration [3,32,86,91]. GA is perceived by its nuclear receptors GA INSENSITIVE DWARF1s (GID1s), which then trigger degradation of downstream repressors DELLAs [92]. In the present study, the expression of GA biosynthesis genes and 3 GID1 genes were up-regulated by waterlogging stress. It is suggested that GA may have positive role in maintaining cotton growth under waterlogging stress.
IAA plays important roles in plant growth and development [93]. In the present study, 3 auxin-related genes were down-regulated in the leaves of waterlogged cotton (Table 4). Interestingly, changes in the transcript levels of these genes were frequently associated with changes in the IAA content in leaves of 15-d waterlogged cotton plant in our previous report [21]. Besides, plant hormones JA, SA and BR are well regulators of plant growth. Nguyen et al. [94] reported that they are possibly involved in a network of signaling cascades to help plants adapt to abiotic stresses. In this study, many differently expressed JA, SA and BR related genes were up-regulated in waterlogged cotton, indicating that JA, SA and BR are possibly involved in waterlogging response (Table 5). Many TFs are differentially regulated under waterlogging stress and have important roles in plant response to waterlogging stress. Thirty-four differentially regulated TFs were found in our experiment, suggesting that transcriptional regulation plays a key role in the waterlogging response in cotton. HRE and RAP type ERF genes were shown to be the important regulators in plant responses to hypoxia and anoxia [33- 36,95]. A more recent study showed that the Nend rule pathway of protein degradation acts as a homeostatic sensor of hypoxia in Arabidopsis through the regulation of key hypoxia-response TFs [33][34]36]. In addition, Licausi et al. [96] identified two Arabidopsis hypoxia-inducible ERF genes, HRE1 and HRE2, which could induce anaerobic gene expression and ethanol fermentation to enhancing plant tolerance of anaerobic stress. The Arabidopsis RAP2.2 and rice SUB1A, which play important roles in survival under hypoxia can be induced by ERF [32,34]. In this study, 13 differently expressed ERFs were identified and 9 of them were up-regulated under waterlogging stress ( Table 6). The increased expression of ERF may have important functions in enhancing plant tolerance of anaerobic stress through inducing some waterlogging tolerance related genes. A further analysis of waterlogging-induced TFs seems required to find more functional candidate genes involved in the adaption to waterlogging stress.
As mentioned above, several published studies of other plant species have provided some fundamental clues to plant response to waterlogging, which identified a common response including the up-regulation of ethanolic fermentative genes and the down-regulation of genes involved in photosynthesis, electron transport chain, and light harvesting. In the present study, cotton showed the up-regulation of genes involved in ethylene production and signaling which was commonly found in response to low oxygen in other plant species. And the alteration of carbohydrate metabolism and the up-regulation of AP2/ERF genes appeared to be a common response to low oxygen stress in plants. Although, plants response to waterlogging had so much in common, various processes were specifically regulated in each species. Besides, cotton is a species with indeterminate growth habit and large compensative ability after the removal of stress, which makes cotton different from other crops. For example, four AP2/ERF genes were up-regulated in waterlogged cotton, as shown in the present results (Fig 2). Although this result was consistent with the up-regulation of two Arabidopsis ERF genes, HRE1 and HRE2, in hypoxia conditions [96], the extent in cotton was far less than Arabidopsis. It indicates that ethylene could play a more important role in the hypoxia response of Arabidopsis than cotton. Additionally, contrary to what has been observed in submerged Jatropha [30], the expression of NR was down-regulated in waterlogged cotton in our experiment suggesting the different-regulation of the NR gene may be species-specific responses to waterlogging stress. Taken together, we propose that the difference between cotton and other plants may also be attributed to the duration or the depth of waterlogging stress in the present study.

Conclusions
Conclusively, waterlogging stress resulted in a number of biological change such as reduced photosynthesis and decreased contents of chlorophyll, NO, total soluble sugar and protein. These biological changes in waterlogged cotton were attributed to the decreased expression of photosynthesis related genes, and increased expression of glycolytic pathway and fermentation genes. Furthermore, the changes in the expression pattern of these genes might be regulated by synthesis and perception of plant hormones like ethylene and GA, and some transcription factors like ERF, MYB and WRKY (Fig 8).