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

Background Cotton (Gossypium hirsutum L.) is one of the world’s most economically-important crops. However, its entire genome has not been sequenced, and limited resources are available in GenBank for understanding the molecular mechanisms underlying leaf development and senescence. Methodology/Principal Findings In this study, 9,874 high-quality ESTs were generated from a normalized, full-length cDNA library derived from pooled RNA isolated from throughout leaf development during the plant blooming stage. After clustering and assembly of these ESTs, 5,191 unique sequences, representative 1,652 contigs and 3,539 singletons, were obtained. The average unique sequence length was 682 bp. Annotation of these unique sequences revealed that 84.4% showed significant homology to sequences in the NCBI non-redundant protein database, and 57.3% had significant hits to known proteins in the Swiss-Prot database. Comparative analysis indicated that our library added 2,400 ESTs and 991 unique sequences to those known for cotton. The unigenes were functionally characterized by gene ontology annotation. We identified 1,339 and 200 unigenes as potential leaf senescence-related genes and transcription factors, respectively. Moreover, nine genes related to leaf senescence and eleven MYB transcription factors were randomly selected for quantitative real-time PCR (qRT-PCR), which revealed that these genes were regulated differentially during senescence. The qRT-PCR for three GhYLSs revealed that these genes express express preferentially in senescent leaves. Conclusions/Significance These EST resources will provide valuable sequence information for gene expression profiling analyses and functional genomics studies to elucidate their roles, as well as for studying the mechanisms of leaf development and senescence in cotton and discovering candidate genes related to important agronomic traits of cotton. These data will also facilitate future whole-genome sequence assembly and annotation in G. hirsutum and comparative genomics among Gossypium species.


Introduction
Cotton (Gossypium spp.) is the world's most important agronomic fiber, as well as a significant oilseed crop. The seed is an important source of feed, foodstuff, and oil. The crop is widely cultivated in more than 80 countries, with China, India, the United States of America, and Pakistan the top four cotton producers (http://www. cotton.org/econ/cropinfo/cropdata/rankings.cfm). China is the largest producer and consumer of raw cotton. Gossypium hirsutum L., or upland cotton, is a primary cultivated species and has an allotetraploid genome (AD; 2n = 4x = 52). Gossypium hirsutum produces over 90% of the world's fibers because of its higher yield and wider environmental adaptability [1,2].
The advent of new molecular genetic technologies and the dramatic increase in plant gene sequence data have provided opportunities to understand the molecular basis of traits important for plant breeding, such as improved yield and plant quality. The entire genomic sequence is not available for G. hirsutum, but a large number of genomic resources have been developed for this species. These include bacterial artificial chromosomes (BACs) [3], polymorphic markers [4], and genome-wide cDNA-based or unigene expressed sequence tag (EST)-based microarrays [5]. A rapid and cost-efficient way to acquire transcriptome data for an organism with a large, complex, and unknown genome is EST sequencing; analysis of ESTs can also complement whole-genome sequencing [6]. ESTs are short, single-pass sequence reads from mRNA (cDNA). Large scale EST data represent a snapshot of genes expressed in a given tissue and/or at a given developmental stage. They are tags of expression for a given cDNA library [7]. Large EST datasets can be used to discover novel genes and carry out functional genetic studies [8,9]. So far, a large number of G. hirsutum ESTs have been produced from cDNA libraries constructed from fibers, ovules, bolls, roots, and stems. The overwhelming majority of these EST resources have focused on fiber-development organs and have been used to explore the key genes involved in fiber development and its mechanism [10,11]. However, large-scale EST data related to leaf development are lacking.
Leaves are specialized photosynthetic organs, and plants harvest energy and nutrients in their production. Leaf development encompasses many distinct stages, from leaf primordium formation to expansion, maturation, and abscission. The onset and progression of leaf senescence, the last phase, is accompanied by changes in expression of a large number of senescence-associated genes (SAGs). Some genes must be newly activated in leaves for the onset of senescence [12,13]. Premature senescence, when the plant drops its leaves too early, has been occurring at an increasing frequency since the introduction of modern, high-yielding cotton cultivars like Bacillus thuringiensis (Bt) cotton. Premature leaf senescence results in reduced lint yield and poor fiber properties in cotton [14]. Understanding the molecular mechanisms of leaf senescence could greatly enhance yield and quality by guiding appropriate management to avoid premature leaf loss.
In recent decades, many advances in the understanding of leaf senescence at the molecular level have been achieved in several species, such as Arabidopsis thaliana and rice, by different experimental methods. Nine yellow-leaf-specific genes (YLS) were isolated, and RNA gel blot analysis revealed that most of them were senescence-up-regulated; the expression characteristics of YLS genes will be useful as potential molecular markers [15]. Transcript abundance in leaves of Populus tremula was studied by microarrays obtained from seven cDNA libraries, and 677 significantly up-regulated genes were identified during leaf senescence. The evidence for increased transcriptional activity before the appearance of visible signs of senescence was also found [16]. In Medicago truncatula leaves, 545 differentially-expressed genes, including 346 senescence-enhanced and 199 repressed genes, were identified by cDNA amplified fragment length polymorphism (AFLP) techniques; comparison with Arabidopsis datasets revealed common physiological events but differences in nitrogen metabolism and transcriptional regulation [17]. In rice, 533 differentially expressed genes were isolated by suppression subtractive hybridization (SSH) from early-senescent flag leaves, 183 had gene ontology (GO) annotations indicating involvement in macromolecule metabolism, protein biosynthesis regulation, energy metabolism, detoxification, pathogenicity and stress, and cytoskeleton organization [18]. A total of 140 annotated upregulated genes in wheat flag leaves were analyzed using an inhouse fabricated cDNA microarray. The results supported a protective role of mitochondria towards oxidative cell damage via the up-regulation of an alternative oxidase and possibly also succinate dehydrogenase [19]. During natural leaf senescence in Arabidopsis, 827 SAGs were identified. Comparison of these genes with artificially-induced senescence suggested that alternative pathways for essential metabolic processes such as nitrogen mobilization were used in different senescent systems [20]. Recently, a high-resolution time-course profile of gene expression during leaf senescence was obtained by microarray analysis. The dynamic changes in transcript levels were identified globally as senescence progresses, and the involvement of metabolic processes, signaling pathways, and specific transcription factors (TFs) were explicitly clarified [21]. Among the SAGs, many TFs, receptors, signaling components for hormones and stress responses, and regulators of metabolism were involved in regulating leaf senescence, indicating that senescence is governed by complex transcriptional regulatory networks.   In this study, a normalized and full-length cDNA library from different developmental stages of G. hirsutum leaves was constructed. Random sequencing of clones from the cDNA library generated a total of 9,874 high-quality ESTs, which were assembled into 5,191 unique sequences, consisting of 1,652 contigs and 3,539 singletons. Several SAGs and TFs were identified. This work will benefit the study of leaf senescence mechanisms of G. hirsutum, form a foundation for cloning the full-length sequences of these genes for genetic engineering, and also provide important resources for comparative genomic studies among closely-related species.

Characterization of cDNA Library and EST Sequencing and Assembly
A normalized full-length cDNA library was constructed using leaves during the plant flowering stage. To evaluate the fullness ratios of the cDNA inserts of the library, 50 clones were randomly selected and fully sequenced; 44 (88%) contained putative fulllength sequences. To assess the normalization efficiency, the relative concentration of 18S ribosomal RNA (18S) and actin in both the non-normalized and normalized cDNA populations were estimated by quantitative real-time PCR (qRT-PCR). The differences in cycle number (DCt) increased by 7.18 and 7.65 after library normalization, respectively. The results showed that the copies of these two genes decreased 145 and 200 fold, respectively, and suggested that the normalization quality of this library was good.
Approximately 11,623 clones were successfully single-pass sequenced from their 39 ends. The insert sizes ranged from 900-3,000 bp, with an estimated average size of 1,200 bp. After removal of vector, poly(A) tails, contaminating microbial sequences, and those shorter than100 bp, 9,874 ESTs were considered high confidence (Q20) and were deposited in the GenBank dbEST database (JZ110066-JZ119939). Clustering and assembly of the ESTs were carried out under stringent conditions to obtain 5,191 putative unigenes, including 1,652 (31.8%) contigs that consisted of two or more ESTs and 3,539 (68.2%) singletons ( Table 1). The EST redundancy of this library was 47.4%, and the unigenes had an average length of 682 bp. The distribution of high-quality EST sequences in the clusters is shown in Figure 1. Of the 1,652 contigs, 885 (53.6%) contained two ESTs, 363 (22.0%) contained three ESTs, 156 (9.4%) contained four ESTs, 86 (5.2%) contained five ESTs, 47 (2.8%) contained six ESTs, and relatively few (7.0%) contained more than six ESTs ( Figure 2). The unigene mean size was only 1.9 sequences, and each contig averaged 3.8 sequences.
These results also suggested that the redundancy rate of this normalized library was relatively low.
The most abundant ESTs are shown in Table 2 (each contig contained $10 EST copies). Some of these genes have important roles in leaf senescence. For example, lipid transfer protein precursors (Contig1797, Contig1794, Contig1792, and Con-tig1778; 97 ESTs) were involved in nutrient recycling for lipid transfer. Cysteine protease (Contig1795), polyubiquitin (Con-tig1784), putative serine carboxypeptidase precursor (Contig1749), and aspartyl protease (Contig1748) play roles in protein degradation. Peroxidase (Contig1766) is an antioxidant important for redox regulation, and metallothionein-like protein (Contig1756) is a low-molecular-weight Cys-rich protein that functions in heavy metal detoxification to remobilize valuable metal ions. Among these unigenes, 4,876 (93.9%) had open reading frames (ORFs) that were longer than 100 bp. The average ORF was 342 bp. The mean G/C content of unigenes was approximately 42%, which was approximately equivalent to that of Arabidopsis (43.2%) and much lower than that of rice (55.2%) [22,23].

Unigene Functional Annotation and Functional Categorization
To annotate the unigenes, all unigenes were used in a blastx search against the NCBI non-redundant (nr) protein database with a cut-off E-value of 10 25 . The nr database is commonly used as the principal target database to search for homologous proteins. Using this approach, most unique sequences (84.4%) had matches in the nr database. However, 808 sequences had no hits. Of the best matches, 1,094 (25%) were to Ricinus, 1011 (23.1%) to Vitis, 975 (22.2%) to Populus, 158 (3.6%) to Arabidopsis, 142 (3.2%) to Glycine, 49 (1.1%) to Medicago, and 27 (0.6%) to Oryza, whereas only 490 (11.2%) of the best matches were to cotton ( Table 3). Comparison of our unigene data set with the NCBI nucleotide database using blastn demonstrated that 4,138 unigenes (79.7%) had significant matches. All unigenes were also blastx searched against the Swiss-Prot database, in which 2,973 (57.3%) unigenes matched. The best hits were mainly to Arabidopsis (1,514 hits, 50.9%) and rice (151 hits, 5.1%).

Comparison with Previous Cotton ESTs
To evaluate potential novel sequences that did not match to sequences from other cotton species in the databases, the unigenes were used as queries in a blastn search against the Dana-Farber Cancer Institute (DFCI) Cotton Gene Index database. Approximately 24.3% of the ESTs and 19.1% unique sequences generated in this study were not highly homologous to known cotton ESTs or unique sequences. Thus, our library provides a valuable new transcript resource, with 2,400 new ESTs and 991 new unique sequences for cotton.

Identification and Analysis of Leaf Senescence-related Protein Families
To identify leaf SAGs, all the unigenes were assessed by blastx against the amino acid sequences of A. thaliana genes from the leaf senescence database (LSD). Of the 5,191 unigenes, 1,339 (25.8%) had homologs matched with 455 (44.6%) SAGs of A. thaliana in the LSD, and could be classified into 29 leaf senescence-related categories ( Table 5). The most abundant leaf senescence-related category was protein degradation/modification, with a total of 199 unigenes. Other highly-abundant leaf SAGs included Nutrient recycling, Lipid/Carbohydrate metabolism, Signal transduction, Transcriptional regulation, Redox regulation, Stress and detoxification, and Hormone response pathway. These functions are all closely involved with leaf senescence.
To study the expression of genes associated with leaf senescence, first, representative leaves were classified as young leaves (Y), mature leaves (M), early-senescent leaves (S1) and late-senescent leaves (S2) by their chlorophyll contents, as shown in Figure 4a. The chlorophyll content in S1 and S2 was 65% and 45% of that in M, respectively. Then, nine putative leaf senescence-related ESTs were randomly selected for qRT-PCR using RNA isolated from leaves of those four stages. Most of these ESTs were up-regulated during senescence, especially Contig773, whose expression level increased significantly in the late senescent leaves (Figure 4b). Only two ESTs, JZ110587 and JZ116048, were down-regulated. Characterization of these potential regulatory genes provided clues to the regulatory mechanism of leaf senescence.
To identify the potential roles of these TFs during leaf senescence, the most abundant MYB family in this normalized library was selected and its expression pattern was analyzed. As shown in Figure 5, several putative cotton MYB orthologs matched AtMYBL (AT1G49010), ZmMYB153 (GRMZM2G050550), AtMYB (AT4G01280), AT3G24860.1, AtMYBR1 (AT5G67300), AT3G52250.1, and AT4G37180.2, which play roles in leaf senescence [26][27][28]. Using qRT-PCR, we confirmed the transcript abundance of selected ESTs encoding putative MYB TFs in leaves at different developmental stages ( Figure 6). JZ118495, JZ116679 and JZ112479 were putative cotton orthologs of AtMYBL (AT1G49010), AT3G24860.1 and AT4G37180.2, respectively, from A. thaliana ( Figure 5). The expression of JZ118495 increased in M stage and reached a maximum in S1 stage, while that of JZ116679 increased gradually during leaf senescence and peaked in S2 stage, and JZ112479 was expressed at high levels in the S1 stage but at reduced levels in the S2 stage ( Figure 6). Six of 11 ESTs were highly expressed in senescent leaves; most increased in the expression level in the S1 stage, including JZ110276, JZ112420, JZ112479, and JZ118495. Other transcripts, such as Con-tig1167, JZ111255, Contig 1171, Contig708 and JZ112513 were down-regulated during leaf senescence. The results indicated that these MYB TFs may be involved in controlling leaf senescence in cotton.

and Expression Analyses
To confirm that our full-length library was an efficient method for rapid functional gene discovery in upland cotton, three A. thaliana homologs of yellow-leaf-specific genes (YLS) were cloned and analyzed. The Arabidopsis YLS proteins were used as queries to search our EST database with tBLASTn. Three unique fulllength sequences were found in upland cotton and named GhYLS5 (JX163920), GhYLS8 (JX163921), and GhYLS9 (JX163922). In A. thaliana, the YLS5 gene encoded a proteaseI (pfpI)-like protein of 398 amino acid residues that was expressed weakly in young leaves and strongly in senescent leaves. This gene can be induced by artificial senescence processes such as darkness, ethylene, and ABA treatment [15]. The GhYLS5 gene had an ORF of 1,188 bp and encoded a protein of 395 amino acid residues. Multiple sequence alignment showed that GhYLS5 proteins were homologous to the glutamine amidotransferase (GAT) of A. thaliana and Theobroma cacao and with YLS5 of A. thaliana, Arabidopsis lyrata and Zea mays with identities of 51-84% (Figure 7a, b). Arabidopsis YLS8 contained an ORF encoding a Dim1 homolog of 142 amino acid residues that had high expression in senescent and virus-infected leaves [15,29]. The GhYLS8 gene had an ORF of 429 bp, encoding a protein of 142 amino acid residues. The protein of GhYLS8 was highly conserved, with very high sequence homology to YLS8 from A. thaliana, Hevea brasiliensis, Matthiola longipetala, Iberis amara, and Lepidium sativum and to thioredoxin-like protein 4A (TRX4A) from A. thaliana, Cucumis sativus, Vitis vinifera and Medicago truncatula (Figure 8a, b). The YLS9 gene (also called NHL10) of Arabidopsis contained an ORF encoding a polypeptide of 227 amino acid residues, whose sequence was similar to tobacco hairpin-induced gene (HIN1) and Arabidopsis non-race specific disease resistance gene (NDR1). Expression of this gene is induced by Cucumber mosaic virus, spermine, and senescence [15,29,30].GhYLS9 gene had an ORF of 669 bp, encoding a protein of 222 amino acid residues. GhYLS9 proteins were homologous to syntaxin (SYP) from Ricinus communis, Cucumis sativus and Glycine max, HIN1 from Casuarina glauca and Nicotiana tabacum, and YLS9 from A. thaliana, with identities of 51-62% (Figure 9a, b). The expression of three GhYLS transcripts were also analyzed using qRT-PCR at different leaf developmental stages (Figures 7c, 8c, 9c). The three genes were all up-regulated in senescent leaves. In particular, expression of GhYLS9 was nearly 400-fold higher than in young leaves. These results suggested that leaf senescence related-genes could be identified from our library using -homologous sequence searches.

Discussion
Gossypium hirsutum is one of the most economically-important species in its genus. Unfortunately, to date, its genome has not been completely sequenced. Recent efforts have demonstrated that EST sequencing is an efficient and relatively low-cost approach for large-scale gene discovery, annotation, and comparative genomics research [31]. In G. hirsutum, although many ESTs are available, the total number is less than that of some field crops and model plants, and most ESTs in GenBank are from fibers or fiber-bearing ovules [11,[32][33][34] and provide little or no information regarding leaf development. Therefore, G. hirsutum leaf ESTs must be sequenced to examine the functional genomics of cotton leaf development. In this study, we produced 9,874 high-quality ESTs that assembled into 5,191 unigenes from a normalized leaf cDNA library. The leaf samples spanned all development stages, including unexpanded young leaves, fullyexpanded mature leaves, and senescent leaves, at the plant blooming stage. This is the first such database and largest number of unique sequences from G. hirsutum leaf tissues to include all developmental periods. This EST resource provides a foundation for molecular control of G. hirsutum leaf growth and development and for future whole-genome sequencing and analysis of the functional genome and gene expression patterns.
Normalized cDNA libraries overcome problems caused by differential expression of genes and are an efficient and costeffective tool for obtaining large-scale unique EST sequences and for gene identification [35]. Our cDNA library was normalized by saturation hybridization with genomic DNA, assuming relatively uniform copy numbers of most of genes in the genome. EST assembly revealed a novelty rate of 52.6%, a redundancy rate of 47.4%, and 68.2% of unigenes that contained only one EST. Thus, there remains considerable potential to discover additional novel sequences by sequencing randomly-selected cDNAs from this library. Alpha-tubulin 10 (TUA10) and ubiquitin (UBI1), the most redundant transcripts in cotton leaves, were represented by only 19 and 17 clones in our ESTs, respectively. Furthermore, the copies of two highly abundant genes actin and 18S, decreased 145 and 200 fold after cDNA library was normalized, respectively. These results reflect the quality of the normalized library and also showed that this approach was an efficient tool for gene identification because it reduced variation among abundant clones and increased the probability of sequencing rare transcripts. The majority of annotated sequences with BLAST hits were transcripts from the rosid clade, to which cotton (Malvales; eurosids II clade) also belongs. Ricinus (25% of the best matches) and Populus (22.2%) belong to the eurosids I clade, while Vitis (23.1%) is a basal rosid [36]. Although A. thaliana and O. sativa are well-studied model systems with completely-sequenced genomes, these organisms were best matches to only 3.6% and 0.6% of our unique sequences, respectively. Yu [37] investigated the conservation of colinearity between cotton BAC sequences and other model plant genomes; on a phylogeny of single-copy orthologous genes from cotton, Arabidopsis, poplar, grape, rice, and maize, poplar was the closest relative to cotton. Arabidopsis thaliana, P. trichocarpa, and G. hirsutum are dicots, while O. sativa is a monocot, which may account for the differences in similarity among their sequences. Only 11.2% of the hits were to cotton sequences already available in GenBank, highlighting the lack of sequence information for this genus and the value of our EST sequences. Clearly, genome sequencing of G. hirsutum represents a vital and urgent need. Furthermore, we discovered 2,400 new cotton ESTs and 991 unique cotton sequences when comparing our data to the DFCI Cotton Gene Index database. Our data will contribute to the enrichment of cotton genetic and physical maps.
In previous studies, much attention was focused on leaf senescence, especially in Arabidopsis and rice [18,[38][39][40][41][42]. Leaf senescence constitutes the last stage of leaf development and strongly affects cotton yield. Currently, however, the dynamic regulatory mechanisms of leaf senescence in cotton remain unclear. A large number of SAGs have been identified in various plants through microarray analyses [21,43]. Some of them have been found to be TFs belonging to several different families, especially NAC, WRKY, C2H2-type zinc finger, AP2/EREBP, and MYB protein families [44,45]. Characterization of these potential regulatory genes led to discovery of a few important senescence regulatory genes and provided some insight into the regulatory mechanism of leaf senescence. Using data from the PlantTFDB 2.0 database, we found 200 unigenes from our library that had high similarity with 163 TFs from 41 families. The most well-represented TF family in our library was the MYB group, followed by the bHLH, bZIP, C3H, NAC, ERF, ARF, C2H2, and WRKY families. Analysis of the expression patterns of several putative MYB transcript families showed that the expression level of some transcripts changed significantly during leaf senescence. This result indicated that some MYB TFs may play roles during leaf senescence. These results also were in accordance with those of previous studies. In addition to these TF families, several others known to be involved in plant development were also present in our data.
Leaf senescence is an integrated response of leaf cells to age and other internal and environmental signals. It is an exceptionally complex and dynamic genetic process [46]. Arabidopsis thaliana is a favorite model for the molecular genetic study of leaf senescence [47][48][49]. The LSD is also a platform to study leaf senescence [50]. Of the unigenes in our library, 1,339 could be classified into 29 SAG categories by a BLAST search against A. thaliana senescence-related proteins (1,021), such as nutrient recycling, Lipid/Carbohydrate metabolism, and hormone response pathway. During leaf senescence, nutrients in the leaf are reallocated to younger leaves, growing seeds, or other growing organs in a process of nutrient salvage, e.g., hydrolysis of macromolecules and subsequent remobilization, which requires complex array of metabolic pathways [51]. Many the genes involved in lipid metabolism function in leaf senescence. Lipid-degrading enzymes, such as lytic acyl hydrolase, phosphatidic acid phosphatase, phospholipase D, and lipoxygenase appear to be involved in hydrolysis and metabolism of the membrane lipid in senescing leaves [51,52]. Changed expression of the Arabidopsis acyl hydrolase gene in transgenic plants led to altered leaf senescence phenotypes [53]. The hormonal pathways appear to affect all stages of leaf senescence. In this work, numerous genes belonging to hormone response pathways were also identified. These results indicated that many previously-known leaf SAGs and pathways were included in this library. Three GhYLS genes were successfully cloned and analyzed. Their expression profiles revealed that their transcripts accumulated in leaves during senescence. Thus, these genes could potentially serve as molecular markers for distinguishing the complex regulatory networks of leaf senescence processes. This library provides a robust sequence resource and will be a useful tool for cloning the full-length sequences of functional genes for further leaf senescence-related analysis in G. hirsutum.

Plant Material
Upland cotton CCRI 36 (a short-season cultivar) was grown on the experimental farm of the Cotton Research Institute of Chinese Academy of Agricultural Sciences, Anyang, Henan Province. At the blooming stage, unexpanded leaves of the same size near the tops of stems were selected and marked. The day when leaves were fully expanded was considered the first day. Leaves were collected every 5 d for 70 d. Samples from each time point were pooled from at least 10 plants, frozen immediately in liquid nitrogen, and stored at -80uC.

RNA Isolation and cDNA Library Construction
Total RNA was isolated by an improved CTAB method [54], and equal amounts of total RNA sampled at different time points were mixed to construct a full-length normalized cDNA library. Purification of mRNA from total RNA was carried out using the FastTrackH 2.0 Kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. cDNAs were synthesized using the Superscript Full-length Library Construction Kit II (Invitrogen) according to the manufacturers' protocols, cloned into a Gateway pDONR222 vector (Invitrogen) by the BP cloning process, and transformed into Escherichia coli strain DH10B competent cells (Invitrogen) through electroporation using an E. coli Pulser (BTX Harvard Apparatus, Holliston, MA, USA). After the full-length library was constructed, plasmid DNA was extracted with the PureLink TM HQ Mini Plasmid DNA Purification Kit (Invitrogen). Normalization was performed by saturation hybridization between genomic DNA and mixed plasmid DNA from the cDNA library [55]. Then, clones were randomly selected and fully sequenced to test fullness ratios of the cDNA inserts of the library. Putative full-length cDNA sequences were identified by comparison with all available ORFcomplete mRNA sequences from the NCBI nr protein database [56]. Finally, qRT-PCR was used to estimate the relative concentration of a highly abundant clone in both the nonnormalized and the normalized cDNA populations.  EST Sequencing, Editing, and Assembly Clones were randomly picked and transferred into 384-well plates. Selected clones were sequenced from the 39 end on an ABI 3730 automatic DNA sequencer (Applied Biosystems, Foster City, CA, USA) using the M13 universal primer (M13R:CAGGAAA-CAGCTATGACC) and the BigDye Terminator Cycle Sequencing Kit (ABI) at the Invitrogen Sequencing Center. All sequences were clustered using the Phred/Phrap/Consed software package [57,58]. The 39 DNA EST sequence chromatogram files were base-called and quality trimmed (low-quality bases with,Q20 and ,99% accuracy were removed) using Phred. Crossmatch (http:// www.phrap.org/) and Repeat-Masker (http://www.repeatmasker. org/) were used to remove vector sequences and to identify and mask repeat sequences. Contaminating microbial sequences were eliminated using VecScreen (http://www.ncbi.nlm.nih.gov/ VecScreen/VecSc-reen.html), and poly(A) tails were deleted. Sequences that passed the quality control screening for highconfidence base calls (Q20) and with lengths longer than 100 bp were defined as high quality EST and deposited into the dbESTs division of GenBank. The processed EST sequence files were combined and assembled into contigs and singlets (unisequences) using Phrap with a high stringency level (95% sequence identity with 20 bp overlap).
To validate potential novel ESTs and unique sequences that did not match any sequences in related cotton species in the existing databases, all the high-quality ESTs and assembled unigenes were compared against ESTs and unigenes already available in the DFCI Cotton Gene Index (http://compbio.dfci.harvard.edu/ cgi-bin/tgi/gimain.pl?gudb = cotton) database, which contains 351,954 cotton ESTs and 2,315 ETs fully assembled into 117,992 unique sequences. With such stringent criteria, an EST was considered as new if it had at least 10% of its sequence with less than 95% of identity to any other EST or unigene in the public EST database.

Prediction of ORFs, Unigene Functional Annotation, and Functional Categorization
All unique sequences were searched for putative ORFs with the Getorf program of EMBOSS-4.1.0 [59], and the longest sequences were used for functional analysis. Unigenes were compared with a variety of databases, including the NCBI nonredundant nucleotide and non-redundant protein databases, and Swiss-Prot, using either blastn (E-value#10 25 ) or blastx (E-value#10 25 ) [60]. To identify putative leaf SAGs and TFs, blastx (E-value#10 25 ) searches against amino acid sequences of A. thaliana genes from a leaf senescence database (LSD) [55,61] and a comprehensive plant TF database (PlantTFDB) [25] were used. Batch searches of the unigenes were performed using the local BLAST tools available at ftp://ftp.ncbi.nlm.nih.gov/blast/ executables/blast+/LATEST/. To assign GO terms, functional annotation was performed using Blast2GO software based on sequence similarity [62][63][64]. Furthermore, to improve annotations, results from an InterProScan search [65] (http://www.ebi. ac.uk/interpro/index.html) were merged with GO annotations and searched in the BlastProDom, FPrint-Scan, HMMPIR, HMMPfam, HMMSmart, HMMTigr, ProfileScan, ScanRegExp, and SuperFamily databases.

Leaf Senescence Related Homolog Identification and Expression Pattern Analysis
To examine gene expressions during leaf development, the leaves used for qRT-PCR were harvested from approximately 10 individual plants for each stage. Total chlorophyll of the samples was measured as described by Lichtenthaler (1987) [66].Homologs of leaf senescence-related protein sequences were identified and randomly selected according to the LSD function annotation. Total RNA was extracted by an improved CTAB method as described above. cDNA was reverse transcribed from RNA by PrimeScriptH RT Reagent Kit with gDNA Eraser (Takara, Otsu, Japan) with an Oligo dT Primer and random six-mers as the RT primer according to the manufacturer's protocol. The specific primer pairs for nine selected genes and the internal control gene actin are listed in Table S1. qRT-PCR was performed with the SYBR Green PCR Master Mix (Takara) as recommended by the manufacturer in an ABI 7500 Real-time PCR System (Applied Biosystems) with three replicates. To analyze changes in gene expression, values from triplicate real-time PCRs were normalized to the expression level of actin and to the Y sample by the 2 -DDCt method [67]. Arabidopsis YLS genes were used as queries to tBLASTn search against the cDNA library. The identified clones were sequenced in both directions with the internal primers. The amino-acid multiple-sequence alignment was analyzed using GeneDoc. Phylogenetic analysis was performed using the neighbor-joining method in MEGA 4 [68]. Expression patterns were detected by qRT-PCR as described above.

Supporting Information
Table S1 Primers used in gene-specific qRT-PCR of leaf senescence related genes.