Uncovering full-length transcript isoforms of sugarcane cultivar Khon Kaen 3 using single-molecule long-read sequencing

Sugarcane sampling, RNA extraction, and cDNA library construction

A leaf sample was collected from a 3-month old plant of the commercial sugarcane cultivar Khon Kaen 3 (Tippayawat, Ponragdee & Sansayawichai, 2012). The plant was cultivated at the National Center for Genetic Engineering and Biotechnology, National Science and Technology Development Agency, Thailand. No specific permits were required and no protected species were involved in this study. Young leaf (3 g) was collected for total RNA extraction using a QIAGEN RNeasy^® Mini Kit. Total RNA was treated with DNaseI (DNA-free™ kit; Ambion, Foster City, CA, USA) to remove genomic DNA contamination. From 100 µg of total RNA, mRNA (600 ng) was isolated using an Absolutely mRNA Purification kit (Agilent Technologies, Santa Clara, CA, USA).

PacBio Iso-Seq libraries were prepared according to the Isoform Sequencing (Iso-Seq) protocol using the Clontech SMARTer PCR cDNA Synthesis Kit and the BluePippin Size Selection System protocol. In brief, a 35 ng sample of mRNA was reverse-transcribed to synthesize first-strand cDNA. The first-strand cDNA was amplified for 12 PCR cycles using KAPA HiFi DNA polymerase (Kapa Biosystems, Wilmington, MA, USA). cDNA was size-selected into 4 bins of size ranges 1–2 kb, 2–3 kb, 3–6 kb and 5–10 kb using a BluePippin Size-Selection system (Sage Science, Beverly, MA, USA). cDNA from each bin was amplified for 15 PCR cycles, using the following cycles: 95 °C for 2 min follows by 15 cycles of 98 °C for 20 s, 65 °C for 15 s, 72 °C for 4 min and a final extension of 72 °C for 5 min. The amplified DNA was subjected to damage and end-repair, and then ligated to a SMRTbell adapter for Iso-Seq SMRTbell library construction. Before sequencing, a sequencing primer was annealed to SMRTbell templates which were then bound to DNA polymerase and prepared for Magbead loading onto SMRT cells. Sequencing was performed on a PacBio RS II system using P6C4 polymerase and chemistry and a 240-min movie time. Libraries composed of transcripts in the size ranges 1–2 kb, 2–3 kb, and 3–6 kb were each sequenced on two cells and transcripts of size 5–10 kb were sequenced on one cell.

100 µg of purified mRNA from the same leaf sample was used to construct an RNA-Seq library using the Ion Total RNA-Seq Kit v2 (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instruction. The mRNA was fragmented using RNaseIII (Thermo Fisher Scientific, Waltham, MA, USA) to obtain fragments with a mean length of around 200 nt. Subsequently, adapters were ligated to fragmented mRNA and first strand cDNA synthesis was performed using SuperScript™ III (Thermo Fisher Scientific, Waltham, MA, USA). The DNA library was sequenced on the Ion Proton using Ion PI™ Chip Kit v3 chips (Thermo Fisher Scientific, Waltham, MA, USA).

Sequencing data processing, isoform clustering and error correction

The raw reads generated from SMRT sequencing were processed separately for each size-selected library using the “RS_IsoSeq” protocol (https://github.com/PacificBiosciences/cDNA_primer/wiki) (Gordon et al., 2015) available in SMRT^® Analysis software v2.3.0 (http://www.pacb.com/support/software-downloads), a free software suite that can be downloaded from the Pacific Biosciences website. The analysis started by generating read of insert (ROI), the representative single sequence of highest quality for an insert generated by self-correction through subread comparison. The ROIs were then classified into full-length non-chimeric and non-full-length reads. The full-length reads with minimum length of 300 bp were clustered and used to construct consensus isoform sequences by the ICE (isoform-level clustering) algorithm. These consensus sequences were then polished using the Quiver algorithm to get high quality isoform sequences with at least 99% accuracy. High quality and low quality isoform sequences from all size fractions were combined into a final transcript sequence dataset for subsequent analyses. To improve the quality of polished PacBio transcripts, short-read error correction was performed on unique transcript isoforms using Ion Proton reads sequenced from the same RNA sample by using LoRDEC software version 0.6 with parameters of -k 21 -s 3 -t 5 -b 200 -e 0.4. Ion proton reads used for correction were preprocessed by sequence trimming of poor-quality base calls (N) and removing reads with ambiguous bases, reads of low overall quality (average Phred quality score < 17), and reads shorter than 50 bp using PRINSEQ version 0.20.4 (Schmieder & Edwards, 2011). BLASTN search against Rfam database version 12.2 (Nawrocki et al., 2015) was performed to remove rRNA sequences from the preprocessed Ion Proton reads before PacBio error correction process. Redundant sequences in the corrected transcripts were removed using CD-HIT (Fu et al., 2012; Li & Godzik, 2006) with 99% sequence identity cutoff. PacBio and Ion Proton raw read data were deposited in the NCBI Sequence Read Archive (SRA) database under accession number SRP103465.

Prediction of coding regions within transcripts

The identification of open reading frame (ORF) within transcript sequences was performed using Transdecoder software version 3.0.1 (Haas et al., 2013) (https://transdecoder.github.io). A Markov model was trained using the top longest ORFs initially identified across all transcripts. The software then identified candidate protein-coding regions on each transcript based on log-likelihood scores calculated from this model. Homology based on BLASTP (version ncbi-blast-2.6.0+) (Altschul et al., 1990) results of all translated ORFs against Phytozome version 12 (Goodstein et al., 2012) plant proteins and the presence of protein domain(s) from the Pfam database (version 31.0) (Finn et al., 2016) were incorporated as additional criteria for ORF prediction. The parameter settings for BLASTP and Pfam were as shown in the Transdecoder GitHub (https://transdecoder.github.io). A single best ORF of minimum length 300 bp was finally selected as the predicted ORF for each transcript. Protein domains in the assigned ORF of constructed transcripts were predicted using InterProScan 5.22–61.0 with default settings which search against protein domains/signatures from a number of different databases (Jones et al., 2014).

Functional annotation and similarity matching of PacBio isoforms

Sequence similarity search was performed on PacBio transcripts against various nucleotide and protein sequence databases using BLASTN and BLASTX (version blast-2.2.26), respectively (Altschul et al., 1990) with an E-value cutoff of 10⁻⁶. The match was assigned to the first hit with the highest bit score. Homology search among previously described sugarcane sequences was performed using the PlantGDB version 157a (http://www.plantgdb.org), SoGI release 3 (http://compbio.dfci.harvard.edu/tgi/) and NCBI dbEST release 130101 (Boguski, Lowe & Tolstoshev, 1993). Homology search among other plant sequences was conducted using the Phytozome version 12 database (Goodstein et al., 2012). Other homologous functions were annotated by searching the NCBI non-redundant protein database (nr) (downloaded in August 2016). Finally, sequence similarity search was performed on PacBio transcripts against long non-coding RNA (lncRNA) sequences downloaded from lncRNAdb v2.0 (Amaral et al., 2011), CANTATAdb v1.0 (Szczesniak, Rosikiewicz & Makalowska, 2016), NONCODE 2016 (Zhao et al., 2016) and Plant non-coding RNA Database (PNRD) (Yi et al., 2015).

To identify PacBio transcripts corresponding to putatively full-length transcripts spanning protein coding sequences (CDS), BLASTN search results of PacBio transcripts against plant CDS of Phytozome version 12 database (Goodstein et al., 2012) were used for calculation of length coverage of PacBio transcripts on the corresponding plant CDS with the highest bit score. Full-length transcripts were defined as PacBio transcripts with at least 90% coverage on matched plant CDS. The length ratios of PacBio transcripts to matched sequences from sugarcane databases (PlantGDB, NCBI dbEST, SoGI), PacBio transcripts from the SUGIT database (Hoang et al., 2017), and Phytozome plant transcript version 12 database were also calculated.

Gene identifiers (e.g., gene names, Gene Ontology (GO) terms, E.C. number, KOG ID and KEGG ID) were assigned to each PacBio transcript based on the hit with the highest bit score from BLASTN search against the Phytozome version 12 plant transcript sequences. The assigned GO terms were summarized in three categories, namely cellular component, molecular function, and biological process, using WEGO software (Ye et al., 2006). The assigned KOG ID is the eukaryotic orthologous group (KOG) which was constructed based on orthologous and paralogous proteins from different eukaryotic genomes to infer possible function of transcript, information of which is available in NCBI COG database (https://www.ncbi.nlm.nih.gov/COG) (Tatusov et al., 2003). To determine metabolic and biological pathways that are related to transcript sequences, KEGG Orthology (KO) associated with orthologous proteins was used in order to link to KEGG pathways (Kanehisa & Goto, 2000; Kanehisa et al., 2016). Only plant-related pathways of the KEGG database were applied.

PacBio isoforms from this study were matched with isoforms in the SUGIT database (Hoang et al., 2017) by BLASTN (version blast-2.2.26) (Altschul et al., 1990) analysis with an E-value cutoff of 10⁻²⁰.

Prediction of alternatively spliced transcripts

PacBio isoforms from this study and in the SUGIT database (Hoang et al., 2017) were aligned against the sorghum genome (Sbicolor_454_v3.0.1) from Phytozome version 12 using GMAP software (Wu & Watanabe, 2005) version 2015-06-23 with the parameter settings recommended for alignment of PacBio Iso-Seq transcripts to the genome (—no-chimeras—cross-species -expand -offsets 1 -B 5 -K 8000 -f samse -n 1 -t 8—min-trimmed-coverage=0.9 -min-identity=0.8) provided in the TAPIS (Transcriptome Analysis Pipeline from Isoform Sequencing) Bitbucket repository (https://bitbucket.org/comp_bio/tapis). The potential alternative splicing events of the PacBio constructed transcripts were predicted using TAPIS software version 1.2.1 (Abdel-Ghany et al., 2016) with default settings. The alternative splicing patterns of each gene was displayed using SpliceGrapher included in the TAPIS softtware.

Repeat sequence analysis

Repetitive sequences on PacBio transcripts were predicted using RepeatMasker version open-4.0.5 (Smit, Hubley & Green, 2013–2015) with default parameters. The transcripts were compared with classified sequences in “grasrep.ref.classified.fa”, which is a database that includes all grass family repeat sequences from RepBase Update 20140131 (Jurka et al., 2005) and RM database version 20140131.

Gene expression analysis

The preprocessed Ion Proton reads were mapped to PacBio transcripts using TMAP software (https://github.com/iontorrent/TMAP) with the ‘mapall’ parameter and output filtered for all best hits (-a 2). The output of TMAP in BAM formatted file was used for calculation of transcript expression levels (RPKM; Reads Per Kilobase of transcript per Million mapped reads) using tigar2 software (Nariai et al., 2014) with default parameters. The transcript sequences were classified into three groups based on expression levels: low (≥0.5 to 5 RPKM), medium (≥5 to 50 RPKM) and high (≥ 50 RPKM). Gene name and KOG ID were assigned to PacBio transcripts based on the top hit from BLASTN search against the Phytozome version 12 plant transcript sequence database.

Identification of sugarcane KK3 transcript isoforms from Iso-Seq data

PacBio transcriptome sequencing from a cDNA library of KK3 leaf sample produced a total of 4,159,174 reads of insert from all RNA size fractions. Out of these, 232,328 full-length, non-chimeric reads with a minimum length of 300 bp were obtained. After clustering of these full-length transcripts and polishing the consensus isoforms, 53,999 high-quality and 75,975 low-quality isoforms were generated. A total of 99,344,049 Ion Proton reads were filtered to 83,393,993 high-quality reads, which were then used for error correction of these Iso-Seq transcripts. After error correction of all 129,974 isoform sequences and removal of redundant sequences, 119,339 transcripts remained. We refer to these error-corrected, non-redundant transcript isoforms as PacBio-isoforms (Table S1). The PacBio-isoforms showed greater representation of long transcripts and ORFs than uncorrected transcripts (Table 1). The length of the PacBio-isoform sequences ranged from 307 to 26,751 bp, with N50, mean length, and L50 of 3,611 bp, 3,099 bp, and 36,456 sequences, respectively.

Annotation of PacBio-isoforms by comparison with available sequence data

The functions of all PacBio-isoforms were annotated by comparison with various nucleotide and protein databases. BLAST hits in the sense direction to these databases are summarized in Table 2. A total of 107,742 sequences (90.28%) showed sequence similarity to previously described sugarcane nucleotide sequences from SoGI, NCBI EST, and PlantGDB databases. We did not include sugarcane transcript sequences in the SUGIT database (Hoang et al., 2017) for this analysis because no annotations are available for these sequences. A total of 5,026 transcript sequences did not match with any sugarcane sequence in public databases are shown in Table S2. Because sugarcane lacks a genome annotation, plant protein sequences from the Phytozome database were also used for sequence annotation. Protein matches in other plant species were found for 109,053 PacBio-isoforms (91.38%) by BLASTX analysis. The majority (68.03%) of sequences have best hits with sorghum, the species closest to sugarcane, as expected (Fig. 1). When performing the database search to the NCBI nr protein database, 95 sequences showed hits with the NCBI nr database but did not give a hit to any plant database. Most hits overlapped among sugarcane, plant protein, and nr protein databases (Fig. 2). 569 transcripts did not match with any of these three databases (Table S2, Fig. 2). A total of 4,902 sequences did not produce any significant hits with known protein sequences in any protein database. 11 of these sequences hit with known long non-coding RNAs (lncRNAs), suggesting that the remaining 4,891 sequences represent sugarcane proteins of unknown function or new lncRNAs.

Prediction of potential protein coding regions

To identify putative protein-coding sequences within the PacBio-isoforms, we used Transdecoder software to predict open reading frames (ORFs). When considering all types of predicted ORFs of at least 300 bp long, including partial ORFs missing start or stop codon(s), a total of 207,027 ORFs with an N50 of 1,074 bp were present among 104,282 PacBio-isoforms (87.38%) (Table 1). When homology to plant proteins and the presence of protein domain(s) were included as additional evidence for ORF prediction, a total of 255,442 ORFs with an N50 of 951 bp were identified in 110,253 ORF-containing transcripts (92.4%). The ORF N50 value increased to 1,416 bp when only the single best ORF of each transcript was used. Most of ORF-containing sequences (104,337; 94.63%) contain at least one protein motif/domain from various databases as identified by InterProScan. P-loop containing nucleoside triphosphate hydrolase and protein kinase were the most frequently found protein domains (Table S3).

73,795 ORF-containing sequences (67%) contain a complete ORF of at least 300 bp, as indicated by the presence of start and stop codons. The ORF length varied from 300 to 6,657 bp, with the N50 of 1,344 bp. The frequency distribution of the predicted ORFs is displayed in Fig. 3. Flanking these complete ORFs, the mean lengths of 5′ and 3′ untranslated regions are 1,249 and 1,187 bp, respectively (Fig. 3).

9,086 PacBio-isoforms lacked a putative ORF longer than 300 bp. The majority of these sequences (6,558 sequences, 72.18%) have partial match with known protein sequences (from NCBI nr and/or Phytozome plant protein databases) in either sense or antisense direction. Of the remaining 2,528 sequences which did not match with known proteins, 1,219 sequences also produced no hits to Phytozome plant CDS/transcript nucleotide sequences and lack InterPro protein domain, which are suggestive of lncRNAs. Among these potential lncRNAs, only three showed significant similarity to lncRNAs in the lncRNA databases used in this study, suggesting that 1,216 sequences are putatively novel lncRNAs (Table S4). Among these putatively novel lncRNAs (average length = 2,410 bp), 173 did not contain any ORF while the remainder contain short ORFs 15 to 297 bp (median = 189 bp). Most of these transcripts are likely to be of low abundance, with average relative expression levels from Ion Proton data of 1.83 RPKM. 483 sequences did not contain any repetitive sequence or transposable element.

Functional annotation of protein-coding transcripts

A total of 1,382 Gene Ontology (GO) terms covering three main categories (178 cellular component, 649 molecular function, and 555 biological process terms) were assigned to 72,281 PacBio-isoforms (60.57%) based on BLAST results against the Phytozome database. The distribution of GO terms found among PacBio-isoforms is displayed in Table S5. Considering the second-level GO terms, the most frequent terms for the “cellular component” category are cell (19.65%), cell part (19.65%), and membrane (18.89%). For the “molecular function” category, the majority of transcript functions are related to binding (48.07%), catalytic activity (41.78%), and transporter activity (4.70%). The top-ranked “biological process” GO terms are metabolic process (37.68%), cellular process (30.70%), and localization (7.54%). Among third-level GO terms, the most frequent terms for “cellular component”, “molecular function” and “biological process” categories are “membrane”, “protein binding”, and “protein phosphorylation”, respectively.

Transcript sequences were also classified to Cluster of Orthologous Groups (COG) based on groups of orthologous and paralogous proteins among eukaryotes. 43,794 transcripts (36.70%) can be assigned to 26 COG categories (Fig. 4, Table S5). The three COG categories with the highest number of related transcripts are signal transduction mechanism (22.59%), general function prediction only (14.64%), and posttranslation modification, protein turnover, chaperones (9.79%). 13,319 PacBio-isoforms were also mapped to plant-related KEGG pathways based on BLAST result search with the Phytozome database. 205 pathways were mapped (Table S5). The number of mapped transcripts per pathway ranges from 1 to 2,231 transcripts with an average of 160 transcripts per pathway. The top three pathways with the highest number of mapped transcripts are pyruvate metabolism (2,231), spliceosome (1,732), and protein processing in endoplasmic rericulum (1,342). A number of transcripts were also mapped to sugar-related pathways such as starch and sucrose metabolism (1,059), fructose and mannose metabolism (83), and plant-pathogen interaction (532).

Analysis of repetitive elements within transcripts

We investigated repetitive sequences in PacBio-isoforms using RepeatMasker software (Smit, Hubley & Green, 2013–2015). 74,650 sequences (62.55%) contain fragments of repetitive sequences, which include interspersed repeats, small RNA, satellites, simple repeats, and low complexity sequences, while 37,758 transcript sequences (31.64%) contain fragments of interspersed repeats. The total bases masked for repetitive sequences are 6.73%, of which 5.34% correspond to interspersed repeats. 57,483 interspersed repeat elements were found, which included short interspersed nuclear elements (SINEs, 2.09%), long interspersed nuclear elements (LINEs, 13.07%), long terminal repeat elements (LTR, 24.30%), DNA transposons (46.07%), and unclassified elements (14.47%). Considering the proportion of total bases masked, LTR/Gypsy is the most represented class in the transcriptome (21.75%), followed by LTR/Copia (17.32%), LINE/L1 (13.12%), DNA/PIF-Harbinger (6.92%), and DNA/CMC-EnSpm (6.84%).

Assessing completeness of PacBio-isoforms

We assessed the completeness of PacBio-isoforms, i.e., whether they represent full-length transcripts and cover the entire CDS by comparison with related sequences in transcript databases. The length distributions of PacBio-isoforms compared with matched sugarcane and Phytozome plant sequences from BLASTN analysis are illustrated in Fig. 5. The majority (93.26%) of PacBio-isoforms are longer than matched sugarcane sequences reported in SoGI, PlantGDB, and dbEST databases. 78,095 PacBio-isoforms (71.88%) are at least 80% of the length of matched transcripts from sorghum and other plant species. We assessed the proportion of PacBio-isoforms corresponding to putatively full-length transcripts spanning CDS. Out of 107,371 PacBio-isoforms with matches in the sense direction, 68,962 (64.23%) have at least 90% coverage on matched plant CDS, suggesting that these isoforms represent full-length transcripts. Out of these full-length transcripts, 67,188 (97.43%) also contained an ORF of at least 300 bp long; however, the actual number of full-length transcripts could be higher. For sequences which have less than 90% coverage on plant CDS and sequences without plant CDS hits, 38,029 are possibly full-length transcripts of at least 1,000 bp and an ORF of at least 300 bp long.

Comparison with previous Iso-Seq study and identification of novel sugarcane full-length transcripts

We compared PacBio-isoforms with those reported in the SUGIT database (Hoang et al., 2017), which were generated using the same sequencing platform. Our study dataset is comparable in size to the SUGIT database (232,328 and 186,999 full-length non-chimeric reads, respectively), and so provides a comparable level of transcriptomic detail. While Hoang et al. explored a pooled transcriptome from various tissues and cultivars, our study was focused on the leaf tissue transcriptome of the KK3 cultivar. The frequency distributions of transcript length for our PacBio-isoforms and isoforms in the SUGIT database are displayed in Fig. 6. The PacBio-isoforms from both studies cover a broad size range including very long sequences (>7,000 bp). However, the PacBio-isoforms generated in this study are more evenly distributed across all size ranges, whereas the SUGIT PacBio-isoforms are more biased towards the size of around 1,200 bp (Fig. 6). As well as greater representation of long transcripts, our dataset contains more unique high quality transcripts than SUGIT (119,339 vs 107,598) with a greater N50 (3,611 vs 1,991 bp) and a greater number of putative full-length transcripts with at least 90% coverage of matched CDS (67,188 vs 39,999). 114,565 PacBio-isoforms (96%) match 91,683 SUGIT sequences (85.21%). 58.14% of the PacBio-isoforms are longer than the corresponding SUGIT sequence. We investigated the 4,774 PacBio-isoforms with no matching sequence in the SUGIT database in detail (Table S6). Out of these, 4,573 sequences hit with publicly available sequences (sugarcane, plant, or nr database), of which 4,163 sequences also contain an ORF. The 201 unknown sequences with no matching sequence in any database have an average length of 3,229 bp. 86 of these novel sequences contain an ORF with an average length of 3,678 bp. 64 of the ORFs are complete and thus code for putatively novel proteins. The remaining 115 novel sequences with no ORF longer than 300 bp are potentially novel lncRNAs (Tables S6, S4). From functional analysis of the 4,774 PacBio-isoforms not present in SUGIT (Table S6), the top COG categories and the second-level GO terms (“molecular function” and “biological process” categories) are the same as for all PacBio-isoforms, whereas the most frequent term for the “cellular component” category (‘membrane’), is lower-ranked among all PacBio-isoforms. Considering third-level GO terms, the top “biological process”, “molecular function”, and “cellular component” are different from the top terms of all PacBio-isoforms (organic substance metabolic process, heterocyclic compound binding, cell part). The top three pathways (transporters, chromosome and associated proteins, and transfer RNA biogenesis) are also different from all PacBio-isoforms.

Investigation of alternatively spliced transcripts

Although a complete sugarcane genome is not yet available, the sorghum genome can be used as a reference for sugarcane sequence analysis owing to the high gene identity between the two species (Jannoo et al., 2007). 72,061 PacBio-isoforms (60.38%) can be mapped to the sorghum genome, which is similar to the previous Iso-Seq study of sugarcane (69.4%) (Hoang et al., 2017). From BLAST results, 74,184 transcripts (68.03%) matched with sorghum transcripts annotated in the Phytozome database as the best hit. Out of these, 65,279 transcripts (88%) are highly similar to sorghum annotated gene sequences and can be mapped to the same genome location of the sorghum best hit genes. The remaining BLAST hits cannot be mapped to the sorghum genome, possibly because of low similarity to the hit gene, or other highly similar regions in the genome that make unique mapping difficult. The mapping statistics are summarized in Table S7.

PacBio-isoforms mapped to the sorghum genome were used to detect potential alternative splicing events represented among PacBio-isoforms. We identified 17,349 alternative splicing events comprising 5,410 intron retentions (31.2%), 1,423 skipped exons (8.2%), 5,695 alternative 3′ splice sites (32.8%) and 4,821 alternative 5′ splice sites (27.8%). The previous Iso-Seq study of sugarcane reported 4,870 alternative splicing events, in which alternative 3′ splice site was the most common and skipped exon the least common event (Hoang et al., 2017), the same pattern as for our data. We re-analyzed alternative splicing events represented among SUGIT sequences and found patterns of shared and distinct alternative transcripts between the two datasets, examples of which are shown in Fig. 7. For the first example, PacBio-isoforms from each dataset were mapped to the sorghum gene Sobic.002G223200 (peptidase S24/S26A/S26B/S26C family protein) (Fig. 7A). PacBio-isoforms from both datasets showed the same intron retention alternative splicing patterns, confirming the reproducibility of these sugarcane isoforms. Another example is shown for isoforms of the Sobic.004G151800 gene (sucrose-phosphatase) (Fig. 7B). Common alternative transcripts were found as shown by the combined results, whereas alternative transcripts unique to each study were also found.

Analysis of transcriptome expression

The expression levels of PacBio transcript isoforms were explored based on the number of mapped Ion Proton reads which were obtained from the same leaf sample. The majority of transcripts have expression levels about 1-2 RPKM (Fig. S1). The numbers of transcripts with low (≥0.5 to 5 RPKM), medium (≥5 to 50 RPKM), and high (≥50 RPKM) expression levels are 50907, 17185, and 1868 transcripts, respectively. Since the lowly expressed genes are the major population, COG analysis of this group (Fig. S1) showed the same top three categories (signal transduction mechanism, general function prediction only, and posttranslation modification, protein turnover, chaperones) as identified in the analysis of all transcripts (“Functional annotation of protein-coding transcripts” subsection). The medium expressed group showed the similar trend to lowly expressed group with the ranks switching between ‘signal transduction mechanism’ and ‘general function prediction only’. The highly expressed genes tend to have different functions in which the top three COG categories are posttranslation modification, protein turnover, chaperones (26.95%), carbohydrate transport and metabolism (10.68%), and translation, ribosomal structure and biogenesis (9.11%). These highly expressed genes were found to be involved in several pathways including photosynthesis, pyruvate metabolism, nitrogen metabolism, ubiquitin system, glycolysis/gluconeogenesis, and RNA transport (Table S8 ).