BACKGROUND
gene expression profiling is accelerating our progress toward a comprehensive understanding of the genetic mechanisms that control responses to environmental stress.
microarray analysis was developed to obtain overall gene expression profiles in various plants.
microarray profiling and the recently introduced tag-based sequencing approaches are proven technologies for estimating gene expression.
however, array-based technologies have critical limitations  <cit> .
as most microarray probes are designed on the basis of gene annotation, arrays are limited to the analysis of transcripts from previously annotated genes of a sequenced accession of a species.
probes are designed to cover only a very small portion of a gene and so do not represent the whole structure of the gene.
moreover, computationally annotated genes have not fully been validated, because ests and full-length cdnas  cannot cover entire transcribed regions.
it is therefore important to identify whole transcripts  for complete gene expression profiling.
there is a need for the development of technologies beyond arrays.
sequencing-based approaches could overcome the limitations of array-based technologies.
following the rapid progress of massive parallel sequencing technology, whole mrna sequencing has been used for gene expression profiling  <cit> .
this sequencing involves mapping of the reads on known annotated gene models but cannot be used to identify novel genes.
recently, a series of programs have been developed for building gene models directly from the piling up of short reads: bowtie efficiently maps short reads on genomic sequences  <cit> ; tophat concatenates adjacent exons and identifies reads that bridge exon junctions  <cit> ; and cufflinks  <cit>  constructs gene models from the exons and bridging sequences predicted by bowtie and tophat and then calculates their abundances of these sequences.
the use of this series of programs has the potential to discover new transcripts from mrna-seq  but has only just begun  <cit> .
in this study, we identified unannotated transcripts in rice on the basis of the piling up of mapped reads.
as a model case, we give examples of salinity stress-inducible unannotated transcripts encoding putative functional proteins.
for these purposes, we performed whole mrna sequencing by using massive parallel sequencing technology.
we also took advantage of various high-quality genomic resources in rice, including the genomic sequence , fl-cdna sequences  <cit> , the rice annotation project database   <cit> , and a rice 44k microarray , in our analysis of rice transcriptomes.
first, to estimate the scale of the transcriptomes in rice, we mapped 36-base-pair  reads from the mrna of salinity stress-treated rice tissues on the rice genome.
the coverage of previously annotated regions or of the rice genome was then calculated.
second, we attempted to identify salinity stress-inducible genes as a model system for gene expression profiling by mrna-seq.
the number of mapped reads was counted and marked on the rice genome.
third, using the mrna-seq data, we used bowtie, tophat, and cufflinks to construct gene models based on the piling up of short reads on the rice genome, and compared these with previous annotations and then characterized the unannotated transcripts.
we conducted a blastx search for the unannotated transcripts, and we discuss the function of the encoded proteins.
fourth, to validate our sequence-based technology, we compared the results of quantification by the array-based and sequence-based approaches, and we discuss the advantages of the latter.
this work contributes to the discovery of whole salinity stress-inducible transcripts without the need to rely on previous annotations.
it should help to establish further sequence-based gene expression profiling in any organism.
RESULTS
mapping of 36-bp reads to the rice genome
we performed rice transcriptome analysis at single-nucleotide resolution by using illumina mrna-seq technology.
briefly, poly rnas from salinity stress-treated rice tissues were reverse-transcribed and sequenced .
millions of 36-bp reads were mapped to the rice genomic sequence , with at most two mismatches or  <dig> bp of indels allowed.
to obtain many kinds of transcripts, data on nine technical replicates of the sequencing run of cdna from the roots after salinity stress were accumulated.
as the number of reads increased, the cumulative coverage of both the genome and the annotated transcribed region gradually approached a plateau .
saturation of sequencing was also estimated on the basis of the fraction of genes that had reached their final rpkm   <cit> .
as the number of reads increased, the fraction of highly expressed genes  close to their final rpkm was almost unchanged, whereas those of genes with relatively low expression  converged more slowly .
with four technical replicates ,  <dig> % of genes with relatively low expression levels  reached to within ± 5% of their final rpkm .
thus, for further analysis, we adopted the summing of four technical replicates after filtration according to their base quality.
rnas were prepared from normal shoot , normal root , shoot with salinity stress , or root with salinity stress .
numbers of 36-bp reads and their percentages of the total number of filtered reads, obtained by the summing of four technical replicates , are shown, as are reads mapped to the rice genome , reads mapped uniquely to a predicted exon-exon bridging sequence , the total number of reads mapped uniquely to the genome and to a predicted exon-exon bridging sequence , reads mapped to multiple loci of the rice genome , and reads unable to be mapped to the rice genome .
rice transcriptome analysis was based on response to salinity stress.
mrnas were prepared from the tissues of normal rice shoots and roots and from those subjected to  <dig> h of salinity stress.
of the  <dig> to  <dig> million quality-evaluated reads ,  <dig> % to  <dig> % were mapped uniquely to the rice genome ;  <dig> % to  <dig> % of the reads bridged flanking exons ;  <dig> % to  <dig> % of the reads were repetitive sequences ; and  <dig> % to  <dig> % had no match in the genome .
thus, a total of  <dig> % to  <dig> % of the reads were mapped uniquely to the rice genome or to exon-exon junctions .
of the unmapped reads,  <dig> % had high levels of identity to sequences derived from sequencing adaptors , contaminating organisms , or ribosomal rna  .
a few transcripts might have been transcribed from unsequenced genomic regions of rice  <cit> .
however, most of the unmapped reads  had no similarity to each other .
our preliminary experiment showed that the ratio of these unmapped reads was higher with mrna-seq  than with genomic sequencing .
thus, part of the random sequences might have come from residual random primers used in cdna synthesis.
the common random sequences might have come from sequencing errors in the use of the illumina sequencing technology.
identification of differentially expressed genes by mrna-seq
mrna-seq quantifies the amount of transcripts on the basis of the number of sequence reads mapped on each gene.
we adopted this method for transcript quantification by rpkm  <cit>  and calculated the rpkm of each gene .
rpkm quantification was distributed from  <dig> to over  <dig>  in shoots under normal conditions, the gene encoding ribulose bisphosphate carboxylase activase  was expressed at extremely high levels .
in roots under normal conditions, the gene for metallothionein  was expressed at extremely high levels .
the statistical mean and median were  <dig>  and  <dig> , respectively, in the shoot, and  <dig>  and  <dig>  in the root under normal conditions.
we then comprehensively compared the rpkm of each gene in response to salinity stress .
we used the g-test with a 1% false discovery rate  and identified  <dig>   and  <dig>   differentially expressed rap <dig> genes.
of these,  <dig>   genes were commonly differentially expressed.
the number of highly differentially expressed genes , such as those encoding bhlh-containing protein  and amino acid transporter , was greater in the root  than in the shoot .
expression of genes previously identified under salinity stress  <cit> --i.e.
ostpp <dig> , lip <dig> , osaba <dig> , osmst <dig> , wsi <dig> , and mybs <dig> --was induced in the root .
for a complete comparison see additional file 2: table s <dig>
the distribution of mapped reads on the rice genome was graphed on a gbrowse  <cit>  .
for example, the ostpp <dig>  gene , which encodes a protein that synthesizes the abiotic stress-protectant trehalose  <cit> , was expressed exclusively in the root after  <dig> h of salinity stress; rcc <dig> , which was previously identified as a root-specific gene  <cit> , was expressed only in the root with and without stress; ak <dig>  was expressed exclusively in the shoot; most of the neighboring genes were expressed evenly in all tissues used .
constructing gene models by mrna-seq
transcribed regions were identified on the basis of the piling up of mapped short reads through the programs bowtie  <cit> , tophat  <cit> , and cufflinks  <cit> .
in the shoot,  <dig>  transcripts were predicted  ;  <dig> %  of the predicted transcripts were mapped on previously annotated loci in rap <dig>  <cit> ; thus, the remaining  <dig>  predicted transcripts were unannotated in rap-db .
in the root,  <dig>  of the  <dig>  predicted transcripts were mapped on unannotated regions .
for example, the previously annotated gene ak <dig>  which is similar to dreb1b in arabidopsis thaliana , was expressed after salinity stress and also predicted by cufflinks ; other exons were also predicted and connected by bridging sequences elucidated by tophat  .
reads were also mapped on the extended parts of the ends of most 5' and 3' exons in previous gene models .
of the transcripts mapped on previously annotated loci,  <dig>   and  <dig>   had not been supported by ests  <cit>  or fl-cdnas  <cit> .
transcripts were predicted on the basis of the piling up of mapped reads by the cufflinks program .
orf prediction was conducted by blastx search against uniprot  and refseq .
orf: open reading frame; a.a: amino acids.
we attempted to predict the functions of unannotated transcripts by blastx search and longest-orf search.
in a blastx search against the uniprot and refseq sequences, of the predicted transcripts,  <dig>  and  <dig>   had orfs similar to those encoding the amino acid sequences of functional proteins .
of the remaining unannotated transcripts,  <dig>   and  <dig>   had orfs encoding at least  <dig> amino acids by longest-orf search .
amino acid length was widely distributed: the mean and median were  <dig> and  <dig> amino acids in the shoot, and  <dig> and  <dig> in the root .
we used the g-test with a 1% fdr and identified  <dig>  and  <dig>  differentially expressed cufflinks transcripts.
even though the lengths of cufflinks transcripts were not completely identical between shoot and root, at least  <dig> differentially expressed transcripts were common to the two tissues.
in response to salinity stress,  <dig>  and  <dig>  unannotated transcripts were upregulated  .
these unannotated transcripts encoded, for example, proteins similar to indole-3-glycerol phosphate lyase and gibberellin 2-beta-dioxygenase .
of the other differentially expressed genes , root_cuff.
<dig>  was upregulated ; it encoded proteins similar to msl <dig>  .
for a complete list of unannotated transcripts see additional file 3: table s <dig>
cufflinks id ; total nucleotide length of each predicted transcript ; rpkm without salinity stress  or with salinity stress ; calculated ratio of rpkm ; number of amino acids encoded by putative orf ; and name of similar protein  are listed.
differentially expressed cufflinks transcripts were identified by the g-test with a 1% false discovery rate.
highly differentially expressed genes  derived from different loci that had orfs predicted by blastx search are listed.
"-" means a calculated ratio of infinity because the rpkm without salinity stress  was  <dig>  detailed data for all unannotated transcripts are listed in additional file 3: table s <dig>
comparison of sequence-based and array-based technologies for gene expression profiling
our sequence-based gene expression profiling was validated against array-based technology.
first, signal intensity and rpkm from the same rna materials were compared.
these two independent measures of transcript abundance were correlated , especially at moderately high signal intensities .
however, the correlation was not as strong at extremely high signal intensities , suggesting that the array signal intensity was saturated but the rpkm was not .
next, the ratios of differentially expressed genes were compared.
the ratio obtained from the array and the corresponding ratio obtained from rpkm was highly correlated over a broad range .
the histogram was highest at log <dig> , suggesting that most genes were expressed evenly both before and  <dig> h after salinity stress .
however, a few discrepancies were found: increased changes in the expression of  <dig> genes were found by using the array , but not by using mrna-seq ; conversely, increased changes in the expression of  <dig> genes were found by using mrna-seq , but not by using the array  .
to further examine these discrepancies, we used quantitative real-time polymerase chain reaction .
the qrt-pcr results suggested that most of the former discrepancy was due to technical inaccuracy in the array experiments.
however, qrt-pcr supported only three of the seven mrna-seq data in the latter discrepancy .
despite these discrepancies, our sequence-based approach was generally valid as a gene expression profiling technology for use with previously annotated genes.
discussion
estimation of variation and abundance of whole transcripts in rice
how many reads are required to cover whole transcripts in the rice cell?
as the number of reads increased, the cumulative coverage approached a plateau .
we summed four technical replicates .
rpkm is widely used to calculate the abundance of each transcript and is linear across a dynamic range  <cit> .
the distribution of rpkm of rice genes ranged from  <dig> to over  <dig> ; genes involved in photosynthesis in the shoot or in regulation of physiological metals in the root were highly expressed, whereas about 30% of genes had rpkm <  <dig> .
the saturation of sequencing in rice  was almost the same as in a previous mammalian analysis  <cit> .
according to that analysis, "one transcript in a cell corresponds to  <dig> to  <dig> rpkm"  <cit> , so genes having rpkm <  <dig> might rarely be expressed.
however, data on the rna content of each rice cell are required to calculate the number of existing molecules of rnas.
as rice tissue contains cells of various sizes and types, the relationship between the number of existing molecules and their rpkm has not yet been accurately determined.
when we used four technical replicates, about 20% of genes expressed at relatively low levels  did not reach their final rpkm , suggesting that these model settings were insufficient for calculating the real rpkm of genes expressed at low levels.
summing of the four technical replicates covered  <dig> % of all annotated regions, corresponding to  <dig> % of  <dig> mb  <cit>  of the rice genome .
this result suggests that these regions were transcriptionally active under the experimental conditions.
even though the cumulative coverage was close to a plateau, the coverage rose gradually; the accumulation of about  <dig> million reads covered  <dig> % of annotated regions , suggesting that some of the reads expressed at low levels were not sequenced.
however, the gradual increase in coverage might have been due to the presence of contaminated genomic dna or a very small amount of partly processed nuclear rnas, because intron retention is the most prevalent alternative splicing form in rice  <cit> , as it is in arabidopsis thaliana  <cit> .
thus, we consider that the summing of four technical replicates of 36-bp reads, corresponding to a total of  <dig> gbp of filtered sequences, covered almost all the transcripts in the rice cell under the experimental conditions, although more reads are required to obtain the final rpkm of genes expressed at relatively low levels.
identification of unannotated transcripts by mrna sequencing
mrna-seq provides information on whole transcribed genes without the need to rely on annotation , whereas array technology is limited to providing data only on those previously annotated genes and on previously identified ests with no known homologies that have corresponding probes on the array.
on the basis of the piling up of mapped reads, we predicted  <dig>   and  <dig>   currently unannotated transcripts in rap-db .
of the rap <dig> unannotated transcripts,  <dig> %  in shoot and  <dig> %  in root had not been annotated by michigan state university  , suggesting that these transcripts were novel transcripts.
unannotated transcripts included extended parts of previously annotated genes .
extension of 5' exons might contribute to the making of a different start codon or the shifting of the reading frame of previously annotated genes.
extension of 3' utrs might contribute to microrna-mediated control of translation or post-transcriptional rna metabolism  <cit> .
for example, mrna-seq provided evidence of the existence of extended parts of previously annotated genes and of the differential regulation of their expression.
ak <dig>  previously annotated as a non-protein-coding transcript, had additional predicted exons distal to the 5' end of the previous gene model, and it encoded an indole-3-glycerol phosphate lyase .
two neighboring genes  were also similar to the indole-3-glycerol phosphate lyase gene, suggesting that all three genes were tandemly duplicated.
although all three genes were upregulated in response to salinity stress, their tissue specificities and expression levels differed , suggesting that their functions diversified after gene duplication.
mrna-seq also provided evidence of expression of computationally predicted genes.
the existence of a number of genes computationally predicted in rap-db  <cit>  has not been supported  <cit>  by ests  <cit>  or fl-cdnas  <cit> .
here,  <dig>   and  <dig>   transcripts identified by mrna-seq have been mapped on computationally predicted genes, the presence of which was not supported by experiments, suggesting the validity of the computationally predicted gene models in rap-db.
we will use these sequence-based transcriptome analyses to improve rap-db.
mrna-seq provided details of the bridging sequences between exons, suggesting the presence of splicing junctions, whereas array technology--including whole-genome tiling arrays  <cit> --provides no information on connecting exons.
because reads that bridge exon boundaries are not mapped directly to the genomic sequence, a mapping technique was required.
as a first step, the enumeration of all theoretical splicing junctions within annotated transcripts allows the mapping of bridging reads  <cit>  by using statistical models  <cit> .
we found that  <dig> % to  <dig> % of reads formed primary bridges with previously annotated exons ; this was not sufficient to discover sequences bridging unannotated transcripts.
programs such as tophat  <cit>  and g-mo.
r-se   <cit>  are designed to align reads to form potential splice junctions without relying on known splice sites.
in this study, sequences flanking potential donor/acceptor splice sites were joined to form canonical  introns between neighboring  islands by using tophat  <cit> .
even though we used tophat for our prediction, some of the predicted transcripts remained to be separated--unlike the case with the fl-cdna sequences--because of the lack of sufficient bridging sequences between the exons , suggesting that more bridging reads should be sequenced to connect predicted exons.
elongation of the length of each read may also enhance the chance to connect predicted exons.
sequence-based transcriptome analysis for capturing salinity stress-inducible genes in rice
mrna-seq comprehensively identified salinity stress-inducible genes.
unannotated transcripts had orfs  with a mean length of  <dig> amino acids  or  <dig> amino acids  , suggesting that these unannotated transcripts could encode functional proteins.
of the unannotated transcripts,  <dig>  and  <dig>  were differentially expressed in response to salinity stress .
these unannotated transcripts encoded proteins associated with functions such as amino acid metabolism  in response to abiotic stress  <cit> , diterpenoid biosynthesis , and mechanosensitive ion channel  function  <cit> .
mechanosensitive ion channels are gated directly by physical stimuli such as osmotic shock and transduce these stimuli into electrical signals  <cit> .
mrna-seq also captured previously identified genes involved in salinity tolerance, namely those associated with trehalose synthesis  , dehydrin , aba synthesis , sugar transport , glycerol transferase , and transcription factors similar to those of the dreb family .
a substantial number of transcripts were exclusively upregulated only in the root .
as only the root was directly exposed to  <dig> h of salinity stress, it might take time to induce the expression of more genes in the shoot; ostpp <dig>  might be expressed in the shoot after  <dig> h of exposure, as has been found in yukihikari rice  <cit> .
with these genes, nipponbare may have the potential to be tolerant to salinity stress.
rice cultivars such as nona bokra and pokkali are substantially more salinity tolerant than nipponbare  <cit> , suggesting that the genuine salinity stress tolerance gene might be missing in nipponbare.
the  <dig> oryza species are geographically, physiologically, and genetically diverse  <cit> , and many of the genes in cultivated rices have been selected by humans under field conditions, not by environmental stress.
these essentially missing genes could serve as potential genetic resources for the improvement of cultivated crops.
sequence-based technology can be used to extract such missing genes by the piling-up of short reads on their own genomes without the need to rely on sequence similarity.
overcoming the technical inaccuracy
microarray technology has been used as a sophisticated platform for the expression profiling of previously annotated genes.
however, as an array-based technology, evaluation of signal intensities close to background levels tends to cause artifacts in array analysis because of high levels of background noise and/or cross-hybridization  <cit> ; moreover, hybridization efficiency might vary with the probes used, suggesting that the calculation of real molar concentrations is inaccurate.
whereas the agilent rice 44k array is designed to quantify 60-mer sequences at the 3'-end of transcripts, mrna-seq quantifies transcript abundance on the basis of the number of mapped sequences on the whole gene model.
in our study, the two measures of transcript abundance  and change ratios  were highly correlated, as in a previous report  <cit> .
moreover, for genes expressed at low or extremely high levels  and for genes differentially expressed in arrays , mrna-seq seemed to be accurate.
therefore, mrna-seq measures the molar concentrations of genes accurately over a broad dynamic range.
biological replication is required for identifying differentially expressed genes through statistical analysis, as in array-based analysis.
unfortunately, with sequence-based transcriptome analysis there are greater costs than with microarrays for cdna preparation and sequencing; this prevented us from performing further experiments.
illumina has improved its sequencing technology.
each read length has been continuously increased.
efficient base calling by using the latest illumina data analysis pipeline software improved the quality and quantity of reads from the same raw image data.
controlled hydrolysis of rna before cdna synthesis substantially improved the uniformity of sequence coverage, as in a previous report  <cit> .
these technical innovations in hardware and software will enable remarkable progress in reducing costs and in increasing the sensitivity of detection of sequences transcribed at low levels, the accuracy of quantification and detection of splice forms, and the prediction of the whole structures of transcripts.
sequence-based transcriptome analysis has recently been applied to various organisms: arabidopsis thaliana  <cit> , yeasts  <cit> , drosophila melanogaster  <cit> , and human  <cit> .
during this study, two types of rice transcriptome analysis were reported, focusing on the transcriptional differences in two rice subspecies and their reciprocal hybrids  <cit>  and in eight organs from different developmental stages of oryza sativa l. ssp.
indica '93-11'  <cit> .
we analyzed salinity stress-inducible transcripts and constructed gene models based on the pilling up of short reads by using the cufflinks program.
this approach should help to discover novel gene models without reliance on gene annotation.
CONCLUSIONS
microarray-based gene expression profiling is limited to the analysis of annotated genes.
in our mrna-seq analysis, unannotated salinity stress-inducible transcripts were identified on the basis of the piling up of mapped reads without reliance on gene annotation or fl-cdna sequences.
some of these novel transcripts had orfs encoding putative functional proteins and were differentially expressed in response to salinity stress.
mrna-seq was valid as a gene expression profiling technology for quantifying the abundance of previously annotated genes.
our findings will contribute to improvement of our rap-db and to further sequence-based gene expression profiling in various organisms.

 gene expression profiling is accelerating our progress toward a comprehensive understanding of the genetic mechanisms that control responses to environmental stress. microarray analysis was developed to obtain overall gene expression profiles in various plants. microarray profiling and the recently introduced tag-based sequencing approaches are proven technologies for estimating gene expression. however, array-based technologies have critical limitations. as most microarray probes are designed on the basis of gene annotation, arrays are limited to the analysis of transcripts from previously annotated genes of a sequenced accession of a species. probes are designed to cover only a very small portion of a gene and so do not represent the whole structure of the gene. moreover, computationally annotated genes have not fully been validated, because ests and full-length cdnas cannot cover entire transcribed regions. it is therefore important to identify whole transcripts for complete gene expression profiling. there is a need for the development of technologies beyond arrays. sequencing-based approaches could overcome the limitations of array-based technologies. following the rapid progress of massive parallel sequencing technology, whole mrna sequencing has been used for gene expression profiling. this sequencing involves mapping of the reads on known annotated gene models but cannot be used to identify novel genes. recently, a series of programs have been developed for building gene models directly from the piling up of short reads: bowtie efficiently maps short reads on genomic sequences; tophat concatenates adjacent exons and identifies reads that bridge exon junctions; and cufflinks constructs gene models from the exons and bridging sequences predicted by bowtie and tophat and then calculates their abundances of these sequences. the use of this series of programs has the potential to discover new transcripts from mrna-seq but has only just begun. in this study, we identified unannotated transcripts in rice on the basis of the piling up of mapped reads. as a model case, we give examples of salinity stress-inducible unannotated transcripts encoding putative functional proteins. for these purposes, we performed whole mrna sequencing by using massive parallel sequencing technology. we also took advantage of various high-quality genomic resources in rice, including the genomic sequence, fl-cdna sequences, the rice annotation project database, and a rice 44k microarray, in our analysis of rice transcriptomes. first, to estimate the scale of the transcriptomes in rice, we mapped 36-base-pair reads from the mrna of salinity stress-treated rice tissues on the rice genome. the coverage of previously annotated regions or of the rice genome was then calculated. second, we attempted to identify salinity stress-inducible genes as a model system for gene expression profiling by mrna-seq. the number of mapped reads was counted and marked on the rice genome. third, using the mrna-seq data, we used bowtie, tophat, and cufflinks to construct gene models based on the piling up of short reads on the rice genome, and compared these with previous annotations and then characterized the unannotated transcripts. we conducted a blastx search for the unannotated transcripts, and we discuss the function of the encoded proteins. fourth, to validate our sequence-based technology, we compared the results of quantification by the array-based and sequence-based approaches, and we discuss the advantages of the latter. this work contributes to the discovery of whole salinity stress-inducible transcripts without the need to rely on previous annotations. it should help to establish further sequence-based gene expression profiling in any organism.  mapping of 36-bp reads to the rice genome we performed rice transcriptome analysis at single-nucleotide resolution by using illumina mrna-seq technology. briefly, poly rnas from salinity stress-treated rice tissues were reverse-transcribed and sequenced. millions of 36-bp reads were mapped to the rice genomic sequence, with at most two mismatches or bp of indels allowed. to obtain many kinds of transcripts, data on nine technical replicates of the sequencing run of cdna from the roots after salinity stress were accumulated. as the number of reads increased, the cumulative coverage of both the genome and the annotated transcribed region gradually approached a plateau. saturation of sequencing was also estimated on the basis of the fraction of genes that had reached their final rpkm. as the number of reads increased, the fraction of highly expressed genes close to their final rpkm was almost unchanged, whereas those of genes with relatively low expression converged more slowly. with four technical replicates, % of genes with relatively low expression levels reached to within ± 5% of their final rpkm. thus, for further analysis, we adopted the summing of four technical replicates after filtration according to their base quality. rnas were prepared from normal shoot, normal root, shoot with salinity stress, or root with salinity stress. numbers of 36-bp reads and their percentages of the total number of filtered reads, obtained by the summing of four technical replicates, are shown, as are reads mapped to the rice genome, reads mapped uniquely to a predicted exon-exon bridging sequence, the total number of reads mapped uniquely to the genome and to a predicted exon-exon bridging sequence, reads mapped to multiple loci of the rice genome, and reads unable to be mapped to the rice genome. rice transcriptome analysis was based on response to salinity stress. mrnas were prepared from the tissues of normal rice shoots and roots and from those subjected to h of salinity stress. of the to million quality-evaluated reads, % to % were mapped uniquely to the rice genome; % to % of the reads bridged flanking exons; % to % of the reads were repetitive sequences; and % to % had no match in the genome. thus, a total of % to % of the reads were mapped uniquely to the rice genome or to exon-exon junctions. of the unmapped reads, % had high levels of identity to sequences derived from sequencing adaptors, contaminating organisms, or ribosomal rna. a few transcripts might have been transcribed from unsequenced genomic regions of rice. however, most of the unmapped reads had no similarity to each other. our preliminary experiment showed that the ratio of these unmapped reads was higher with mrna-seq than with genomic sequencing. thus, part of the random sequences might have come from residual random primers used in cdna synthesis. the common random sequences might have come from sequencing errors in the use of the illumina sequencing technology. identification of differentially expressed genes by mrna-seq mrna-seq quantifies the amount of transcripts on the basis of the number of sequence reads mapped on each gene. we adopted this method for transcript quantification by rpkm and calculated the rpkm of each gene. rpkm quantification was distributed from to over in shoots under normal conditions, the gene encoding ribulose bisphosphate carboxylase activase was expressed at extremely high levels. in roots under normal conditions, the gene for metallothionein was expressed at extremely high levels. the statistical mean and median were and, respectively, in the shoot, and and in the root under normal conditions. we then comprehensively compared the rpkm of each gene in response to salinity stress. we used the g-test with a 1% false discovery rate and identified and differentially expressed rap genes. of these, genes were commonly differentially expressed. the number of highly differentially expressed genes, such as those encoding bhlh-containing protein and amino acid transporter, was greater in the root than in the shoot. expression of genes previously identified under salinity stress --i.e. ostpp, lip, osaba, osmst, wsi, and mybs --was induced in the root. for a complete comparison see additional file 2: table s the distribution of mapped reads on the rice genome was graphed on a gbrowse. for example, the ostpp gene, which encodes a protein that synthesizes the abiotic stress-protectant trehalose, was expressed exclusively in the root after h of salinity stress; rcc, which was previously identified as a root-specific gene, was expressed only in the root with and without stress; ak was expressed exclusively in the shoot; most of the neighboring genes were expressed evenly in all tissues used. constructing gene models by mrna-seq transcribed regions were identified on the basis of the piling up of mapped short reads through the programs bowtie, tophat, and cufflinks. in the shoot, transcripts were predicted; % of the predicted transcripts were mapped on previously annotated loci in rap; thus, the remaining predicted transcripts were unannotated in rap-db. in the root, of the predicted transcripts were mapped on unannotated regions. for example, the previously annotated gene ak which is similar to dreb1b in arabidopsis thaliana, was expressed after salinity stress and also predicted by cufflinks; other exons were also predicted and connected by bridging sequences elucidated by tophat. reads were also mapped on the extended parts of the ends of most 5' and 3' exons in previous gene models. of the transcripts mapped on previously annotated loci, and had not been supported by ests or fl-cdnas. transcripts were predicted on the basis of the piling up of mapped reads by the cufflinks program. orf prediction was conducted by blastx search against uniprot and refseq. orf: open reading frame; a.a: amino acids. we attempted to predict the functions of unannotated transcripts by blastx search and longest-orf search. in a blastx search against the uniprot and refseq sequences, of the predicted transcripts, and had orfs similar to those encoding the amino acid sequences of functional proteins. of the remaining unannotated transcripts, and had orfs encoding at least amino acids by longest-orf search. amino acid length was widely distributed: the mean and median were and amino acids in the shoot, and and in the root. we used the g-test with a 1% fdr and identified and differentially expressed cufflinks transcripts. even though the lengths of cufflinks transcripts were not completely identical between shoot and root, at least differentially expressed transcripts were common to the two tissues. in response to salinity stress, and unannotated transcripts were upregulated. these unannotated transcripts encoded, for example, proteins similar to indole-3-glycerol phosphate lyase and gibberellin 2-beta-dioxygenase. of the other differentially expressed genes, root_cuff. was upregulated; it encoded proteins similar to msl. for a complete list of unannotated transcripts see additional file 3: table s cufflinks id; total nucleotide length of each predicted transcript; rpkm without salinity stress or with salinity stress; calculated ratio of rpkm; number of amino acids encoded by putative orf; and name of similar protein are listed. differentially expressed cufflinks transcripts were identified by the g-test with a 1% false discovery rate. highly differentially expressed genes derived from different loci that had orfs predicted by blastx search are listed. "-" means a calculated ratio of infinity because the rpkm without salinity stress was detailed data for all unannotated transcripts are listed in additional file 3: table s comparison of sequence-based and array-based technologies for gene expression profiling our sequence-based gene expression profiling was validated against array-based technology. first, signal intensity and rpkm from the same rna materials were compared. these two independent measures of transcript abundance were correlated, especially at moderately high signal intensities. however, the correlation was not as strong at extremely high signal intensities, suggesting that the array signal intensity was saturated but the rpkm was not. next, the ratios of differentially expressed genes were compared. the ratio obtained from the array and the corresponding ratio obtained from rpkm was highly correlated over a broad range. the histogram was highest at log, suggesting that most genes were expressed evenly both before and h after salinity stress. however, a few discrepancies were found: increased changes in the expression of genes were found by using the array, but not by using mrna-seq; conversely, increased changes in the expression of genes were found by using mrna-seq, but not by using the array. to further examine these discrepancies, we used quantitative real-time polymerase chain reaction. the qrt-pcr results suggested that most of the former discrepancy was due to technical inaccuracy in the array experiments. however, qrt-pcr supported only three of the seven mrna-seq data in the latter discrepancy. despite these discrepancies, our sequence-based approach was generally valid as a gene expression profiling technology for use with previously annotated genes. discussion estimation of variation and abundance of whole transcripts in rice how many reads are required to cover whole transcripts in the rice cell? as the number of reads increased, the cumulative coverage approached a plateau. we summed four technical replicates. rpkm is widely used to calculate the abundance of each transcript and is linear across a dynamic range. the distribution of rpkm of rice genes ranged from to over; genes involved in photosynthesis in the shoot or in regulation of physiological metals in the root were highly expressed, whereas about 30% of genes had rpkm <. the saturation of sequencing in rice was almost the same as in a previous mammalian analysis. according to that analysis, "one transcript in a cell corresponds to to rpkm", so genes having rpkm < might rarely be expressed. however, data on the rna content of each rice cell are required to calculate the number of existing molecules of rnas. as rice tissue contains cells of various sizes and types, the relationship between the number of existing molecules and their rpkm has not yet been accurately determined. when we used four technical replicates, about 20% of genes expressed at relatively low levels did not reach their final rpkm, suggesting that these model settings were insufficient for calculating the real rpkm of genes expressed at low levels. summing of the four technical replicates covered % of all annotated regions, corresponding to % of mb of the rice genome. this result suggests that these regions were transcriptionally active under the experimental conditions. even though the cumulative coverage was close to a plateau, the coverage rose gradually; the accumulation of about million reads covered % of annotated regions, suggesting that some of the reads expressed at low levels were not sequenced. however, the gradual increase in coverage might have been due to the presence of contaminated genomic dna or a very small amount of partly processed nuclear rnas, because intron retention is the most prevalent alternative splicing form in rice, as it is in arabidopsis thaliana. thus, we consider that the summing of four technical replicates of 36-bp reads, corresponding to a total of gbp of filtered sequences, covered almost all the transcripts in the rice cell under the experimental conditions, although more reads are required to obtain the final rpkm of genes expressed at relatively low levels. identification of unannotated transcripts by mrna sequencing mrna-seq provides information on whole transcribed genes without the need to rely on annotation, whereas array technology is limited to providing data only on those previously annotated genes and on previously identified ests with no known homologies that have corresponding probes on the array. on the basis of the piling up of mapped reads, we predicted and currently unannotated transcripts in rap-db. of the rap unannotated transcripts, % in shoot and % in root had not been annotated by michigan state university, suggesting that these transcripts were novel transcripts. unannotated transcripts included extended parts of previously annotated genes. extension of 5' exons might contribute to the making of a different start codon or the shifting of the reading frame of previously annotated genes. extension of 3' utrs might contribute to microrna-mediated control of translation or post-transcriptional rna metabolism. for example, mrna-seq provided evidence of the existence of extended parts of previously annotated genes and of the differential regulation of their expression. ak previously annotated as a non-protein-coding transcript, had additional predicted exons distal to the 5' end of the previous gene model, and it encoded an indole-3-glycerol phosphate lyase. two neighboring genes were also similar to the indole-3-glycerol phosphate lyase gene, suggesting that all three genes were tandemly duplicated. although all three genes were upregulated in response to salinity stress, their tissue specificities and expression levels differed, suggesting that their functions diversified after gene duplication. mrna-seq also provided evidence of expression of computationally predicted genes. the existence of a number of genes computationally predicted in rap-db has not been supported by ests or fl-cdnas. here, and transcripts identified by mrna-seq have been mapped on computationally predicted genes, the presence of which was not supported by experiments, suggesting the validity of the computationally predicted gene models in rap-db. we will use these sequence-based transcriptome analyses to improve rap-db. mrna-seq provided details of the bridging sequences between exons, suggesting the presence of splicing junctions, whereas array technology--including whole-genome tiling arrays --provides no information on connecting exons. because reads that bridge exon boundaries are not mapped directly to the genomic sequence, a mapping technique was required. as a first step, the enumeration of all theoretical splicing junctions within annotated transcripts allows the mapping of bridging reads by using statistical models. we found that % to % of reads formed primary bridges with previously annotated exons; this was not sufficient to discover sequences bridging unannotated transcripts. programs such as tophat and g-mo. r-se are designed to align reads to form potential splice junctions without relying on known splice sites. in this study, sequences flanking potential donor/acceptor splice sites were joined to form canonical introns between neighboring islands by using tophat. even though we used tophat for our prediction, some of the predicted transcripts remained to be separated--unlike the case with the fl-cdna sequences--because of the lack of sufficient bridging sequences between the exons, suggesting that more bridging reads should be sequenced to connect predicted exons. elongation of the length of each read may also enhance the chance to connect predicted exons. sequence-based transcriptome analysis for capturing salinity stress-inducible genes in rice mrna-seq comprehensively identified salinity stress-inducible genes. unannotated transcripts had orfs with a mean length of amino acids or amino acids, suggesting that these unannotated transcripts could encode functional proteins. of the unannotated transcripts, and were differentially expressed in response to salinity stress. these unannotated transcripts encoded proteins associated with functions such as amino acid metabolism in response to abiotic stress, diterpenoid biosynthesis, and mechanosensitive ion channel function. mechanosensitive ion channels are gated directly by physical stimuli such as osmotic shock and transduce these stimuli into electrical signals. mrna-seq also captured previously identified genes involved in salinity tolerance, namely those associated with trehalose synthesis, dehydrin, aba synthesis, sugar transport, glycerol transferase, and transcription factors similar to those of the dreb family. a substantial number of transcripts were exclusively upregulated only in the root. as only the root was directly exposed to h of salinity stress, it might take time to induce the expression of more genes in the shoot; ostpp might be expressed in the shoot after h of exposure, as has been found in yukihikari rice. with these genes, nipponbare may have the potential to be tolerant to salinity stress. rice cultivars such as nona bokra and pokkali are substantially more salinity tolerant than nipponbare, suggesting that the genuine salinity stress tolerance gene might be missing in nipponbare. the oryza species are geographically, physiologically, and genetically diverse, and many of the genes in cultivated rices have been selected by humans under field conditions, not by environmental stress. these essentially missing genes could serve as potential genetic resources for the improvement of cultivated crops. sequence-based technology can be used to extract such missing genes by the piling-up of short reads on their own genomes without the need to rely on sequence similarity. overcoming the technical inaccuracy microarray technology has been used as a sophisticated platform for the expression profiling of previously annotated genes. however, as an array-based technology, evaluation of signal intensities close to background levels tends to cause artifacts in array analysis because of high levels of background noise and/or cross-hybridization; moreover, hybridization efficiency might vary with the probes used, suggesting that the calculation of real molar concentrations is inaccurate. whereas the agilent rice 44k array is designed to quantify 60-mer sequences at the 3'-end of transcripts, mrna-seq quantifies transcript abundance on the basis of the number of mapped sequences on the whole gene model. in our study, the two measures of transcript abundance and change ratios were highly correlated, as in a previous report. moreover, for genes expressed at low or extremely high levels and for genes differentially expressed in arrays, mrna-seq seemed to be accurate. therefore, mrna-seq measures the molar concentrations of genes accurately over a broad dynamic range. biological replication is required for identifying differentially expressed genes through statistical analysis, as in array-based analysis. unfortunately, with sequence-based transcriptome analysis there are greater costs than with microarrays for cdna preparation and sequencing; this prevented us from performing further experiments. illumina has improved its sequencing technology. each read length has been continuously increased. efficient base calling by using the latest illumina data analysis pipeline software improved the quality and quantity of reads from the same raw image data. controlled hydrolysis of rna before cdna synthesis substantially improved the uniformity of sequence coverage, as in a previous report. these technical innovations in hardware and software will enable remarkable progress in reducing costs and in increasing the sensitivity of detection of sequences transcribed at low levels, the accuracy of quantification and detection of splice forms, and the prediction of the whole structures of transcripts. sequence-based transcriptome analysis has recently been applied to various organisms: arabidopsis thaliana, yeasts, drosophila melanogaster, and human. during this study, two types of rice transcriptome analysis were reported, focusing on the transcriptional differences in two rice subspecies and their reciprocal hybrids and in eight organs from different developmental stages of oryza sativa l. ssp. indica'93-11'. we analyzed salinity stress-inducible transcripts and constructed gene models based on the pilling up of short reads by using the cufflinks program. this approach should help to discover novel gene models without reliance on gene annotation.  microarray-based gene expression profiling is limited to the analysis of annotated genes. in our mrna-seq analysis, unannotated salinity stress-inducible transcripts were identified on the basis of the piling up of mapped reads without reliance on gene annotation or fl-cdna sequences. some of these novel transcripts had orfs encoding putative functional proteins and were differentially expressed in response to salinity stress. mrna-seq was valid as a gene expression profiling technology for quantifying the abundance of previously annotated genes. our findings will contribute to improvement of our rap-db and to further sequence-based gene expression profiling in various organisms. 
microarray technology is limited to monitoring the expression of previously annotated genes that have corresponding probes on the array. of these unannotated transcripts, and had orfs similar to those encoding the amino acid sequences of functional proteins in a blastx search against uniprot and refseq databases. thirty-six-base-pair reads from whole mrnas were mapped to the rice genomic sequence: % to % were mapped uniquely to the genome, and % to % bridged exons. computationally annotated genes have not fully been validated, because ests and full-length cdnas cannot cover entire transcribed regions. some of these unannotated transcripts encoding putative functional proteins were expressed differentially in response to salinity stress. here, mrna-seq was used to monitor whole mrnas of salinity stress-treated rice tissues. among the unannotated genes, and were differentially expressed in response to salinity stress.
 microarray technology is limited to monitoring the expression of previously annotated genes that have corresponding probes on the array. computationally annotated genes have not fully been validated, because ests and full-length cdnas cannot cover entire transcribed regions. here, mrna-seq was used to monitor whole mrnas of salinity stress-treated rice tissues.  thirty-six-base-pair reads from whole mrnas were mapped to the rice genomic sequence: % to % were mapped uniquely to the genome, and % to % bridged exons. from the piling up of short reads mapped on the genome, a series of programs comprehensively predicted and transcripts, including and currently unannotated in the rice annotation project database. of these unannotated transcripts, and had orfs similar to those encoding the amino acid sequences of functional proteins in a blastx search against uniprot and refseq databases. among the unannotated genes, and were differentially expressed in response to salinity stress. sequence-based and array-based measurements of the expression ratios of previously annotated genes were highly correlated.  unannotated transcripts were identified on the basis of the piling up of mapped reads derived from mrnas in rice. some of these unannotated transcripts encoding putative functional proteins were expressed differentially in response to salinity stress. 