Skip to main content
  • Research article
  • Open access
  • Published:

Global analysis of the Gossypium hirsutum L. Transcriptome during leaf senescence by RNA-Seq



Leaf senescence is an important developmental programmed degeneration process that dramatically affects crop quality and yield. The regulation of senescence is highly complex. Although senescence regulatory genes have been well characterized in model species such as Arabidopsis and rice, there is little information on the control of this process in cotton. Here, the senescence process in cotton (Gossypium hirsutum L.) leaves was investigated over a time course including young leaf, mature leaf and leaf samples from different senescence stages using RNA-Seq.


Of 24,846 genes detected by mapping the tags to Gossypium genomes, 3,624 genes were identified as differentially expressed during leaf senescence. There was some overlap between the genes identified here and senescence-associated genes previously identified in other species. Most of the genes related to photosynthesis, chlorophyll metabolism and carbon fixation were downregulated; whereas those for plant hormone signal transduction were upregulated. Quantitative real-time PCR was used to evaluate the results of RNA-Seq for gene expression profiles. Furthermore, 519 differentially expressed transcription factors were identified, notably WRKY, bHLH and C3H. In addition, 960 genes involved in the metabolism and regulation of eight hormones were identified, of which many genes involved in the abscisic acid, brassinosteroid, jasmonic acid, salicylic acid and ethylene pathways were upregulated, indicating that these hormone-related genes might play crucial roles in cotton leaf development and senescence. However, most auxin, cytokinin and gibberellin pathway-related genes were downregulated, suggesting that these three hormones may act as negative regulators of senescence.


This is the first high-resolution, multiple time-course, genome-wide comprehensive analysis of gene expression in cotton. These data are the most comprehensive dataset currently available for cotton leaf senescence, and will serve as a useful resource for unraveling the functions of many specific genes involved in cotton leaf development and senescence.


Leaf senescence is the ultimate phase of plant leaf development, and involves coordinated action at the cell, tissue, organ and organism levels. Senescence is controlled by a highly regulated genetic program and leads to death or the end of the life span [1]. Leaf senescence is not a passive and unregulated degeneration process; during senescence leaf cells undergo dramatically coordinated changes in cell structure, metabolism and gene expression [2]. The earliest and most significant cellular degeneration process begins with the chloroplast, which contains most of the proteins in a leaf cell. The mitochondria and nucleus remain intact until the final stages of leaf senescence. Metabolically, carbon assimilation, such as photosynthesis, is replaced by catabolism of chlorophyll and macromolecules, such as proteins, membrane lipids and nucleic acids [3]. Although leaf senescence is a deleterious process, it is crucial for the fitness of plants; efficient senescence is essential to maximize viability during the plant’s life cycle. Senescence can also occur prematurely when plants are stressed; it is a protective mechanism, leading to decreased yield and quality in crop plants by limiting the growth phase, which has become an increasing concern because of global climate change in recent years [4].

Recently, many advances have been made in understanding leaf senescence at the molecular level by the identification and characterization of a large number of senescence-associated genes (SAGs) and senescence-related mutants in various plant species, including plants such as Arabidopsis thaliana, Oryza sativa, and Medicago truncatula [5-10]. Among these SAGs, numerous transcription factors (TFs) such as NAC, WRKY, MYB [11], kinases and receptor-like kinases [12], signal transduction-related proteins, and regulators of metabolism are involved in regulating leaf senescence, indicating that senescence is an integrated response to many signals that are governed by highly complex transcriptional regulatory networks.

Cotton (Gossypium hirsutum L.) is one of most important economic crops, and is widely cultivated for the value of its fiber. In cotton, leaf senescence may occur too early or too late under the influence of certain internal factors or uneven environmental stresses [13]. Late senescence would affect remobilization of nutrients from the leaf to newly formed sink organs, such as developing boll. Whereas, premature senescence has occurred with increasing frequency since the introduction of modern, high-yielding cotton cultivars, such transgenic Bacillus thuringiensis (Bt) cotton, where premature senescence frequently developed during the period of rapid boll filling. Premature senescence results in reduced lint yield and poor fiber properties [14,15]. Previous studies on leaf senescence in cotton were limited mainly to physiological and cytological studies. The initiation and progression of leaf senescence in cotton can be affected by internal factors, such as phytohormones, mainly cytokinin (CK) and abscisic acid (ABA) [16], or various environmental factors such as drought [17], nutrient deficiency [15,18], salinity [19], UV-B radiation and elevated CO2 [20,21]. However, the mechanisms and the global transcriptional regulation of leaf senescence in cotton remain poorly understood.

To date, many molecular approaches have been used to understand the biochemical pathways and regulatory mechanisms associated with leaf senescence. Many novel SAGs were identified with distinct expression profiles during Arabidopsis developmental leaf senescence and induced senescence; comparison of changes in gene expression patterns in different senescent systems indicated differences in the genetic programs, such as nitrogen mobilization [12,22]. A DNA microarray comprising 13,490 expressed sequence tags was used to investigate gene expression during leaf senescence in Populus; a major shift in gene expression from photosynthetic competence to energy generation was observed [23]. Recently, analysis of DNA microarrays analysis showed a high-resolution time-course profile of gene expression during Arabidopsis leaf senescence [24]. The next-generation high-throughput sequencing technology, referred to as RNA-Seq, has emerged as a revolutionary tool to better understand complex eukaryote transcriptome [25]. In this approach, there is no strict requirement for a reference genome sequence [26]; therefore, it is suitable for the study of non-model species whose entire genomic sequences are not available. RNA-Seq has been used to investigate gene expression in many plants [27-29]. In the present study, RNA-Seq was performed on samples from different stages of leaf development, and the global changes in the developing leaf transcriptome were analyzed. Based on extensive data analyses, many genes and pathways were identified and characterized, which involved in leaf development and senescence in cotton. Analysis of TFs and hormone-related genes revealed an integrative image of gene regulatory networks. This work will serve as the foundation for further study of the molecular regulation mechanisms of leaf senescence in G. hirsutum and represents an invaluable resource for further identification of genes involved in leaf senescence.


Biochemical changes during senescence and selection of the analysis period

In this study, newly emerging cotton fourth leaves of the same size were tagged; the day when the newly emerging leaf expanded was considered as the first day. As shown in Figure 1A, leaves started to show yellowing at the tip at around 35 d and leaves were harvested at defined time points until 85 d, when they were visibly senescent (about 85% of the leaf area was yellow). Malondialdehyde (MDA) and chlorophyll levels are often used as markers for the progression of senescence because the levels of both change significantly during senescence. Levels of total chlorophyll and MDA were measured (Figure 1B,C). Chlorophyll levels increased from 5 d until 25 d, reaching a maximum at 25 d, and then declined. The MDA level did not change significantly until 55 d and then continued to increase. Based on the phenotype and physiological tests (Figure 1), leaf samples at 15, 25, 35, 45, 55 and 65 d were selected for digital gene expression (DGE) sequencing analysis.

Figure 1
figure 1

Plant phenotypes, and chlorophyll and MDA measurements. A. Appearance of G. hirsutum at eight different stages of development. B. Chlorophyll levels were measured in leaf samples at each stage of development. FW, fresh weight. Error bars represent the SD from three independent experiments. C. MDA contents were measured in leaf samples at each stage of development. FW, fresh weight. Error bars represent the SD from three independent experiments.

RNA-Seq analysis over multiple time points and the identification of differentially expressed genes (DEGs)

Six DGE libraries were sequenced, representing time points at 15, 25, 35, 45, 55 and 65 d, which generated approximately 4.8 million raw tags in each library (SRA submission number: SRR656612). After filtering out the low quality tags, the total number of clean tags per library ranged from 4.3 to 4.7 million. To reveal the molecular events underlying the DGE profiles, the tag sequences of these six DGE libraries were mapped to non-redundant reference sequence (the reference sequence is composed of 47% of G. raimondii genome and 53% of G. arboretum genome). Approximately 30% of all distinct tags could be unambiguously mapped to a gene in the reference database and 60% of the distinct clean tags were mapped to the Gossypium genome database. 24,846 genes were detected as expressed during the cotton leaf development; among these, 3,624 DEGs were obtained (Additional file 1).

Confirmation of tag-mapped genes by qRT-PCR

To validate the results of the gene expression analysis obtained using DGE data, qRT-PCR was performed on selected genes with different expression levels and functional assignments. Transcriptional regulation revealed by RNA-Seq was confirmed in a biologically independent experiment using qRT-PCR. Of the 25 candidates, four genes showed either no specific amplification or unexpected size amplification. The Pearson correlation coefficient (calculated by SAS software) was used to assess the correlation between different platforms. Overall, the outcome of qRT-PCR was in very good agreement (the mean correlation coefficient of 0.86) with the DGE results (Figure 2). Thus, our data demonstrated that the DGE technique for counting transcripts gives an accurate reflection of transcript abundance, and can be used for gene expression analysis in an organism lacking genome information.

Figure 2
figure 2

qRT-PCR validation of selected transcripts in cotton leaves. The columns represent the expression level of the transcript (left y-axis); the lines represent the relative intensity of real-time qRT-PCR from independent biological replicates (right y-axis).

Comparison of SAGs between cotton and other plants

To analyze SAG orthologs between cotton and other plants, identified DEGs were compared with those in other species in the leaf senescence database (LSD, June, 2013; using BLASTX searches (e-value ≤10−10). These genes were grouped into several categories according to the function and classification of the most closely related homolog. 2, 414 out of the 3,624 DEGs have homologous genes in other plants in the LSD. Most of the genes were homologous to genes in Arabidopsis, closely followed by Banana. Interestingly, these orthologs were mainly involved in transcriptional regulation, protein degradation/modification, signal transduction, redox regulation, lipid/carbohydrate metabolism, nutrient recycling, hormone response pathway and transport (Additional file 2).

Clustering of DEGs by expression pattern illustrated metabolic changes occurring during leaf senescence

Genes with similar expression patterns are often functionally correlated. To screen novel candidate genes whose functions correlate with leaf senescence, the 3,624 DEGs were clustered by MultiExperiment Viewer (MeV, v4.7.4), based on the K-means method and hierarchical clustering, respectively. Forty-two clusters were obtained for the data for the six time points (Figure 3A), the cluster number and corresponding expression values for each DEG in the six time point clusters are shown in Additional file 1. A heat map (Figure 3B) was generated to illustrate the changes in gene expression at each time point. There are several time points at which an obvious step change in the transcriptome occurs. Overall, the major switch in gene expression, both in upregulated and downregulated genes, occurs between 25 d and 35 d. Similarly, the genes identified as differentially expressed can be divided into two major groups based on their expression patterns (Figure 3). Type I includes 1,847 genes from clusters 1 to 25 that were downregulated during this period, and Type II includes 1,434 genes from clusters 26 to 39 that were upregulated. Some of the clusters showed a more complex pattern; for example, genes in clusters 40, 41, and 42 showed an initial increase followed by a decrease in expression (Figure 3).

Figure 3
figure 3

MeV cluster analysis of differentially expressed genes. A. Cluster analysis by the K-means method from the gene expression profiles. Selected enriched gene ontology (GO) terms (data shown in Additional file 3) are indicated; green indicates downregulated and red indicates upregulated genes. B. Heat map illustrating the expression profiles of the genes in each cluster.

To permit the global analysis of gene expression profiles and reveal the major functional categories represented in the DEGs during cotton leaf development, the two groups of genes showing either decreasing (clusters 1–25) or increasing (clusters 26–39) expression during leaf development were analyzed by GO enrichment tool BiNGO [30]. The BiNGO-derived graph (Figure 4) indicated the most highly significant enrichment of specific functions. Downregulated genes were significantly enriched for genes linked to cellular components such as thylakoid, plastid, vacuole, membrane, cell wall and extracellular region, and there was significant enrichment for photosynthesis, biosynthetic process, generation of precursor metabolites and energy in the biological process term. Within the molecular function category, catalytic activity and structural molecule activity were overrepresented. All these functions are essential for leaf growth, but are downregulated during senescence. The loss of integrity of cellular structures, such the membrane, thylakoid, plastid and cell wall leads to disruption of cellular homeostasis, ending the life of a cell in senescing leaves. The upregulated genes showed a very different picture, the most significantly overrepresented terms focused on the defense and protective steps that the plant takes to respond to the stress generated by a series of degradative processes and nutrition recycling. Within the biological process, response to endogenous stimulus, response to biotic stimulus, cell death, response to stress, protein modification process and signal transduction were overrepresented. For the molecular function, kinase activity, transferase activity, transcription factor activity and protein binding are overrepresented, whereas few cellular component terms are overrepresented. Meanwhile, to investigate the dynamic changes in gene expression and associated metabolic process, enriched GO terms in individual clusters were also analyzed (Figure 3A and Additional file 3). Unsurprisingly, chlorophyll biosynthesis genes and photosynthesis-related genes were significantly enriched; particularly photosynthesis-related genes, which were downregulated in many clusters. Clusters 18 and 19 are enriched genes for ribosome and cytosolic ribosome, cluster 21 contained genes involved in chlorophyll binding. Meanwhile, cluster 24 and 25 also contained many genes associated with cell wall. For the upregulated clusters, enrichment was seen for many genes involved in response to carbohydrate stimulus, defense response and protein kinase activity.

Figure 4
figure 4

GO terms enriched for genes upregulated and downregulated during leaf Senescence. The network graphs show BiNGO visualization of the overrepresented GO terms for the combined clusters of genes either upregulated or downregulated during senescence. Colored nodes represent GO terms that are significantly overrepresented; the colors are shaded according to the significance level as shown in the color bar.

Several biochemical pathways are significantly enriched during leaf senescence

All DEGs were analyzed using the KEGG Orthology Based Annotation System (KOBAS) to identify the metabolic pathways in which they function. A total of 1, 769 downregulated genes mapped to 100 KEGG pathways and 1, 357 upregulated genes mapped to 88 KEGG pathways (Additional file 4). Four of these pathways were significantly downregulated, whereas five of these pathways were significantly upregulated (corrected p-value ≤ 0.05) during leaf senescence (Table 1). Many typical senescence symptoms such as photosynthesis-antenna proteins, photosynthesis, porphyrin and chlorophyll metabolism, and carbon fixation were significantly overrepresented in the downregulated pathways. Certain protective activities pathways were significantly enriched in the upregulated genes, such as plant-pathogen interaction, fatty acid elongation, glutathione metabolism, taurine and hypotaurine metabolism. Meanwhile, plant hormone signal transduction pathway was obvious activated in the upregulated genes. These results suggested that these pathways might play an important role in leaf senescence.

Table 1 Significantly enriched pathways in cotton identified by KOBAS during leaf senescence (Corrected P-value ≤ 0.05)

TF families are active at different times during senescence

TFs are key regulatory proteins that are essential for regulation of gene expression. TFs are known to play significant roles in leaf senescence. Among the 3,624 DEGs, 519 TF transcripts from 50 TF families were identified (Additional file 5). The WRKY (54), bHLH (44), C3H (43), NAC (42), AP2-EREBP (25), FAR1 (24), DBP (23), SET (22), MYB (21) and HB (17) families were the top ten largest families active during leaf development, some of which are critical components of plant adaptive response to biotic, abiotic stresses and senescence. A heat map depicting the overall trend of the differential expression profiles of the TF genes during leaf senescence was constructed using MeV. As shown in Figure 5, many TF families were significantly upregulated. For example, the large WRKY and C3H family were significantly upregulated, with many members of these families being upregulated early, from around 25 d. Another significantly overrepresented upregulated family was the large NAC family, with over 40 of the members of this family showing altered expression at various times during senescence. Many members of the NAC family have been implicated functionally in senescence and in a variety of stress-related programs [11,23,31-34]. In addition, bHLH, GRAS, MYB, DBP and AP2-EREBP TFs are differentially regulated during senescence and were also upregulated during the early stages in our data [11,22-24,32].

Figure 5
figure 5

Hierarchical cluster analyses for cotton leaf putative TFs differentially expressed during leaf senescence. Various transcription factor families showing differential expression are shown on the right side. The x-axis refers to each stage of leaf development. The y-axis represents the expression levels of the genes. Numbers in parentheses represent the number of transcription factor family members.

Differential expression of hormone regulation-related genes during leaf Senescence

To identify the contribution of hormone-mediated transcriptional regulation to changes in patterns of gene expression during leaf senescence, we mapped the DGE transcripts to eight hormone-related pathways in the Arabidopsis Hormone Database, including ABA, auxin, ethylene (ET), brassinosteroid (BR), jasmonic acid (JA), CK, GA and salicylic acid (SA). 960 genes associated with various aspects of hormone homeostasis, such as biosynthesis, response, signaling, receptors and metabolism exhibited significant differential expression during senescence. (Table 2 and Additional file 6).

Table 2 Hormone-related genes differentially expressed in leaves during leaf senescence

The largest group that showed significant differential expression was related to ABA (Figure 6A), which is involved extensively in response to environmental stresses and leaf senescence [35,36]. Several genes that are involved in ABA synthesis, response, metabolism, receptors or signaling showed differential expression. Two UDP-glucosyltransferase (UGT) 73B1 genes, three UGT73B3 genes, three UGT74D1 genes, one UGT75B1 and one UGT75B2 gene were significantly upregulated during senescence. Meanwhile, two CYP707A3 homologs and one CYP707A4 homolog were upregulated at all leaf senescence stages. However, four homologs of ABA-inactivating UGT71B6 were downregulated during senescence. Investigation of the expressions of ABA signaling related genes showed that many genes were upregulated. For example, among protein phosphatase genes, Arabidopsis thaliana ABA deficient 2 (ABA2), which is implicated in early events of ABA signalling [37], had four homologs that were upregulated. All three ABA receptor pyrabactin resistance (PYR)/PYR1-like (PYL) homologs were upregulated. Of the fourteen homologs of protein phosphatase 2C (PP2C), which negatively regulates the ABA response, most homologs were upregulated during senescence. The ARIA-interacting double AP2 domain protein (ADAP) gene is involved in the ABA response, where it interacts with ARIA, which was also identified in this study. One ADAP homolog was upregulated. Plant U-box 44 (PUB44) genes (also called senescence-associated E3 ubiquitin ligase 1 homologs), which have important roles in the regulation of ABA signaling process, were also significantly upregulated.

Figure 6
figure 6

Hierarchical cluster analysis of differentially expressed transcripts involved in each hormone pathway during senescence. A. abscisic acid, ABA; B. auxin; C. cytokinin, CK; D. brassinosteroid, BR; E. ethylene, ET; F. Gibberellic acid, GA; G. salicylic acid, SA; H. jasmonic acid, JA. The left heat map illustrates the expression profiles for transcripts in the hormone pathway; the right heat map illustrates the expression profiles of transcripts mentioned in the context of hormone. The x-axis refers to each stage of leaf development. The y-axis represents the expression levels of the genes.

Auxin plays important roles in the regulation of root formation, apical dominance, fruit development and senescence [38-40]. For auxin biosynthesis related genes (Figure 6B), one cytochrome P450 75B1 (CYP75B1) homolog and one CYP79B2 homolog were downregulated. P-glycoprotein 1 (PGP1) is associated with auxin transport and is probably involved in the correct polar distribution of auxin during plant development. PGP1 homologs were upregulated during the early senescence stage. The homolog of PGP4, which presumably functions as an auxin importer in root caps [41], appeared to be weakly upregulated during later stages of senescence. In addition, the homologs of putative auxin efflux carrier component 8 (PIN8) from the PIN family (involved in auxin efflux) were initially upregulated at around 25 d and then declined. The crucial factors in the auxin response pathway are the ARF and Aux/IAA proteins, and a previous study showed that Arabidopsis auxin response factor 1 (AtARF1) and AtARF2 were weakly transcriptionally repressed or induced during senescence [34,39]. In our data, the expression of the ARF1 and ARF6 homologs decreased and that of the ARF2 and ARF4 homologs showed increased early and then declined during senescence. In addition, some other auxin response and signal DEGs were downregulated during senescence; e.g., homologs of small auxin upregulated RNA 12 (SAUR12), glycoside hydrolase family 3 (GH3) and auxin signaling f-box 2 (AFB2). These observations indicated a trend of decreasing auxin response at the gene transcription level as a result of changes in signaling proteins.

CKs have long been known as senescence-delaying hormones and are thought to be one of the key signals for the initiation of senescence [42-44]. Among genes related to CK signaling and responses, very few were differentially expressed in our dataset (Figure 6C). Among the CK-degrading oxidase/dehydrogenases [45], the cytokinin oxidase 5 (CKX5) homolog showed a slight decrease. CK levels could also be reduced during senescence by translocation into the vascular system. The phosphorus uptake 1 (purine permease 1, PUP1) homologs, members of the purine TP family, are involved in mediating CK uptake [46]. Most of the PUP1 homologs were upregulated. For the CK signal transduction pathway, the homolog of Arabidopsis histidine phosphotransfer protein 2 (AHP2) was downregulated. Arabidopsis response regulators (ARRs), which are nuclear response regulators, act in the CK signaling pathway [47]. Of five B-type ARRs in our dataset, four ARR12 homologs were downregulated during later stages of senescence. Interestingly, ARR4 and ARR8, which belong to the A-type ARRs, were upregulated. These data suggested that CK might act as a negative regulator during cotton leaf senescence.

BRs play an essential role in several aspects of plant development and previous studies showed that BRs appear to promote developmental senescence [48,49]. Among BRs-biosynthesis genes, two cytochrome P450 72C1 (CYP72C1) homologs and one CYP90A1 homolog were slightly upregulated (Figure 6D). This was consistent with earlier studies in Arabidopsis leaf senescence. For the BR signals, receptors and responses, the majority of altered transcripts were upregulated, including three BRI1-Associated Receptor Kinase 1 (BAK1) homologs, which were involved in BRs signaling, fifteen brassinosteroid insensitive 1 (BRI1) homologs, acts as BRs receptors, and two touch 4 (TCH4) homologs that are associated with BRs responses. BKI1 acts negatively in the BR signaling by inhibiting BRI1’s activity [50,51], which was downregulated in our study.

ET has long been known as a major hormone that promotes leaf senescence [52]. In our study, a large proportion of ET-related key transcripts were upregulated (Figure 6E). In the ET-biosynthesis process, the homologs of 1-aminocyclopropane-1-carboxylic acid oxidase 4 (ACO4) were significantly upregulated at all leaf senescence stages, the 1-aminocyclopropane-1-carboxylic acid synthase 6 (ACS6) was significantly upregulated, while the ACS10 gene was slightly upregulated at the early senescence stage. Among ET signal transduction-related genes, the ET receptor gene ethylene response 1 (ETR1) and ethylene insensitive 2 (EIN2) were not regulated, but several other components of ET signal transduction were significantly upregulated. The ET receptor gene EIN4 that acts upstream of constitutive triple response 1 (CTR1) was continuously upregulated. CTR1, interacts with the ET receptors ETR1 and ethylene response sensor (ERS), was also significantly upregulated. The TP-like RAN1 protein, which might be required to form functional ET receptors [53], and the downstream component EIN3 were also up-regulated during leaf senescence. Meanwhile, several ET signal-related genes were upregulated, including the homologs of EIN3-binding F-box protein 1 (EBF1), ethylene response factor 2 (ERF2), ERF3, ERF4, ERF5, ERF9, ERF12, ERF104, and mitogen-activated protein kinase 4 (MPK4).

GAs comprises a group of growth hormones that control diverse processes, such as cell elongation, reproductive growth and senescence [54]. Previous reports showed that GAs are inhibitors of leaf senescence [55,56]. In our study, a large proportion of GA-related transcripts were downregulated (Figure 6 F), including two gibberellic acid-stimulated in Arabidopsis (GASA1), two GASA14 homologs, and seven alpha-expansin 1 (EXPA1) homologs that involved in GA-signaling, eight aquaporin tip growth defective1 (TIP1) homologs that belong to GA-induced genes and two GA receptors GA-insensitive dwarf 1 (GID1) homologs. In addition, many homologs of GA 2-oxidases and GA 3-oxidases were significantly upregulated during senescence.

SA has been well studied in stress responses, including disease and systemic acquired resistance. It was also reported that the SA has a role in leaf senescence [22,57]. A total of 154 genes associated with various aspects of SA-related pathways were identified (Figure 6G). The majority of transcripts were upregulated, including homologs of syntaxin of plants 121 (SYP121), necrotic spotted lesions 1 (NSL1), enhanced disease susceptibility 1 (EDS1), nonexpressor of PR genes 1 (NPR1), phytoalexin-deficient 4 (PAD4) and NDR1/HIN1-like 25 (NHL25). Differential upregulation was most prominent for PAD4 and NHL25. It was reported that EDS1 and PAD4 interact with each other, and that both genes are transducers of redox signals and important activators of SA signaling [58]. These data indicated that GA and SA might act as modulators of gene expression in senescing leaves.

JA regulates a range of plant growth and development processes, including seed germination, leaf senescence and plant defense against biotic and abiotic stresses [59,60]. Fourteen of 21 JA-biosynthesis genes were upregulated, including the homologs of defective in anther dehiscence1 (DAD1), lipoxygenase 3 (LOX3), allene oxide synthase (AOS), allene oxide cyclase 4 (AOC4), fatty acid desaturase 3 (FAD3) and OPC-8:0 CoA Ligase1 (OPCL1) (Figure 6H). Several genes that are supposedly involved in JA signaling or response were up- or down-regulated during senescence. For example, the homologs of jasmonate resistant 1 (JAR1), Jasmonate-ZIM domain proteins 6 (JAZ6) and coronatine insensitive 1 (COI1) were upregulated, while the homolog of JAZ3 was downregulated. Collectively, these data indicated a coordinated upregulation of JA biosynthesis genes during leaf senescence in cotton.


The utilization of diploid cotton genome is a good supplement for leaf senescence work in cotton

The whole genome sequence of G. hirsutum is currently not completed. The emergence and development of high-throughput sequencing technologies has helped the study of species with complex genomes (such as tetraploids), and have provided an unprecedented opportunity to examine the expression patterns of transcripts in specific tissues of interest. In this study, high-throughput Illumina sequencing technology combined with an advanced mapping strategy was used for expression profiling analysis of different developmental stages in cotton leaves. When mapping the DGE tags, we carefully considered the reference sequences, because the mapping percentages of tags were critical for downstream analysis G. hirsutum is a primary cultivated allotetraploid species and has a tetraploid genome (AD; 2n = 4 × = 52)], whereas, G. raimondii and G. arboreum are diploid cotton species and have D genome and A genome, respectively (D/A; 2n = 2× = 26) [61,62]. Initially, we just used the G. raimondii genome sequence for mapping tag sequences, about 23.6% of all distinct tags could be unambiguously mapped to a gene in the G. raimondii genome sequence, and about 45.2% of all distinct tags mapped to the genome (Additional file 7). It means only small part of tags could be assigned to genes. However, subsequently, we used two databases described in the methods to map the tag sequences, which obviously increased the mapping percentages of tag sequences. Thus, using two databases, but not only the G. raimondii genome sequence, for mapping tag sequences benefitted the downstream analysis.

The initiation phase of leaf senescence in cotton

To obtain an insight into the timing and potential co-regulation of changes in genes and pathways in the complex process of leaf senescence in cotton, it was essential to observe sufficient leaf development stages. In this article, nine time points over 2 months of leaf development from young to full senescence were used. It has been reported that the well-established senescence markers include chlorophyll content, membrane ion leakage and gene expression, and that senescence symptoms usually start from the tip [1,33]. In the experiments reported here, the leaves started to show yellowing at the tip at around 35 d (Figure 1A); the chlorophyll levels reached a maximum at 25 d and then declined (Figure 1B). Chloroplast degeneration was accompanied by chlorophyll degradation and the progressive loss of proteins in the chloroplast, such as ribulose bisphosphate carboxylase (Rubisco) and chlorophyll a/b binding protein (CAB). Thus, leaf senescence is accompanied by decreased expression of genes related to photosynthesis (e.g., CAB) and protein synthesis (e.g., ribosomal protein small subunit (RPS), Rubisco (RBC), and by increased expression of SAGs. Some RBC and CAB genes were downregulated, and WRKY and NAC TFs were upregulated early, at around 25 d (Additional file 8). Considered together, according to the physiological characteristics and the expression pattern of these genes, we hypothesized that the initiation phase of leaf senescence in cotton under our growth condition was between 25 and 35 d.

Comparative analysis of the SAGs with other plants

Previous studies have identified many SAGs during plant senescence. The most recent version of the leaf senescence database (LSD, June, 2013; documented 5,356 SAGs from 44 plant species; the majority (over 3,745) of these were identified in Arabidopsis. These SAGs encode a wide variety of proteins with diverse functions, including macromolecule degradation, nutrient recycling, transcriptional regulation, and hormone response pathways. In our experiment, 3,624 DEGs were identified during developmental senescence in cotton. However, the extent of the similarity or difference of the molecular mechanisms underlying senescence between cotton and other species is unknown. To address this question and analyze SAG orthologs between cotton and other plants, all identified DEGs were compared with those in other species in the leaf senescence database. Based on this database, 2,414 (66.6%) out of 3,624 DEGs have homologous genes in other plants in the LSD, whereas, 1, 210 DEGs generated in this study were not highly homologous to known SAGs in the LSD. This implied that these genes might be specific to cotton during leaf senescence. These novel SAGs were searched against a variety of databases, and approximately 50% of the genes remained of unknown function, which could be considered as putative cotton-specific sequences [63-67]. However, many matched known function genes are involved in photosynthesis (e.g., photosystem I/ II reaction center subunit, light-harvesting complex II protein), chloroplast degeneration (e.g., chlorophyll a/b binding protein), ribosomal protein, metabolic process protein, and zinc finger protein. Therefore, these novel SAGs are candidate target genes for future modelling of senescence transcriptome in cotton.

The TFs are highly enriched during leaf senescence

Several previous studies have highlighted the crucial role of TFs in plant leaf senescence. In the present study, 519 TFs were significantly enriched among the DEGs identified during cotton leaf senescence (Figure 5 and Additional file 5). Among the various TFs, WRKY, bHLH, C3H, NAC, and AP2-EREBP were the top five significantly upregulated TFs during leaf senescence of cotton. WRKY TFs have been reported to be important for senescence [68,69], and are downstream of defense signaling mitogen-activated protein kinase pathways, which are involved in the regulation of JA and SA dependent defense signaling pathways [70,71]. Several AP2-EREBP families modulate responses to leaf senescence-associated signaling molecules, such as ROS, ethylene, JA, ABA, and CK [72,73]. For some of these, there is a striking overlap with the senescence-associated TFs in Arabidopsis: the WRKY, NAC, C3H and AP2-EREBP families. The apparent conservation of the regulatory mechanisms behind the process of senescence across distantly related species is rather encouraging. It suggested that these TFs may lead to fundamental clues about the regulation of senescence. However, there were considerable differences. The expression patterns of these TFs gene were different. For example, some MYB TFs have been shown to participate in leaf senescence, drought stress and disease resistance [74-76], which were significantly upregulated during leaf senescence in our study, but there were no significant changes in Arabidopsis. These data suggested that different regulatory patterns existed in cotton leaf senescence.

Regulation of plant hormone related pathways during senescence

For many years, plant hormones have been thought to play important regulatory roles in either promoting or delaying leaf senescence [77]. For example, ET, ABA and JA are regarded as promoters of senescence, whereas CKs have been observed to delay the progress of senescence. In this study, many genes related to hormone pathways showed differential expression during senescence (Figure 6). The results gave a clear picture with regard to the up or downregulation of genes associated with various hormones. As one of the main plant hormones, ABA has been implicated in the regulation of environmental stresses and senescence [35]. In Arabidopsis, the genes encoding the key enzymes in ABA biosynthesis were upregulated during senescence, such as AtNCED, AtAAO1 and AtAAO3 [12,78]. In our data, key ABA-biosynthesis genes showed no significant changes; only the homolog of AtNCED5 was slightly upregulated at the early senescence stage, around 35 d (Additional file 6). However, many ABA glucosyltransferases, which are involved in ABA catabolism, were significantly upregulated, such as UGT73B3, UGT74D1, UGT75B1 and UGT75B2; while the ABA-inactivating glucosylaseUGT71B6 were downregulated (Figure 6A). In addition, the key ABA signal transduction related genes, including the core group of ABA receptor PYR/PYL family proteins and the type 2C protein phosphatase, were significantly upregulated, which indicated that the ABA signaling pathway is active during leaf senescence. The ABA inducible receptor-like kinase gene of Arabidopsis, RPK1, as a positive regulator of senescence in Arabidopsis Leaves [79], was also significantly upregulated during leaf senescence. Therefore, we predicted that ABA may be regarded as a hormonal trigger of leaf senescence, by maintaining the certain ABA level and upregulating ABA signaling.

Although a delaying effect on senescence by auxin was reported many years ago [35], the function of this hormone in senescence is poorly understood compared with other hormones, such as ABA, ET, JA and CK. Previous studies showed that indole acetic acid (IAA) biosynthetic genes encoding tryptophan synthase (TSA1), IAAld oxidase (AO1), and nitrilases (NIT1-3) are upregulated during age-dependent leaf senescence [12], whereas these genes showed no significant changes in our data (Additional file 6). Many auxin efflux carrier family genes were strongly upregulated, and some genes that play essential roles in auxin signal transduction were significantly differentially expressed (Figure 6B). This suggested that changes in auxin signaling and transport, rather than the endogenous auxin level itself, could be important in modulating senescence. Information on CKs’ involvement in the regulation of senescence in cotton is sparse. A few genes involved in CK transport was upregulated, however, several CK signaling genes, such as the homologs of AHP2 and ARR12, appeared to be downregulated. Thus we presumed that leaf senescence is not so closely related to CK regulation under our growth conditions as other hormones.

Few previous studies have demonstrated the role of GA in leaf senescence. Some bioactive GAs, such as GA1, GA4, and GA53, can be deactivated by GA 2-oxidases, e.g. Gibberellin 2-oxidase 7 (GA2ox7) and GA2ox8 in plants [80-83]. Thus, GA 2-oxidation is a well-characterized type of GA catabolism that is involved in the deactivation of GA. In our data, many homologs of GA 2-oxidases and GA 3-oxidases were significantly upregulated, such as GA2ox1, GA2ox2, GA2ox8, GA3ox1 and GA3ox3 (Figure 6F and Additional file 6). This suggested that during senescence at least some biologically active gibberellins are deactivated. These results are consistent with those obtained from Arabidopsis [12]. The GA signal is received and transduced by the GID1 GA receptor/DELLA repressor pathway [84]. RGL1 and RGL2 belong to the DELLA group of proteins, which are negative regulators of the GA signaling pathway [85,86]. Two homolog of RGL2 were significantly upregulated. This indicated that DELLAs repressed the GA signal. In addition, several other GA signaling proteins were significantly downregulated, such as the homolog of GASA1. Some GA-responsive genes also displayed senescence-associated regulation. A well-known example is the GA-induced aquaporin, TIP1.1 [87], which was downregulated. However, the homologs of GA-responsive xyloglucan endotransglycosylase MERI5 were upregulated during leaf senescence. Studies on TIP1.1 and MERI5 revealed that there is cross talk between the GA and BR pathways. GA and BR regulate AtTIP1.1 antagonistically, whereas MERI5 was is positively and additively regulated by either GA or BR treatment [88]. During leaf senescence, GA activity appears to be reduced; therefore, a different regulator, possibly BR, effects the senescence-associated induction of MERI5.

BRs play an essential role in diverse developmental programs, including leaf senescence, and several evidences have revealed their role in promoting leaf senescence [49]. Some BR biosynthesis-related genes were slightly upregulated, and some others were downregulated (Figure 6D and Additional file 6), suggesting that BR biosynthesis does not significantly increase during leaf senescence. In addition, BIN2, BES1 and BZR1 showed no significant change in expression, which act as the central transmitters of the BR signal from the receptors to the target genes. BRL1 and BRL3 act as BRs receptors, supposedly functioning specifically in provascular differentiation [89]. However, most of the BRL1 and BRL3 homologs were upregulated during leaf senescence, which suggested that BRs might modulate the properties of the vascular tissue in senescing leaves toward efficient export of metabolites.

ET acts as a major hormone in hastening senescence [52]. The ET biosynthetic genes encoding ACC synthase (ACS), ACC oxidase (ACO), and nitrilase are upregulated in senescing leaves [12]. During the process of ethylene metabolism, Arabidopsis ethylene overproducer 1 (ETO1) affected the posttranscriptional regulation of ACS by inhibiting ACS enzyme activity and targeting it for protein degradation. This permits rapid modulation of the concentration of ET [90]. Furthermore, ACS is substrates of MPK3/MPK6, and phosphorylation of ACS2 and ACS6 by MPK3/MPK6 leads to the accumulation of ACS protein, and thus, elevated levels of cellular ACS activity and ethylene production [91,92] In the present study, the homolog of ACS10, which catalyze the rate-limiting step in the ET biosynthesis pathway, was slightly upregulated at the early stage of leaf senescence (Figure 6E and Additional file 6). Coupled with the significant upregulation of the homologs of the MPK3/MPK6 genes, the data suggested that the ACS activity was maintained at a high level during leaf senescence in cotton. In the subsequent step, ACC is converted to ET by ACO, which was also upregulated. Thus, the expression profiles of this set of proteins imply an increase in ET production in senescing leaves. ET signal transduction supposedly follows a “linear” pathway, with membrane-bound receptors at the beginning, multiple positive and negative regulators in between, and TFs at the end of the chain. The expression of ET receptor ETR1, ETR2 and ERS1 were not significantly regulated. EIN3/EIN4, as the downstream components of the ET receptor showed increased expression. Meanwhile, one immediate target of EIN3 is ERF1 [93], which appeared to be upregulated. Several other ERF proteins were also upregulated but with different expression profiles. Moreover, three EIN3-binding F box proteins were also upregulated. Taken together, the data suggested that ET may act as a positive regulator in leaf senescence in cotton.

SA is involved in pathogen response and pathogen-mediated cell death. It was also reported that the SA-signaling pathway is active in the control of gene expression during developmental senescence [57]. In our study, many SA signaling genes were significantly upregulated. For example, two homologs of NHL25, one homolog of NPR1 and two homologs of PAD4 were significantly upregulated (Figure 6G). Scarecrow-like 14 (SCL14), a member of the GRAS family of regulatory proteins, interacts with the TGA2 transcription factor and affects the transcription of stress-responsive genes [94]. We observed that all four homologs of this gene were upregulated. In addition, the expression profiles of many SA-responsive genes were also upregulated, such as the homolog of the AF2 gene, accelerated cell death 6 (ACD6). Moreover, as mentioned before, many WRKY TFs involved in the regulation of SA-dependent defense signaling pathways were upregulated, which revealed that SA has important roles in positively regulating leaf senescence in cotton by interacting with the WRKY family.

JA functions in the induction of leaf senescence in many plant species [60,95] by regulating the expression of various SAGs [22]. Here, JA biosynthesis genes showed a significantly increase in expression almost the entire process of senescence, including LOX3, AOS, AOC4 and OPR3. The COI1, an F-box protein, as a key component of the receptor complex in the JA signaling pathway [96,97], were also up-regulated during senescence. This protein forms part of an enzyme complex called SCFCOI1. JAZ proteins are ubiquitinated via SCFCOI1 in response to JA, they function as repressors of the JA signaling pathway, and recent investigations revealed that COI1 (in the SCFCOI1 complex) and several JAZ proteins form a coreceptor complex [98,99]. The genes implicated in controlling JA responses are upregulated, including JAZ6, JAZ7 and MYC2. JAZ proteins have been shown to bind to MYC2 directly and prevent its action [100]. Thus, it indicated that the JA signaling was active at the very early stage during leaf senescence, and these data also provide additional evidence for JA as a modulator of gene expression in developmental leaf senescence.

Transcriptome dynamics during leaf senescence

As expected, we identified groups of genes that are co-expressed and hence may be co-regulated across the entire time series; however, it is very important to identify an ordering of events across the leaf developmental stages. It is clear that there are progressive changes in genes being upregulated or downregulated as senescence progresses, and these are highly informative in indicating changes in metabolic pathways. Genes upregulated from the first time point (15 d), before the leaf is fully expanded, were enriched for response to biotic stimulus, response to endogenous stimulus, response to abiotic stimulus, response to stress. In addition, genes involved in both JA and SA responses were enhanced in expression before 25 d, indicating that there may be a role for these hormones even before the leaf is fully expanded. JA-related genes upregulated at this time point include genes required for JA biosynthesis such as two LOX, two AOC, OPR, and OPCL (Additional file 6). Genes implicated in JA responses and signalling were upregulated, including JAZ6, JAZ7, and JAR1. SA biosynthetic genes showed no significantly increase in senescence, only some SA metabolism genes such as UDP-glycosyltransferase 74F, UGT74F1 and UGT74F2 were regulated. The SA glycosides are actively transported from the cytosol to the vacuole as an inactive storage precursor that can later be converted back to SA [101]. Besides, several JA-signaling genes (e.g., NHL25 and PAD4) were upregulated before 25 d. This suggested that these two hormones modulated leaf senescence primarily by responding to environmental cues and stress at the early stage of senescence. Tetrapyrroles are essential components of critical biological processes, including respiration and photosynthesis; many tetrapyrrole or chlorophyll biosynthesis genes are first downregulated at 25 d. Thus, the requirement for de novo chlorophyll biosynthesis appears to cease in stage of leaf maturity, indicating that all chloroplasts are fully developed. Genes involved in ribosome and ribosomal protein genes were downregulated at 25 d, indicating that expansion of the ribosomal content of the cells has slowed down. Also, downregulated genes were significantly enriched for genes linked to the plastid, vacuole and thylakoid, and with functions in metabolic processes, particularly photosynthesis, lipid biosynthetic process and fatty acid biosynthetic process at this stage. It indicated that photosynthetic activity probably started to drop at this stage. Down-regulation of these groups of genes reflected the shut down of cellular biosynthetic activity at the early senescence stage. Moreover, at 25 d, within the molecular function annotations, transferase activity, catalytic activity, kinase activity, oxygen binding, calmodulin binding and transcription regulator activity were overrepresented, whereas there was significant enrichment for protein metabolic process and protein modification process in the biological process terms, reflecting kinase-signalling cascades are likely to function and the considerable degradation that is underway.

At the next time point (35 d), there is an overrepresentation of genes involved in macromolecule metabolic process, response to ABA stimulus, regulation of defense response. Similarly, genes involved in ribonucleotide binding, ATP binding, protein binding and lipase activity are up regulated. The previous investigations showed that ABA is a key hormone mediating plant responses to environmental stresses, and ABA induces the expression of antioxidant genes and enhances the activities of antioxidative enzymes that play a role in protecting the cellular functions [102]. ABA-related genes upregulated at this time point include genes required for ABA signal transduction such as PP2C and RPK1. It appeared that ABA may control activities of both the cellular protection activities and senescence activities at the early stages of senescence. By the later stage of senescence, many genes encoding TFs are obviously upregulated, including WRKY, NAC and C3H family. Genes for response to chemical stimulus also upregulated, some genes annotated ET related pathway, i.e., ACS6, ACO4, ERF12 and ERF104 being upregulated late, indicating that the ET regulation of senescence may have a significant role at this time. Once the senescence process is initiated, massive degradation and recycling of macromolecule occurs during this stage. The leaf was converted from a source of photosynthetic carbon to mobilizable nutrients, as cellular components become degraded and mobilized from the leaf. These results are of great importance for future modeling of senescence transcriptional networks. However, the exact role of these regulators and the cross-talk between them to fine-tune the leaf senescence in cotton remains undiscovered.


To the best of our knowledge, this study is the first to reveal genome-wide multiple time-course transcriptome dynamics during leaf senescence in cotton. In total, 3,624 DEGs were identified and analyzed for their potential role in cotton leaf development using clustering, GO enrichment and KOBAS analysis. Furthermore, 519 differentially expressed TF transcripts from 50 TF families were identified; some of these TFs have been reported to be critical components of plant adaptive response to biotic, abiotic stresses and senescence. Eight hormone-related differentially expressed genes were identified, of which many genes involved in the ET, BR, JA, SA and ABA pathways were upregulated, indicating that these hormone-related genes might act as positive regulators in cotton leaf development and senescence. However, most auxin and GA pathway-related genes were downregulated, suggesting that these two hormones may act as negative regulators of senescence. Thus, our study revealed a large set of candidate genes that may play important roles in leaf development and senescence in cotton. Meanwhile, the detailed investigations of the pathways and candidate genes identified in this study would not only be useful for enhancing our understanding of the molecular mechanism involved in leaf development and senescence of G. hirsutum, but also may provide potential targets for crop improvement via breeding and reverse genetics approaches.


Plant material preparation and RNA extraction

G. hirsutum cv. CCRI 36 was grown in a cotton climatic chamber at 30/25 ± 3°C under a light/dark cycle of 12/12 h. Similar sized emerging fourth true leaves were tagged. The day when emerging true leaf 4 was expanded was considered the first day. Sampling of leaves started on the 5th day and continued every 10 days until full senescence was reached. At each time point, three biological replicates were designed; the samples of each biological replicate were pooled from 10 plants, the plants being randomly selected to avoid any potential effects of position within the growth room. The three biological replicates were mixed divided into two parts, one was used for RNA isolation for Illumina sequencing (only one mixed sample was used in each developmental stage in sequencing.), and the other was used for chlorophyll and malondialdehyde (MDA) measurements. An RNA prep Pure Plant Kit (TIANGEN BIOTECH) was used to isolate the total RNA, according to the manufacturer's protocol.

MDA and chlorophyll measurements

Chlorophyll was extracted using 80% (v/v) acetone, vortexed, and then stored at 4°C overnight in the dark. The samples were then centrifuged for 2 min at 12,500 × g, and the absorbance of 1 mL of supernatant was measured at 663 and 645 nm using the method described by Lichtenthaler [103]. MDA contents were measured by the 5% thiobarbituric acid reaction, as described by the procedure of Saher [104].

DGE sequencing and analysis

Approximately, 8 μg of total RNA was used. Tag libraries were prepared using the Illumina Gene Expression Sample Prep Kit, following the manufacturer’s protocol, as described by Luan et al. [105]. The libraries were then sequenced using an Illumina HiSeq™ 2000 with 50-bp single-end (SE) reads each. To increase the mapping percentages of tag sequences, two different reference databases were employed, which included G. raimondii genome and G. arboretum genome. In order to remove redundancy, these two genome sequences were merged into one reference database according to following procedure: all CATG + 17-nucleotide tag sequences from the 3’-most CATG to the poly(A)-tail of the predicted coding sequences from the genome data of G. raimondii and G. arboretum were created, all these tags were allowed no more than one nucleotide mismatch to remove redundancy; if two or more tags from different genome sequences are the same, G. arboretum genome sequence were selected with higher priority. The remaining genome sequences were used as reference sequence for DGE tag mapping. Then, the preprocessed database of all possible CATG + 17-nucleotide tag sequences was generated using reference sequence. All tags were allowed no more than one nucleotide mismatch. Clean tags mapping to reference sequences from multiple genes were filtered out, and the remaining clean tags were designated as unambiguous clean tags. For gene expression analysis, the number of unambiguous clean tags for each gene was calculated and then normalized to TPM (number of transcripts per million clean tags) [106,107].

Validation of gene expression profile using qRT-PCR

To validate the results of the DGE-based analyses, the expressions of 25 genes that were randomly selected based on different expression levels in DGE-based data were measured by qRT-PCR. Importantly, the samples used for RNA isolation in qRT-PCR experiments were different from the samples used in RNA-Seq analysis. The RNA samples were prepared using the same method mentioned above. cDNA was reverse transcribed from RNA by PrimeScript® RT reagent Kit with gDNA Eraser (Takara, Biotechnology Co., Ltd.) using an Oligo dT Primer and random hexamer as RT primers, according to the manufacturer’s protocol. The specific primers for the selected genes and internal control gene (actin) are listed in Additional file 9. Samples were run in technical triplicate on each plate using the SYBR Green PCR Master Mix (Takara) on an ABI 7500 Real-time PCR System (Applied Biosystems, USA) following the manufacturer’s instructions. The relative expression levels were calculated by the comparative 2Ct method [108]. The results were normalized to the expression level of actin and relative to the 15 d sample. The Pearson correlation coefficients of the expression patterns of selected genes between qRT-PCR and DGE were calculated using the SAS software.

Identification and analysis of DEGs

To identify DEGs across the six different time points, pair wise comparisons among the six samples were performed using a rigorous algorithm method based on a previous method [109]. The DEGs were obtained after filtering using a threshold false discovery rate (FDR) of ≤0.001 and an absolute value of log2Ratio ≥2. Two-fold change has been generally accepted as a significant or meaningful change in DGE data analyses; however, we used four-fold change as the cut-off for more rigorous analysis in this study. Additionally, one of the two samples was required to have a raw intensity of gene (raw mapped tags) ≥10 in the pairwise comparison. MultiExperiment Viewer (MeV, v4.7.4) [110] produced the cluster analysis of the DEGs, based on the K-means method [111]. MeV performed the clustering analysis of co-regulated genes expression patterns based on hierarchical clustering. For GO annotation enrichment, the best hit corresponding to each cotton transcript was identified using a BLAST search. The BiNGO tool preformed the GO enrichment analysis of gene clusters, as described by Maere [30]. KOBAS was then used to identify biochemical pathways involved in cotton leaf senescence and to calculate the statistical significance of each pathway with the default criteria [112]. KOBAS uses Fisher’s exact test to identify significant pathways and then performs an FDR correction to reduce type-1 errors. Pathways with corrected P values <0.05 were considered statistically significant.

To identify genes encoding for TFs among the DEGs, all DEGs were compared with protein sequences download from a comprehensive plant TF database, PlnTFDB [113] using Blastx with an E-value ≤ 10−10. The HMMER (E-value ≤ 10−10, program was used to search against the BLAST search results using the Pfam program ( Finally, we manually checked the search results to reduce false positives. For analysis and classification of hormone-related genes, all DEGs were compared with the protein sequences from the Arabidopsis Hormone Database ( using Blastx (E-value ≤ 10−10).

Availability of supporting data

Raw Illumina reads have been deposited into NCBI’s SRA (sequence read archive), the accession number is SRR656612.



Senescence-associated genes


Transcription factors


Bacillus thuringiensis




Abscisic acid




Digital gene expression


Differentially expressed genes


Leaf senescence database


MultiExperiment Viewer


Kyoto Encyclopedia of Genes and Genomes


KEGG Orthology Based Annotation System






Jasmonic acid


Salicylic acid




ABA deficient 2


Pyrabactin resistance




Protein phosphatase 2C


ARIA-interacting double AP2 domain protein


Plant U-box 44


Cytochrome P450 75B1


P-glycoprotein 1


Putative auxin efflux carrier component 8


Auxin response factor 1


Small auxin upregulated RNA 12


Glycoside hydrolase family 3


Auxin signaling f-box 2


Cytokinin oxidase 5


Purine permease 1


Arabidopsis histidine phosphotransfer protein 2


Arabidopsis response regulators


Brassinosteroid insensitive 1


BRI1-Associated Receptor Kinase 1


Two touch 4


1-aminocyclopropane-1-carboxylic acid oxidase 4


1-aminocyclopropane-1-carboxylic acid synthase 6


Ethylene response 1


Ethylene insensitive 2


Constitutive triple response 1


Ethylene response sensor


EIN3-binding F-box protein 1


Ethylene response factor 2


Mitogen-activated protein kinase 4


Gibberellic acid-stimulated in Arabidopsis


Alpha-expansin 1


Tip growth defective1


GA-insensitive dwarf 1


Syntaxin of plants 121


Necrotic spotted lesions 1


Enhanced disease susceptibility 1


Nonexpressor of PR genes 1


Phytoalexin-deficient 4


NDR1/HIN1-like 25


Defective in anther dehiscence1


Lipoxygenase 3


Allene oxide synthase


Allene oxide cyclase 4


Fatty acid desaturase 3


OPC-8:0 CoA Ligase1


Jasmonate resistant 1


Jasmonate-ZIM domain proteins 6


Coronatine insensitive 1


Ribulose bisphosphate carboxylase


Chlorophyll a/b binding protein


Ribosomal protein small subunit




Indole acetic acid


Tryptophan synthase


IAAld oxidase




Gibberellin 2-oxidase 7


Ethylene overproducer 1


EIN3-binding F box proteins


Scarecrow-like 14


Accelerated cell death 6


  1. Lim PO, Kim HJ, Gil Nam H. Leaf Senescence. Annu Rev Plant Biol. 2007;58(1):115–36.

    CAS  PubMed  Google Scholar 

  2. Gan S, Amasino RM. Making sense of senescence (molecular Genetic regulation and manipulation of leaf senescence). Plant Physiol. 1997;113(2):313–9.

    PubMed Central  CAS  PubMed  Google Scholar 

  3. Nam HG. The molecular genetic analysis of leaf senescence. Curr Opin Biotech. 1997;8(2):200–7.

    CAS  PubMed  Google Scholar 

  4. Hortensteiner S, Feller U. Nitrogen metabolism and remobilization during senescence. J Exp Bot. 2002;53(370):927–37.

    CAS  PubMed  Google Scholar 

  5. Desclos M, Etienne P, Coquet L, Jouenne T, Bonnefoy J, Segura R, et al. A combined 15 N tracing/proteomics study in Brassica napus reveals the chronology of proteomics events associated with N remobilisation during leaf senescence induced by nitrate limitation or starvation. Proteomics. 2009;9(13):3580–608.

    CAS  PubMed  Google Scholar 

  6. De Michele R, Formentin E, Todesco M, Toppo S, Carimi F, Zottini M, et al. Transcriptome analysis of Medicago truncatula leaf senescence: similarities and differences in metabolic and transcriptional regulations as compared with Arabidopsis, nodule senescence and nitric oxide signalling. New Phytol. 2009;181(3):563–75.

    PubMed  Google Scholar 

  7. Liu L, Zhou Y, Zhou G, Ye R, Zhao L, Li X, et al. Identification of early senescence-associated genes in rice flag leaves. Plant Mol Biol. 2008;67(1–2):37–55.

    CAS  PubMed  Google Scholar 

  8. Rivero RM, Kojima M, Gepstein A, Sakakibara H, Mittler R, Gepstein S, et al. Delayed leaf senescence induces extreme drought tolerance in a flowering plant. Proc Natl Acad Sci U S A. 2007;104(49):19631–6.

    PubMed Central  CAS  PubMed  Google Scholar 

  9. Guo Y, Gan S. Leaf senescence: signals, execution, and regulation. Curr Top Dev Biol. 2005;71:83–112.

    CAS  PubMed  Google Scholar 

  10. Lim PO, Nam HG. The molecular and genetic control of leaf senescence and longevity in Arabidopsis. Curr Top Dev Biol. 2005;67:49–83.

    CAS  PubMed  Google Scholar 

  11. Balazadeh S, Riano-Pachon DM, Mueller-Roeber B. Transcription factors regulating leaf senescence in Arabidopsis thaliana. Plant Biol (Stuttg). 2008;10 Suppl 1:63–75.

    Google Scholar 

  12. van der Graaff E. Transcription analysis of Arabidopsis membrane transporters and hormone pathways during developmental and induced leaf senescence. Plant Physiol. 2006;141(2):776–92.

    PubMed Central  PubMed  Google Scholar 

  13. Guinn G. Abscisic acid and cutout in cotton. Plant Physiol. 1985;77(1):16–20.

    PubMed Central  CAS  PubMed  Google Scholar 

  14. Wright PR. Premature senescence of cotton (gossypium hirsutum L.) – predominantly a potassium disorder caused by an imbalance of source and sink. Plant Soil. 1999;211(2):231–9.

    CAS  Google Scholar 

  15. Wright PR. Research into early senescence syndrome in cotton. Better Crops Int. 1998;2(12):14–6.

    Google Scholar 

  16. Dong H, Niu Y, Li W, Zhang D. Effects of cotton rootstock on endogenous cytokinins and abscisic acid in xylem sap and leaves in relation to leaf senescence. J Exp Bot. 2008;59(6):1295–304.

    CAS  PubMed  Google Scholar 

  17. Radin JW. Water relations of cotton plants under nitrogen deficiency. Physiol Plantarum. 1981;51(1):145–9.

    CAS  Google Scholar 

  18. Dong H, Li W, Eneji AE, Zhang D. Nitrogen rate and plant density effects on yield and late-season leaf senescence of cotton raised on a saline field. Field Crop Res. 2012;126:137–44.

    Google Scholar 

  19. Ganieva RA, Allahverdiyev SR, Guseinova NB, Kavakli HI, Nafisi S. Effect of salt stress and synthetic hormone polystimuline K on the photosynthetic activity of cotton (Gossypium hirsutum). Turk J Bot. 1998;4(22):217–21.

    Google Scholar 

  20. Chang CW. Carbon dioxide and senescence in cotton plants. Plant Physiol. 1975;55(3):515–9.

    PubMed Central  CAS  PubMed  Google Scholar 

  21. Kakani VG, Reddy KR, Zhao D, Gao W. Senescence and hyperspectral reflectance of cotton leaves exposed to ultraviolet-B radiation and carbon dioxide. Physiol Plantarum. 2004;121(2):250–7.

    CAS  Google Scholar 

  22. Buchanan-Wollaston V, Page T, Harrison E, Breeze E, Lim PO, Nam HG, et al. Comparative transcriptome analysis reveals significant differences in gene expression and signalling pathways between developmental and dark/starvation-induced senescence in Arabidopsis. Plant J. 2005;42(4):567–85.

    CAS  PubMed  Google Scholar 

  23. Andersson A, Keskitalo J, Sjodin A, Bhalerao R, Sterky F, Wissel K, et al. A transcriptional timetable of autumn senescence. Genome Biol. 2004;5(4):R24.

    PubMed Central  PubMed  Google Scholar 

  24. Breeze E, Harrison E, McHattie S, Hughes L, Hickman R, Hill C, et al. High-resolution temporal profiling of transcripts during Arabidopsis leaf senescence reveals a distinct chronology of processes and regulation. Plant Cell. 2011;23(3):873–94.

    PubMed Central  CAS  PubMed  Google Scholar 

  25. Ansorge WJ. Next-generation DNA sequencing techniques. N Biotechnol. 2009;25(4):195–203.

    CAS  PubMed  Google Scholar 

  26. Lister R, Gregory BD, Ecker JR. Next is now: new technologies for sequencing of genomes, transcriptomes, and beyond. Curr Opin Plant Biol. 2009;12(2):107–18.

    PubMed Central  CAS  PubMed  Google Scholar 

  27. Tang Q, Ma X, Mo C, Wilson IW, Song C, Zhao H, et al. An efficient approach to finding Siraitia grosvenorii triterpene biosynthetic genes by RNA-seq and digital gene expression analysis. BMC Genomics. 2011;12(1):343.

    PubMed Central  CAS  PubMed  Google Scholar 

  28. Hao DC, Ge G, Xiao P, Zhang Y, Yang L. The first insight into the tissue specific taxus transcriptome via illumina second generation sequencing. PLoS One. 2011;6(6):e21220.

    PubMed Central  CAS  Google Scholar 

  29. Wang Z, Fang B, Chen J, Zhang X, Luo Z, Huang L, et al. De novo assembly and characterization of root transcriptome using Illumina paired-end sequencing and development of cSSR markers in sweetpotato (Ipomoea batatas). BMC Genomics. 2010;11(1):726.

    PubMed Central  CAS  PubMed  Google Scholar 

  30. Maere S, Heymans K, Kuiper M. BiNGO: a cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics. 2005;21(16):3448–9.

    CAS  PubMed  Google Scholar 

  31. Hickman R, Hill C, Penfold CA, Breeze E, Bowden L, Moore JD, et al. A local regulatory network around three NAC transcription factors in stress responses and senescence in Arabidopsis leaves. Plant J. 2013;75(1):26–39.

    PubMed Central  CAS  PubMed  Google Scholar 

  32. Gregersen PL, Holm PB. Transcriptome analysis of senescence in the flag leaf of wheat (Triticum aestivum L.). Plant Biotechnol J. 2007;5(1):192–206.

    CAS  PubMed  Google Scholar 

  33. Guo Y, Gan S. AtNAP, a NAC family transcription factor, has an important role in leaf senescence. Plant J. 2006;46(4):601–12.

    CAS  PubMed  Google Scholar 

  34. Lin J, Wu S. Molecular events in senescing. Plant J. 2004;39(4):612–28.

    CAS  PubMed  Google Scholar 

  35. Noodén LD. Abscisic Acid, Auxin, and Other Regulators of Senescence. In: Noodén LD, Leopold A, editors. Senescence and Aging in Plants. San Diego: Academic; 1988. p. 329–68.

    Google Scholar 

  36. Hung KT, Kao CH. Hydrogen peroxide is necessary for abscisic acid-induced senescence of rice leaves. J Plant Physiol. 2004;161(12):1347–57.

    CAS  PubMed  Google Scholar 

  37. Leung J. The Arabidopsis ABSCISIC ACID-INSENSITIVE2 (ABI2) and ABI1 genes encode homologous protein phosphatases 2C involved in abscisic acid signal transduction. Plant Cell. 1997;9(5):759–71.

    PubMed Central  CAS  PubMed  Google Scholar 

  38. Kim JI, Murphy AS, Baek D, Lee SW, Yun DJ, Bressan RA, et al. YUCCA6 over-expression demonstrates auxin function in delaying leaf senescence in Arabidopsis thaliana. J Exp Bot. 2011;62(11):3981–92.

    PubMed Central  CAS  PubMed  Google Scholar 

  39. Ellis CM. AUXIN RESPONSE FACTOR1 and AUXIN RESPONSE FACTOR2 regulate senescence and floral organ abscission in Arabidopsis thaliana. Development. 2005;132(20):4563–74.

    CAS  PubMed  Google Scholar 

  40. Woodward AW, Bartel B. Auxin: regulation, action, and interaction. Ann Bot. 2005;95(5):707–35.

    PubMed Central  CAS  PubMed  Google Scholar 

  41. Terasaka K. PGP4, an ATP binding cassette P-glycoprotein, catalyzes auxin transport in Arabidopsis thaliana roots. Plant Cell. 2005;17(11):2922–39.

    PubMed Central  CAS  PubMed  Google Scholar 

  42. McCabe MS, Garratt LC, Schepers F, Jordi WJ, Stoopen GM, Davelaar E, et al. Effects of P(SAG12)-IPT gene expression on development and senescence in transgenic lettuce. Plant Physiol. 2001;127(2):505–16.

    PubMed Central  CAS  PubMed  Google Scholar 

  43. Gan S, Amasino RM. Inhibition of leaf senescence by autoregulated production of cytokinin. Science. 1995;270(5244):1986.

    CAS  PubMed  Google Scholar 

  44. Richmond AE, Lang A. Effect of kinetin on protein content and survival of detached xanthium leaves. Science. 1957;125(3249):650–1.

    CAS  Google Scholar 

  45. Schmülling T, Werner T, Riefler M, Krupková E, Bartrina Y, Manns I. Structure and function of cytokinin oxidase/dehydrogenase genes of maize, rice, Arabidopsis and other species. J Plant Res. 2003;116(3):241–52.

    PubMed  Google Scholar 

  46. Cedzich A, Stransky H, Schulz B, Frommer WB. Characterization of cytokinin and adenine transport in Arabidopsis cell cultures. Plant Physiol. 2008;148(4):1857–67.

    PubMed Central  CAS  PubMed  Google Scholar 

  47. Hill K, Mathews DE, Kim HJ, Street IH, Wildes SL, Chiang YH, et al. Functional characterization of type-B response regulators in the Arabidopsis cytokinin response. Plant Physiol. 2013;162(1):212–24.

    PubMed Central  CAS  PubMed  Google Scholar 

  48. Yin Y, Wang ZY, Mora-Garcia S, Li J, Yoshida S, Asami T, et al. BES1 accumulates in the nucleus in response to brassinosteroids to regulate gene expression and promote stem elongation. Cell. 2002;109(2):181–91.

    CAS  PubMed  Google Scholar 

  49. Clouse SD, Sasse JM. BRASSINOSTEROIDS: essential regulators of plant growth and development. Annu Rev Plant Physiol Plant Mol Biol. 1998;49(1):427–51.

    CAS  PubMed  Google Scholar 

  50. Li J, Jin H. Regulation of brassinosteroid signaling. Trends Plant Sci. 2007;12(1):37–41.

    CAS  PubMed  Google Scholar 

  51. Wang X, Chory J. Brassinosteroids regulate dissociation of BKI1, a negative regulator of BRI1 signaling, from the plasma membrane. Science. 2006;313(5790):1118–22.

    CAS  PubMed  Google Scholar 

  52. Abeles FB, Dunn LJ, Morgens P, Callahan A, Dinterman RE, Schmidt J. Induction of 33-kD and 60-kD peroxidases during ethylene-induced senescence of cucumber cotyledons. Plant Physiol. 1988;87(3):609–15.

    PubMed Central  CAS  PubMed  Google Scholar 

  53. Woeste KE, Kieber JJ. A strong loss-of-function mutation in RAN1 results in constitutive activation of the ethylene response pathway as well as a rosette-lethal phenotype. Plant Cell. 2000;12(3):443–55.

    PubMed Central  CAS  PubMed  Google Scholar 

  54. Rodrigues C, Vandenberghe LP, de Oliveira J, Soccol CR. New perspectives of gibberellic acid production: a review. Crit Rev Biotechnol. 2012;32(3):263–73.

    CAS  PubMed  Google Scholar 

  55. Yu K, Wang Y, Wei J, Ma Q, Yu D, Li J. Improving rhizome yield and quality of Paris polyphylla through gibberellic acid-induced retardation of senescence of aerial parts. Plant Signal Behav. 2009;4(5):413–5.

    PubMed Central  CAS  PubMed  Google Scholar 

  56. Schippers JHM, Jing H, Hille J, Dijkwel PP. Developmental and Hormonal Control of Leaf Senescence. Oxford: Blackwell Publishing Ltd; 2007. p. 145–70.

    Google Scholar 

  57. Morris K, MacKerness SA, Page T, John CF, Murphy AM, Carr JP, et al. Salicylic acid has a role in regulating gene expression during leaf senescence. Plant J. 2000;23(5):677–85.

    CAS  PubMed  Google Scholar 

  58. Wiermer M, Feys BJ, Parker JE. Plant immunity: the EDS1 regulatory node. Curr Opin Plant Biol. 2005;8(4):383–9.

    CAS  PubMed  Google Scholar 

  59. Wasternack C. Jasmonates: an update on biosynthesis, signal transduction and action in plant stress response, growth and development. Ann Bot. 2007;100(4):681–97.

    PubMed Central  CAS  PubMed  Google Scholar 

  60. He Y, Fukushige H, Hildebrand DF, Gan S. Evidence supporting a role of jasmonic acid in Arabidopsis leaf senescence. Plant Physiol. 2002;128(3):876–84.

    PubMed Central  CAS  PubMed  Google Scholar 

  61. Mei M, Syed NH, Gao W, Thaxton PM, Smith CW, Stelly DM, et al. Genetic mapping and QTL analysis of fiber-related traits in cotton (Gossypium). Theor Appl Genet. 2004;108(2):280–91.

    CAS  PubMed  Google Scholar 

  62. Han ZG, Guo WZ, Song XL, Zhang TZ. Genetic mapping of EST-derived microsatellites from the diploid Gossypium arboreum in allotetraploid cotton. Mol Genet Genomics. 2004;272(3):308–27.

    CAS  PubMed  Google Scholar 

  63. Samuel YS, Cheung F, Lee JJ, Ha M, Wei NE, Sze SH, et al. Accumulation of genome-specific transcripts, transcription factors and phytohormonal regulators during early stages of fiber cell development in allotetraploid cotton. Plant J. 2006;47(5):761–75.

    Google Scholar 

  64. Yuan D, Tu L, Zhang X. Generation, annotation and analysis of first large-scale expressed sequence tags from developing fiber of Gossypium barbadense L. PLoS One. 2011;6(7):e22758.

    PubMed Central  CAS  PubMed  Google Scholar 

  65. Shi YH, Zhu SW, Mao XZ, Feng JX, Qin YM, Zhang L, et al. Transcriptome profiling, molecular biological, and physiological studies reveal a major role for ethylene in cotton fiber cell elongation. Plant Cell. 2006;18(3):651–64.

    PubMed Central  CAS  PubMed  Google Scholar 

  66. Yoo MJ, Wendel JF. Comparative evolutionary and developmental dynamics of the cotton (Gossypium hirsutum) fiber transcriptome. PLoS Genet. 2014;10(1):e1004073.

    PubMed Central  PubMed  Google Scholar 

  67. Nigam D, Kavita P, Tripathi RK, Ranjan A, Goel R, Asif M, et al. Transcriptome dynamics during fibre development in contrasting genotypes of Gossypium hirsutum L. Plant Biotechnol J. 2014;12(2):204–18.

    CAS  PubMed  Google Scholar 

  68. Miao Y, Laun T, Zimmermann P, Zentgraf U. Targets of the WRKY53 transcription factor and its role during leaf senescence in Arabidopsis. Plant Mol Biol. 2004;55(6):853–67.

    CAS  PubMed  Google Scholar 

  69. Robatzek S, Somssich IE. A new member of the Arabidopsis WRKY transcription factor family, AtWRKY6, is associated with both senescence- and defence-related processes. Plant J. 2001;28(2):123–33.

    CAS  PubMed  Google Scholar 

  70. Eulgem T, Somssich IE. Networks of WRKY transcription factors in defense signaling. Curr Opin Plant Biol. 2007;10(4):366–71.

    CAS  PubMed  Google Scholar 

  71. Ulker B, Somssich IE. WRKY transcription factors: from DNA binding towards biological function. Curr Opin Plant Biol. 2004;7(5):491–8.

    PubMed  Google Scholar 

  72. Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. AP2/ERF family transcription factors in plant abiotic stress responses. Biochim Biophys Acta. 2012;1819(2):86–96.

    CAS  PubMed  Google Scholar 

  73. Nakano T, Suzuki K, Fujimura T, Shinshi H. Genome-wide analysis of the ERF gene family in Arabidopsis and rice. Plant Physiol. 2006;140(2):411–32.

    PubMed Central  CAS  PubMed  Google Scholar 

  74. Lin M, Lai D, Pang C, Fan S, Song M, Yu S. Generation and analysis of a large-scale expressed sequence Tag database from a full-length enriched cDNA library of developing leaves of Gossypium hirsutum L. PLoS One. 2013;8(10):e76443.

    PubMed Central  CAS  PubMed  Google Scholar 

  75. Dai X, Xu Y, Ma Q, Xu W, Wang T, Xue Y, et al. Overexpression of an R1R2R3 MYB gene, OsMYB3R-2, increases tolerance to freezing, drought, and salt stress in transgenic Arabidopsis. Plant Physiol. 2007;143(4):1739–51.

    PubMed Central  CAS  PubMed  Google Scholar 

  76. Liu Y, Schiff M, Dinesh-Kumar SP. Involvement of MEK1 MAPKK, NTF6 MAPK, WRKY/MYB transcription factors, COI1 and CTR1 in N-mediated resistance to tobacco mosaic virus. Plant J. 2004;38(5):800–9.

    CAS  PubMed  Google Scholar 

  77. Thomas H, Stoddart JL. Leaf senescence. Annual Rev Plant Physiol. 1980;31:83–111.

    CAS  Google Scholar 

  78. Seo M, Peeters AJM, Koiwai H, Oritani T, Marion-Poll A, Zeevaart JAD, et al. The Arabidopsis aldehyde oxidase 3 (AAO3) gene product catalyzes the final step in abscisic acid biosynthesis in leaves. Proc Natl Acad Sci U S A. 2000;97(23):12908–13.

    PubMed Central  CAS  PubMed  Google Scholar 

  79. Lee IC, Hong SW, Whang SS, Lim PO, Nam HG, Koo JC. Age-dependent action of an ABA-inducible receptor kinase, RPK1, as a positive regulator of senescence in Arabidopsis leaves. Plant Cell Physiol. 2011;52(4):651–62.

    CAS  PubMed  Google Scholar 

  80. Lee DJ, Zeevaart JA. Molecular cloning of GA 2-oxidase3 from spinach and its ectopic expression in Nicotiana sylvestris. Plant Physiol. 2005;138(1):243–54.

    PubMed Central  CAS  PubMed  Google Scholar 

  81. Schomburg FM, Bizzell CM, Lee DJ, Zeevaart JA, Amasino RM. Overexpression of a novel class of gibberellin 2-oxidases decreases gibberellin levels and creates dwarf plants. Plant Cell. 2003;15(1):151–63.

    PubMed Central  CAS  PubMed  Google Scholar 

  82. Yamaguchi SKY. Gibberellin biosynthesis: its regulation by endogenous and environmental signals. Plant Cell Physiol. 2000;41(3):251–7.

    CAS  PubMed  Google Scholar 

  83. Thomas SG, Phillips AL, Hedden P. Molecular cloning and functional expression of gibberellin 2- oxidases, multifunctional enzymes involved in gibberellin deactivation. Proc Natl Acad Sci U S A. 1999;96(8):4698–703.

    PubMed Central  CAS  PubMed  Google Scholar 

  84. Ueguchi-Tanaka M, Nakajima M, Motoyuki A, Matsuoka M. Gibberellin receptor and its role in gibberellin signaling in plants. Annu Rev Plant Biol. 2007;58:183–98.

    CAS  PubMed  Google Scholar 

  85. Silverstone AL, Ciampaglio CN, Sun T. The Arabidopsis RGA gene encodes a transcriptional regulator repressing the gibberellin signal transduction pathway. Plant Cell. 1998;10(2):155–69.

    PubMed Central  CAS  PubMed  Google Scholar 

  86. Peng J, Harberd NP. Gibberellin deficiency and response mutations suppress the stem elongation phenotype of phytochrome-deficient mutants of Arabidopsis. Plant Physiol. 1997;113(4):1051–8.

    PubMed Central  CAS  PubMed  Google Scholar 

  87. Phillips AL, Huttly AK. Cloning of two gibberellin-regulated cDNAs from Arabidopsis thaliana by subtractive hybridization: expression of the tonoplast water channel, gamma-TIP, is increased by GA3. Plant Mol Biol. 1994;24(4):603–15.

    CAS  PubMed  Google Scholar 

  88. Kauschmann A, Jessop A, Koncz C, Szekeres M, Willmitzer L, Altmann T. Genetic evidence for an essential role of brassinosteroids in plant development. Plant J. 1996;9(5):701–13.

    CAS  Google Scholar 

  89. Cano-Delgado A. BRL1 and BRL3 are novel brassinosteroid receptors that function in vascular differentiation in Arabidopsis. Development. 2004;131(21):5341–51.

    CAS  PubMed  Google Scholar 

  90. Wang KL, Yoshida H, Lurin C, Ecker JR. Regulation of ethylene gas biosynthesis by the Arabidopsis ETO1 protein. Nature. 2004;428(6986):945–50.

    CAS  PubMed  Google Scholar 

  91. Liu Y, Zhang S. Phosphorylation of 1-aminocyclopropane-1-carboxylic acid synthase by MPK6, a stress-responsive mitogen-activated protein kinase, induces ethylene biosynthesis in Arabidopsis. Plant Cell. 2004;16(12):3386–99.

    PubMed Central  CAS  PubMed  Google Scholar 

  92. Han L, Li GJ, Yang KY, Mao G, Wang R, Liu Y, et al. Mitogen-activated protein kinase 3 and 6 regulate Botrytis cinerea-induced ethylene production in Arabidopsis. Plant J. 2010;64(1):114–27.

    CAS  PubMed  Google Scholar 

  93. Solano R, Stepanova A, Chao Q, Ecker JR. Nuclear events in ethylene signaling: a transcriptional cascade mediated by ETHYLENE-INSENSITIVE3 and ETHYLENE-RESPONSE-FACTOR1. Genes Dev. 1998;12(23):3703–14.

    PubMed Central  CAS  PubMed  Google Scholar 

  94. Fode B, Siemsen T, Thurow C, Weigel R, Gatz C. The Arabidopsis GRAS protein SCL14 interacts with class II TGA transcription factors and is essential for the activation of stress-inducible promoters. Plant Cell. 2008;20(11):3122–35.

    PubMed Central  CAS  PubMed  Google Scholar 

  95. Jiang Y, Liang G, Yang S, Yu D. Arabidopsis WRKY57 functions as a node of convergence for jasmonic acid- and auxin-mediated signaling in jasmonic acid-induced leaf senescence. Plant Cell. 2014;26(1):230–45.

    PubMed Central  CAS  PubMed  Google Scholar 

  96. Yan J, Zhang C, Gu M, Bai Z, Zhang W, Qi T, et al. The Arabidopsis CORONATINE INSENSITIVE1 protein is a jasmonate receptor. Plant Cell. 2009;21(8):2220–36.

    PubMed Central  CAS  PubMed  Google Scholar 

  97. Xie DX, Feys BF, James S, Nieto-Rostro M, Turner JG. COI1: an Arabidopsis gene required for jasmonate-regulated defense and fertility. Science. 1998;280(5366):1091–4.

    CAS  PubMed  Google Scholar 

  98. Thines B, Katsir L, Melotto M, Niu Y, Mandaokar A, Liu G, et al. JAZ repressor proteins are targets of the SCF(COI1) complex during jasmonate signalling. Nature. 2007;448(7154):661–5.

    CAS  PubMed  Google Scholar 

  99. Sheard LB, Tan X, Mao H, Withers J, Ben-Nissan G, Hinds TR, et al. Jasmonate perception by inositol-phosphate-potentiated COI1-JAZ co-receptor. Nature. 2010;468(7322):400–5.

    PubMed Central  CAS  PubMed  Google Scholar 

  100. Staswick PE. JAZing up jasmonate signaling. Trends Plant Sci. 2008;13(2):66–71.

    CAS  PubMed  Google Scholar 

  101. Dean JV, Mohammed LA, Fitzpatrick T. The formation, vacuolar localization, and tonoplast transport of salicylic acid glucose conjugates in tobacco cell suspension cultures. Planta. 2005;221(2):287–96.

    CAS  PubMed  Google Scholar 

  102. Hung KT, Kao CH. Nitric oxide counteracts the senescence of rice leaves induced by abscisic acid. J Plant Physiol. 2003;160(8):871–9.

    CAS  PubMed  Google Scholar 

  103. Lichtenthaler HK. Chlorophylls and carotenoids: pigments of photosynthetic biomembranes. Method Enzymol. 1987;148:350–82.

    CAS  Google Scholar 

  104. Saher S, Piqueras A, Hellin E, Olmos E. Hyperhydricity in micropropagated carnation shoots: the role of oxidative stress. Physiol Plant. 2004;120(1):152–61.

    CAS  PubMed  Google Scholar 

  105. Luan JB, Li JM, Varela N, Wang YL, Li FF, Bao YY, et al. Global analysis of the transcriptional response of whitefly to tomato yellow leaf curl china virus reveals the relationship of coevolved adaptations. J Virol. 2011;85(7):3330–40.

    PubMed Central  CAS  PubMed  Google Scholar 

  106. Morrissy AS, Morin RD, Delaney A, Zeng T, McDonald H, Jones S, et al. Next-generation tag sequencing for cancer gene expression profiling. Genome Res. 2009;19(10):1825–35.

    PubMed Central  CAS  PubMed  Google Scholar 

  107. Hoen PAC T, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RHAM, de Menezes RX, et al. Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res. 2008;36(21):e141.

    Google Scholar 

  108. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2− CT method. Methods. 2001;25(4):402–8.

    CAS  PubMed  Google Scholar 

  109. Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997;7(10):986–95.

    CAS  PubMed  Google Scholar 

  110. Saeed AI, Bhagabati NK, Braisted JC, Liang W, Sharov V, Howe EA, et al. TM4 microarray software suite. Methods Enzymol. 2006;411:134–93.

    CAS  PubMed  Google Scholar 

  111. Sturn A, Quackenbush J, Trajanoski Z. Genesis: cluster analysis of microarray data. Bioinformatics (Oxford, England). 2002;18(1):207.

    CAS  Google Scholar 

  112. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, et al. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39(suppl):W316–22.

    PubMed Central  CAS  PubMed  Google Scholar 

  113. Zhang H, Jin J, Tang L, Zhao Y, Gu X, Gao G, et al. PlantTFDB 2.0: update and improvement of the comprehensive plant transcription factor database. Nucleic Acids Res. 2011;39(Database issue):D1114–7.

    PubMed Central  CAS  PubMed  Google Scholar 

Download references


We thank Beijing Genomics Institute (BGI) at Shenzhen for assistance with sequencing and original data processing; Deyong Lai for help with bioinformatics analysis and for critically reading earlier versions of the manuscript. This work was supported by the National Basic Research Program of China (No. 2010CB126006) and the earmarked fund for the China Agriculture Research System (CARS-18).

Author information

Authors and Affiliations


Corresponding author

Correspondence to Shuxun Yu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

ML, CP, and SY designed and oversaw the research. ML and CP performed the research. ML, CP and HW performed the bioinformatics analysis. SF, MS and SY provided helpful suggestion during the data analysis. ML and CP analyzed the results and wrote the manuscript. All authors read and approved the final manuscript.

Additional files

Additional file 1:

List of gene differentially expressed during cotton development and senescence.

Additional file 2:

List of leaf senescence orthologs between cotton and other species.

Additional file 3:

Significantly enriched GO terms in each cluster.

Additional file 4:

All enriched pathways in cotton identified by KOBAS during leaf development and senescence.

Additional file 5:

List of transcription factors differentially expressed during leaf senescence.

Additional file 6:

List of eight plant hormone pathway-related genes.

Additional file 7:

Tag analysis statistics by using G. raimondii genome sequence.

Additional file 8:

List of some senescence-associated genes (SAGs) whose differential expression started at 25 d during leaf senescence.

Additional file 9:

List of primers used for qRT-PCR.

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit

The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lin, M., Pang, C., Fan, S. et al. Global analysis of the Gossypium hirsutum L. Transcriptome during leaf senescence by RNA-Seq. BMC Plant Biol 15, 43 (2015).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: