BACKGROUND
as a key component of plant cell molecules, phosphorus  is an essential macronutrient for plant growth.
large quantities are used in fertiliser, but worldwide p resources will be exhausted by the end of this century  <cit> .
phosphate  starvation can generally be observed throughout an afflicted field.
visual symptoms of pi starvation  are the development of dark-green leaf colour and a reduction in shoot elongation and leaf size.
as −p progresses in wheat , the oldest leaves become chlorotic and show signs of desiccation  <cit> .
wheat is a major staple food crop in many parts of the world in terms of both cultivation area and prevalence as a food source.
to meet the increasing global demand for wheat, this crop’s exploitation of nutrients must be made more efficient and its requirement for nutritional fertilisers reduced.
because wheat is primarily grown on substrates with low p levels, such as the acidic soils of tropical and subtropical regions and the calcareous soils of temperate regions, an important constraint to wheat production is its lack of tolerance to −p.
various genetic approaches have been used to understand genetic control of −p tolerance in wheat; these include aneuploid analyses of the nulli-tetrasomic series and wheat alien chromosome addition lines of the cultivar ‘chinese spring’  and quantitative trait locus  mapping  <cit> .
qtl analyses using –p-sensitive cs and the tolerant variety ‘lovrin 10’ indicated that cs possesses positive alleles of the major qtls for p use efficiency on chromosomes 3b, 4b, and 5a  <cit> .
in another study, seven and six qtls were repeatedly detected controlling p uptake and use efficiency  <cit> .
a large number of qtls for agronomic trait changes under low or high p concentrations have been detected on all chromosomes in the hexaploid wheat genome, implying that −p tolerance is controlled by polygenes  <cit> .
however, the studies are few in number; a reverse genetic approach could help characterise genes that potentially contribute to complex multilocus traits and their global transcriptional networks in pi-starved wheat.
several technologies, including massively parallel sequencing and microarray analysis, have recently been used to simultaneously catalogue the effects of −p on the expressions of thousands of genes in model species  <cit> .
transcriptome sequencing using next generation sequencing  technology provides high-resolution data and is a powerful tool for studying global transcriptional networks.
the evaluation of sequence-based expression profiles can identify responsive genes and provide functional annotation for genes underlying complex and multilocus traits under −p in wheat.
in model species, transcriptome profiling and the quantification of gene expression levels are generally performed by mapping reads from the ngs analysis to a reference genome sequence and annotating genes.
the strategies for model species are not feasible in wheat, as its reference genome sequence and gene annotation are still incomplete; an international project to achieve these goals is currently making progress .
this project may take considerable time, because of the difficulties involved in sequencing the huge , highly-repetitive hexaploid genome of wheat.
the database of putative full-length cdnas for wheat, trifldb, has released approximately  <dig>  full-length cdnas   <cit> .
although this dataset is a useful reference for transcript mapping, it is incomplete, because  <dig>  to  <dig>  genes have been estimated per diploid genome based on the 3b chromosome of hexaploid wheat  <cit> .
recently, de novo transcript assembly analysis has made possible comprehensive analyses of transcriptomes, and several studies have detailed the transcriptome sequencing of various non-model species, including wheat, using massively parallel sequencing technology  <cit> .
de novo assembly of short sequences of transcripts enables researchers to reconstruct the sequences of entire transcriptomes, identify and catalogue all of the expressed genes, separate isoforms, and capture transcript expression levels.
although computer-based de novo assembly tools   <cit>  have been developed in conjunction with massively parallel sequencing, their usefulness in transcriptome assembly is not yet well demonstrated, and improvements can still be achieved using recent advances in bioinformatics.
some studies have used short-read sequence data obtained with the illumina sequencer for de novo assembly; others have used the relatively long-read sequence data obtained with the roche  <dig> pyrosequencing system or have adopted a hybrid approach of both short and long reads.
in addition, contig construction is greatly affected by sequence read quality  and quantity.
furthermore, the cdna library construction methods, sequencing technologies, and data pre-treatment techniques chosen influence the quality of the assembled transcriptomes  <cit> .
consequently, a comparison of several assembly programs is needed to determine the best combination of parameters, which can then serve as a guideline for sequence assembly performance.
in this study, we verified the de novo assembly approach by comparing analyses from several programs using short-read sequence data obtained from wheat cultivar cs seedlings under –p.
we constructed a wheat transcript dataset for de novo assembly and quantified gene expression.
as a reference in the gene expression analysis, we used a non-redundant set of transcripts generated from the de novo assembly and full-length cdnas.
this dataset was also used to assess transcripts, to investigate sequence similarity, for conservation analyses among several plant species, and for comparison with our previous report on rice transcript profiling under −p conditions  <cit> .
we demonstrated that an overall mechanism regulating gene expression of −p-responsive genes in wheat could be effectively characterised using short-read ngs data.
a comparison of gene expression profiles in wheat and rice revealed the presence of conserved gene expression systems, which appear to be essential to adaptation to –p conditions.
finally, we described an effective method for assembly of short transcript sequences to discover novel functional genes in the absence of a reference genome.
the transcript assembly generated in this study should serve as a useful resource for wheat genomics and genetics.
RESULTS
construction of the wheat transcript dataset
a set of transcripts consisting of cdnas from trifldb and contigs from rna-seq reads was constructed as shown in figure  <dig> .
a total of  <dig> libraries  from seedlings subjected to –p  and control  treatments were used in the rna-seq analysis.
overall,  <dig> million paired-end short-read sequences were produced by an illumina highseq  <dig> system  and used for de novo assembly after the removal of low-quality segments.
the choice of data pre-treatment technique, in which data are pre-processed to remove sequencing errors and other artefacts, influences the accuracy and precision of the gene expression analyses  <cit> .
trans-abyss, velvet-oases , and trinity were used for sequence assembly in this experiment .
each contig included various transcript isoforms, which cause redundancy.
because their inclusion could present difficulties in a quantitative analysis, the redundant contigs were removed.
the removal of these contigs generated a non-redundant set of contigs, which were assessed using various output parameters—such as number of contigs and mean contig length as a function of k-mer length—to select the appropriate assembly program.
the largest assembly, comprising  <dig>  contigs at least  <dig> bp long, was produced by trinity, while the smallest, consisting of  <dig>  contigs of at least  <dig> bp in length, was generated by oases.
however, oases yielded the best results of the three assembly programs, with a maximum contig length of  <dig>  bp, a mean contig length of  <dig>  bp, and a median contig length of  <dig> bp .
contig quantity and mean length are often used in de novo assembly assessment, but these values may not reflect contig accuracy.
we therefore assessed accuracy using the three assembly programs discussed above.
we counted the number of contigs that could be aligned with 90% or greater coverage to their corresponding full-length cdnas in trifldb .
of the three assembly programs, oases assembled the longest transcripts and had the highest percentage of alignments to the wheat full-length cdnas.
although it generated the smallest number of contigs, oases was able to reconstruct 87% of the trifldb cdnas with identity ≥90% and coverage ≥50%, or 75% of the trifldb cdnas with identity ≥90% and coverage ≥90% .
conversely, trinity produced the largest number of contigs, but could only reconstruct 16% of the trifldb cdnas .
trans-abyss could reconstruct just 6% of the trifldb cdnas .
thus, the reconstruction percentages of the latter two programs were considerably worse than those of oases, and we concluded that oases was the most appropriate program for our analyses.
of the  <dig>  contigs generated from oases,  <dig>  could not be aligned to the cdnas using blat.
when the longest contigs per locus were selected,  <dig>  of them could not be aligned with the cdnas.
because the  <dig>  non-redundant contigs could indicate homology to either plants or non-plants, we treated them as contaminants and removed them from the analysis, except when a blastx hit matched a land plant.
the contigs were analysed for similarity/sequence conservation against nr  using a blastx search .
a total of  <dig>  contigs were found in land plants , with  <dig> % in monocots, such as hordeum vulgare ssp.
vulgare , brachypodium distachyon , oryza sativa ssp.
japonica group , sorghum bicolor , o. sativa ssp.
indica group , t. aestivum , zea mays , and other h. vulgare subspecies  .
notably,  <dig>  of the  <dig>  contigs could not be aligned to full-length cdnas in the barley full-length cdna database  using blat, suggesting that most of the contigs were distinct from the barley transcripts.
finally, we obtained  <dig>  wheat transcripts by combining the  <dig>  non-redundant contigs with  <dig>  non-redundant full-length cdnas .
this dataset is available upon request.
mapping reads to the wheat transcript dataset
characterisations of contigs and cdnas from the  <dig>  transcripts are presented in table  <dig>  to determine whether our contigs reflected the effect of −p in wheat,  <dig> , <dig> of  <dig> , <dig> quality-controlled paired-end reads were aligned back to the wheat transcript dataset  using bowtie , which is designed for short-read mapping onto the genome.
we used  <dig> – <dig> % of the uniquely-mapped reads in the expression analysis of each treatment.
replicates in all treatments were highly correlated , and reads from the same treatment were merged to calculate reads per kilobase of exon model per million mapped reads , which indicates the relative transcription amount.
we observed that 70%  of the  <dig>  non-redundant cdnas were associated with at least one unique read.
on average, 92% of transcript length was covered by reads, with coverage depth of each transcript averaging  <dig> reads.
these statistics suggested that nearly complete coverage of the entire cdna length was obtained in this analysis.
many reads were impossible to align back to the contigs, possibly because the redundant contigs had been removed in the previous analysis .
trifldb: http://trifldb.psc.riken.jp/index.pl.
distribution of rpkm values was assessed by comparing contigs and cdnas .
average rpkm values from each treatment ranged from  <dig> – <dig>  for contigs and  <dig> – <dig>  for cdnas.
the distributions of maximum and median rpkm values were lower for contigs than cdnas, suggesting that most cdnas were highly expressed, probably reflecting the overall distribution of gene expression in this study.
we identified  <dig> contigs, however, with higher expressions than the average observed for cdnas under one or more of the four tissue/treatment combinations.
some of these contigs may correspond to transcripts rarely represented in the full-length cdna library or transcripts that were expressed at low levels in other conditions.
identification of differentially-expressed wheat transcripts under −p
a g-test  of the rpkm-derived read counts was performed to detect differences in gene expression between control and –p-treated plants, and to identify responsive wheat transcripts under –p conditions .
in roots,  <dig>  transcripts were upregulated and  <dig> were downregulated; in shoots,  <dig>  and  <dig>  transcripts were upregulated and downregulated, respectively.
the most numerous responsive transcripts were those upregulated in shoots, which were more than twice as abundant as those upregulated in roots .
this result may be a consequence of the fact that most of the cdnas in trifldb were constructed from shoot samples.
on average,  <dig> % of the transcripts in the responsive transcript set were contigs not covered by cdnas in trifldb .
the wheat transcript set contained more contigs than cdnas, but there were more responsive cdnas than contigs.
the low-expression contigs might reflect small differences in expression between the two conditions that may not be detectable statistically, making the characterization of these transcripts difficult .
we were able to demonstrate, nonetheless, that the de novo assembly strategy can improve transcriptome analysis of a non-model species and that these upregulated or downregulated contigs could be functionally annotated as pi responses in wheat under –p.
functional annotation of responsive transcripts under –p
to obtain a functional annotation of upregulated wheat transcripts under –p conditions, we used gene ontology  biological process categories.
go annotations were assigned to  <dig> %  of the root transcripts and  <dig> %  of the shoot transcripts.
we identified the top  <dig> go categories into which upregulated transcripts in roots and shoots were distributed .
although the numbers of transcripts differed, the overall categorisations were very similar between roots and shoots.
sixteen of the  <dig> go categories were represented in both roots and shoots: oxidation-reduction process, protein phosphorylation, transmembrane transport, metabolic process, regulation of transcription , lipid metabolic process, proteolysis, transport, translation and carbohydrate metabolic process, defence response, response to oxidative stress, biosynthetic process, cation transport, response to stress, and cell redox homeostasis.
transcripts from both roots and shoots most belonged to the oxidation-reduction process category.
this category was represented by many instances of cytochrome p <dig> and oxidoreductase 2og-feii, both of which were reported to be upregulated in response to oxidative stress in arabidopsis <cit> .
the protein phosphorylation category included several protein kinases, which may participate in signal transduction during many cellular processes under –p conditions, including metabolism, transcription, and cell cycle progression and differentiation.
the lipid metabolic process category had many instances of lipase, which can alter membrane lipid composition to maintain internal pi homeostasis under −p  <cit> .
analysis of ips <dig> expression by rna-seq and qrt-pcr
because rpkm values obtained from rna-seq indicated significant upregulation of taips <dig>  the wheat ips <dig> homologue, in root and shoot samples after  <dig> d of –p conditions, its expression levels were further analysed.
ips <dig>  a non-protein coding gene that includes highly conserved sequences near mir <dig> complementary regions in different plant species  <cit> , has been found to be strongly upregulated under −p in rice  <cit> .
the homologous sequences taips <dig>  , taips <dig>  , and taips <dig>   have been isolated and deposited in genbank, but molecular studies, such as expression analyses, have not been reported.
based on rna-seq, expression levels increased 341-fold in roots and 13-fold in shoots under –p conditions.
in this study, we did not distinguish homeologous forms, and we present the results for each gene as an integration of counts from the three genomes of hexaploid wheat as reported by pellny et al.
<cit> .
the primary goal of our study to identify –p-responsive wheat genes, but we noted that using an integration of counts from the three homeologous forms for each gene could cause problems with transcript abundance when reads were aligned back to the wheat transcript set using bowtie.
a quantitative real-time pcr  analysis was used to confirm these results under identical conditions.
expression levels increased 368-fold in roots and 17-fold in shoots according to qrt-pcr, confirming the results obtained using rna-seq analysis .
these data demonstrated that these transcriptomes accurately reflected the response of wheat to pi starvation and suggest that the ips1-mediated signalling cascade may also function in wheat.
a highly conserved phr1-ips1-mir399-ubc24/pho <dig> signalling cascade
members of the ips <dig> gene family from different plant species, including rice, have a highly-conserved 23-nt-long motif that exhibits complementarity with the mir <dig> involved in pi response  <cit> .
because ips <dig> is strongly upregulated in rice  <cit> , we investigated the conserved responsive genes by comparing the wheat transcripts with rice transcripts used in our previous rna-seq study  <cit> .
to assign rice annotation project  ids to responsive transcripts in wheat under –p conditions, we searched for homology to rap  proteins using blastp.
of the  <dig>  wheat transcripts in our constructed dataset, open reading frames  were predicted for  <dig> .
among these,  <dig>  transcripts with reciprocal top hits  were identified.
a total of  <dig>  transcripts were found to be orthologous to their corresponding rice genes based on a threshold of ≥ 60% identity and ≥ 60% coverage.
in addition, blastn was used to identify rice orthologues of wheat non-protein coding genes.
orthologues of the responsive wheat transcripts were verified for their responses by comparing their data with those of the responsive rice transcripts.
because most well-characterised transcripts have been found to be upregulated under −p conditions in plants, we focused on the upregulated transcripts .
in wheat,  <dig> of the  <dig>  upregulated transcripts in roots and  <dig>  of the  <dig>  upregulated transcripts in shoots possessed orthologues, which were then compared with those found in the responsive rice transcripts.
thirty-nine transcripts in roots and  <dig> transcripts in shoots were upregulated in both wheat and rice, suggesting their importance in adaptation to –p conditions in poaceae species .
the  <dig> transcripts upregulated in both roots and shoots are of particular interest.
the pi regulatory mechanism has been elucidated through arabidopsis mutant analyses and involves phr <dig> , a myb-type transcription factor, acting as a key factor in regulating downstream –p-responsive gene expression, including that of ips <dig>  through the p1bs cis-acting element  of the phr <dig> binding site .
this mechanism is conserved in vascular plants and unicellular algae  <cit> .
as in conserved genes in wheat and rice, the p1bs cis-element exists in the promoter of wheat-orthologous rice genes; ips <dig> is upregulated primarily in the roots of wheat  and rice  <cit> , suggesting conservation of the ips1-mediated signalling cascade under the control of phr <dig> .
in arabidopsis, siz <dig>   <cit> , mir <dig>  and ubc24/pho <dig>  are also involved in the cascade   <cit> .
according to genome sequencing information, siz <dig> is conserved in rice  <cit> ; mir <dig> expression has been confirmed in rice, and a potential orthologue of ubc24/pho <dig> has been identified by assembly using genomic dna and expressed sequence tags   <cit> .
transcribed sequences encoding mir <dig>  mir399-binding sites, or protein sequences homologous to n- or c-terminal extensions of ubc24/pho <dig> have been observed in wheat  <cit> .
we also detected rice-orthologous wheat transcripts of phr <dig>  ubc24/pho <dig>  and siz <dig>  although these transcripts were not highly responsive to −p in wheat  or rice  <cit> .
mature mirnas shorter than the sequencing length of the short reads could be difficult to detect in our rna-seq libraries, however.
based on the identification of hvmir399s  <cit>  and the upregulation of hvips <dig> <cit> , the ips1-mir399-ubc24/pho <dig> signalling cascade appears to be conserved in barley  <cit> .
consequently, the phr1-ips1-mir399-ubc24/pho <dig> signalling cascade is probably also conserved in wheat .
in silico comparative analysis  <cit>  suggested the presence of conserved cis-acting elements, such as the p1bs-like motif, and trans-acting factors that are capable of regulating the sole putative wheat high-affinity phosphate transporter tapt <dig>  which is expressed in a tissue-specific and pi-dependent fashion in both monocots and dicots.
we also studied the upregulated expressions of other responsive rice-orthologous wheat transcripts .
spx <dig>  upregulated under –p conditions, acts as a negative regulator for the expression of phr <dig> and, consequently, pi transporters in rice  <cit> .
rns1- <dig>   <cit> , pap   <cit> , and pho <dig>  may function in the remobilisation of pi.
in arabidopsis, phr <dig> sumoylated by siz <dig> positively regulates rns <dig> expression  <cit> .
pht1; <dig>   <cit>  may function in the uptake of pi.
atpht1; <dig> repression and atpht1; <dig> expression are attenuated in the pho <dig> mutant of arabidopsis <cit> .
genes involved in lipid metabolism, such as those that function in the synthesis of galactolipids and sulfolipids, were strongly induced under −p conditions.
these include sqd <dig>   <cit> , mgd   <cit> , and gdpd   <cit> , which may work to re-route pi during lipid metabolism.
udp-glucose pyrophosphorylase may act as a sugar-signalling network under –p conditions  <cit> .
these factors may be related to the phr1-ips1-mir399-ubc24/pho <dig> signalling cascade, as pi remobilisation, pi uptake, and pi-related metabolism are under the control of that cascade .
upregulation of several transcripts in this cascade was also confirmed by qrt-pcr .
these results suggested that the phr1-ips1-mir399-ubc24/pho <dig> signalling cascade functions as a potent adaptation to −p in wheat.
transcriptional control of other −p responsive transcripts
we identified several upregulated regulatory wheat transcripts involved in −p response from top hits to arabidopsis proteins using blastp.
wrky <dig>  and pho <dig>  transcripts were upregulated in roots.
in arabidopsis, enhanced wrky <dig> mutants are more sensitive to −p conditions, having lower pi contents in shoots compared with wild-type seedlings  <cit> .
wrky <dig> is involved in responses to –p by regulating phosphate <dig>  expression  <cit> .
pho <dig> plays a role in pi translation from roots to shoots, aiding plant adaptation to –p  <cit> .
salt and mannitol stress-inducible wrky <dig>  and skip  <cit> , and wrky <dig>   <cit> , which have pivotal roles in determining the balance between salicylic acid- and jasmonic acid-dependent defence pathways, were also upregulated in shoots.
swi3c , swi3d , and gcn5-related n-acetyltransferase  were upregulated in shoots.
swi/snf  chromatin-remodelling complexes mediate atp-dependent alterations of dna-histone contacts.
histone h2a.z regulates the expression of several classes of phosphate starvation response genes  <cit> .
although mechanisms of transcription rate modulation entailing chromatin structure alteration have not been fully elucidated under –p in arabidopsis, swi/snf complexes and gcn <dig> histone acetyltransferase are necessary for full induction of several phosphatase genes and pho <dig> high-affinity phosphate transporter gene in yeast  <cit> .
finally, we investigated downregulated wheat transcripts under –p and found that psbq-like  <dig> , psbq-like  <dig> , psbp- <dig> , psbq- <dig> , lhcb <dig>  , lhca <dig> , lhcb <dig> , and photosystem i reaction centre subunit n  were downregulated in shoots.
these genes are involved in photosynthesis.
in brassica nigra leaf petiole suspension cells, the rate of photosynthesis and photosynthetic product partitioning were altered under –p  <cit> .
glyceraldehyde-3-phosphate dehydrogenase  and sedoheptulose- <dig> -bisphosphatase   were also downregulated in shoots under –p.
these are key regulatory enzymes in co <dig> reduction and the regeneration phase of the calvin cycle for carbon fixation pathways, respectively.
sbpase, which catalyses the dephosphorylation of sedoheptulose- <dig> -bisphosphate into sedoheptulose-7-phosphate and pi, is specific to the eukaryotic calvin cycle and plays vital roles in regulating the pathways in the cycle  <cit>  and improving photosynthetic capacity  <cit> .
downregulation under –p of transcripts associated with photosynthesis and alteration of the balance of carbon metabolism are well documented in arabidopsis <cit> .
CONCLUSIONS
in this study, we demonstrated the use of short-read sequence data to rapidly characterise a wheat transcriptome and have contributed significantly to the corpus of wheat transcript data.
differentially-expressed transcripts under –p in wheat included, on average,  <dig> % of the contigs not covered by cdnas in trifldb.
the induction of these responses requires a sophisticated regulatory system, however, and details of this regulation in wheat have only recently begun to be elucidated.
comparison of upregulated wheat and rice transcripts revealed that the phr1-ips1-mir399-ubc24/pho <dig> signalling cascade is conserved in both crop species.
data from our previous study of rice transcripts and other supportive studies in other species confirmed that our analysis captured the transcriptome for –p response in wheat.
this study thus represents a genomic approach to discovering wheat transcripts when genome sequences are unavailable.
this contribution is significant to the development of genomic resources for wheat and other species and should accelerate the progress of functional genomic studies and breeding programs.

 as a key component of plant cell molecules, phosphorus is an essential macronutrient for plant growth. large quantities are used in fertiliser, but worldwide p resources will be exhausted by the end of this century. phosphate starvation can generally be observed throughout an afflicted field. visual symptoms of pi starvation are the development of dark-green leaf colour and a reduction in shoot elongation and leaf size. as −p progresses in wheat, the oldest leaves become chlorotic and show signs of desiccation. wheat is a major staple food crop in many parts of the world in terms of both cultivation area and prevalence as a food source. to meet the increasing global demand for wheat, this crop’s exploitation of nutrients must be made more efficient and its requirement for nutritional fertilisers reduced. because wheat is primarily grown on substrates with low p levels, such as the acidic soils of tropical and subtropical regions and the calcareous soils of temperate regions, an important constraint to wheat production is its lack of tolerance to −p. various genetic approaches have been used to understand genetic control of −p tolerance in wheat; these include aneuploid analyses of the nulli-tetrasomic series and wheat alien chromosome addition lines of the cultivar ‘chinese spring’ and quantitative trait locus mapping. qtl analyses using –p-sensitive cs and the tolerant variety ‘lovrin 10’ indicated that cs possesses positive alleles of the major qtls for p use efficiency on chromosomes 3b, 4b, and 5a. in another study, seven and six qtls were repeatedly detected controlling p uptake and use efficiency. a large number of qtls for agronomic trait changes under low or high p concentrations have been detected on all chromosomes in the hexaploid wheat genome, implying that −p tolerance is controlled by polygenes. however, the studies are few in number; a reverse genetic approach could help characterise genes that potentially contribute to complex multilocus traits and their global transcriptional networks in pi-starved wheat. several technologies, including massively parallel sequencing and microarray analysis, have recently been used to simultaneously catalogue the effects of −p on the expressions of thousands of genes in model species. transcriptome sequencing using next generation sequencing technology provides high-resolution data and is a powerful tool for studying global transcriptional networks. the evaluation of sequence-based expression profiles can identify responsive genes and provide functional annotation for genes underlying complex and multilocus traits under −p in wheat. in model species, transcriptome profiling and the quantification of gene expression levels are generally performed by mapping reads from the ngs analysis to a reference genome sequence and annotating genes. the strategies for model species are not feasible in wheat, as its reference genome sequence and gene annotation are still incomplete; an international project to achieve these goals is currently making progress. this project may take considerable time, because of the difficulties involved in sequencing the huge, highly-repetitive hexaploid genome of wheat. the database of putative full-length cdnas for wheat, trifldb, has released approximately full-length cdnas. although this dataset is a useful reference for transcript mapping, it is incomplete, because to genes have been estimated per diploid genome based on the 3b chromosome of hexaploid wheat. recently, de novo transcript assembly analysis has made possible comprehensive analyses of transcriptomes, and several studies have detailed the transcriptome sequencing of various non-model species, including wheat, using massively parallel sequencing technology. de novo assembly of short sequences of transcripts enables researchers to reconstruct the sequences of entire transcriptomes, identify and catalogue all of the expressed genes, separate isoforms, and capture transcript expression levels. although computer-based de novo assembly tools have been developed in conjunction with massively parallel sequencing, their usefulness in transcriptome assembly is not yet well demonstrated, and improvements can still be achieved using recent advances in bioinformatics. some studies have used short-read sequence data obtained with the illumina sequencer for de novo assembly; others have used the relatively long-read sequence data obtained with the roche pyrosequencing system or have adopted a hybrid approach of both short and long reads. in addition, contig construction is greatly affected by sequence read quality and quantity. furthermore, the cdna library construction methods, sequencing technologies, and data pre-treatment techniques chosen influence the quality of the assembled transcriptomes. consequently, a comparison of several assembly programs is needed to determine the best combination of parameters, which can then serve as a guideline for sequence assembly performance. in this study, we verified the de novo assembly approach by comparing analyses from several programs using short-read sequence data obtained from wheat cultivar cs seedlings under –p. we constructed a wheat transcript dataset for de novo assembly and quantified gene expression. as a reference in the gene expression analysis, we used a non-redundant set of transcripts generated from the de novo assembly and full-length cdnas. this dataset was also used to assess transcripts, to investigate sequence similarity, for conservation analyses among several plant species, and for comparison with our previous report on rice transcript profiling under −p conditions. we demonstrated that an overall mechanism regulating gene expression of −p-responsive genes in wheat could be effectively characterised using short-read ngs data. a comparison of gene expression profiles in wheat and rice revealed the presence of conserved gene expression systems, which appear to be essential to adaptation to –p conditions. finally, we described an effective method for assembly of short transcript sequences to discover novel functional genes in the absence of a reference genome. the transcript assembly generated in this study should serve as a useful resource for wheat genomics and genetics.  construction of the wheat transcript dataset a set of transcripts consisting of cdnas from trifldb and contigs from rna-seq reads was constructed as shown in figure. a total of libraries from seedlings subjected to –p and control treatments were used in the rna-seq analysis. overall, million paired-end short-read sequences were produced by an illumina highseq system and used for de novo assembly after the removal of low-quality segments. the choice of data pre-treatment technique, in which data are pre-processed to remove sequencing errors and other artefacts, influences the accuracy and precision of the gene expression analyses. trans-abyss, velvet-oases, and trinity were used for sequence assembly in this experiment. each contig included various transcript isoforms, which cause redundancy. because their inclusion could present difficulties in a quantitative analysis, the redundant contigs were removed. the removal of these contigs generated a non-redundant set of contigs, which were assessed using various output parameters—such as number of contigs and mean contig length as a function of k-mer length—to select the appropriate assembly program. the largest assembly, comprising contigs at least bp long, was produced by trinity, while the smallest, consisting of contigs of at least bp in length, was generated by oases. however, oases yielded the best results of the three assembly programs, with a maximum contig length of bp, a mean contig length of bp, and a median contig length of bp. contig quantity and mean length are often used in de novo assembly assessment, but these values may not reflect contig accuracy. we therefore assessed accuracy using the three assembly programs discussed above. we counted the number of contigs that could be aligned with 90% or greater coverage to their corresponding full-length cdnas in trifldb. of the three assembly programs, oases assembled the longest transcripts and had the highest percentage of alignments to the wheat full-length cdnas. although it generated the smallest number of contigs, oases was able to reconstruct 87% of the trifldb cdnas with identity ≥90% and coverage ≥50%, or 75% of the trifldb cdnas with identity ≥90% and coverage ≥90%. conversely, trinity produced the largest number of contigs, but could only reconstruct 16% of the trifldb cdnas. trans-abyss could reconstruct just 6% of the trifldb cdnas. thus, the reconstruction percentages of the latter two programs were considerably worse than those of oases, and we concluded that oases was the most appropriate program for our analyses. of the contigs generated from oases, could not be aligned to the cdnas using blat. when the longest contigs per locus were selected, of them could not be aligned with the cdnas. because the non-redundant contigs could indicate homology to either plants or non-plants, we treated them as contaminants and removed them from the analysis, except when a blastx hit matched a land plant. the contigs were analysed for similarity/sequence conservation against nr using a blastx search. a total of contigs were found in land plants, with % in monocots, such as hordeum vulgare ssp. vulgare, brachypodium distachyon, oryza sativa ssp. japonica group, sorghum bicolor, o. sativa ssp. indica group, t. aestivum, zea mays, and other h. vulgare subspecies. notably, of the contigs could not be aligned to full-length cdnas in the barley full-length cdna database using blat, suggesting that most of the contigs were distinct from the barley transcripts. finally, we obtained wheat transcripts by combining the non-redundant contigs with non-redundant full-length cdnas. this dataset is available upon request. mapping reads to the wheat transcript dataset characterisations of contigs and cdnas from the transcripts are presented in table to determine whether our contigs reflected the effect of −p in wheat,, of, quality-controlled paired-end reads were aligned back to the wheat transcript dataset using bowtie, which is designed for short-read mapping onto the genome. we used – % of the uniquely-mapped reads in the expression analysis of each treatment. replicates in all treatments were highly correlated, and reads from the same treatment were merged to calculate reads per kilobase of exon model per million mapped reads, which indicates the relative transcription amount. we observed that 70% of the non-redundant cdnas were associated with at least one unique read. on average, 92% of transcript length was covered by reads, with coverage depth of each transcript averaging reads. these statistics suggested that nearly complete coverage of the entire cdna length was obtained in this analysis. many reads were impossible to align back to the contigs, possibly because the redundant contigs had been removed in the previous analysis. trifldb: http://trifldb.psc.riken.jp/index.pl. distribution of rpkm values was assessed by comparing contigs and cdnas. average rpkm values from each treatment ranged from – for contigs and – for cdnas. the distributions of maximum and median rpkm values were lower for contigs than cdnas, suggesting that most cdnas were highly expressed, probably reflecting the overall distribution of gene expression in this study. we identified contigs, however, with higher expressions than the average observed for cdnas under one or more of the four tissue/treatment combinations. some of these contigs may correspond to transcripts rarely represented in the full-length cdna library or transcripts that were expressed at low levels in other conditions. identification of differentially-expressed wheat transcripts under −p a g-test of the rpkm-derived read counts was performed to detect differences in gene expression between control and –p-treated plants, and to identify responsive wheat transcripts under –p conditions. in roots, transcripts were upregulated and were downregulated; in shoots, and transcripts were upregulated and downregulated, respectively. the most numerous responsive transcripts were those upregulated in shoots, which were more than twice as abundant as those upregulated in roots. this result may be a consequence of the fact that most of the cdnas in trifldb were constructed from shoot samples. on average, % of the transcripts in the responsive transcript set were contigs not covered by cdnas in trifldb. the wheat transcript set contained more contigs than cdnas, but there were more responsive cdnas than contigs. the low-expression contigs might reflect small differences in expression between the two conditions that may not be detectable statistically, making the characterization of these transcripts difficult. we were able to demonstrate, nonetheless, that the de novo assembly strategy can improve transcriptome analysis of a non-model species and that these upregulated or downregulated contigs could be functionally annotated as pi responses in wheat under –p. functional annotation of responsive transcripts under –p to obtain a functional annotation of upregulated wheat transcripts under –p conditions, we used gene ontology biological process categories. go annotations were assigned to % of the root transcripts and % of the shoot transcripts. we identified the top go categories into which upregulated transcripts in roots and shoots were distributed. although the numbers of transcripts differed, the overall categorisations were very similar between roots and shoots. sixteen of the go categories were represented in both roots and shoots: oxidation-reduction process, protein phosphorylation, transmembrane transport, metabolic process, regulation of transcription, lipid metabolic process, proteolysis, transport, translation and carbohydrate metabolic process, defence response, response to oxidative stress, biosynthetic process, cation transport, response to stress, and cell redox homeostasis. transcripts from both roots and shoots most belonged to the oxidation-reduction process category. this category was represented by many instances of cytochrome p and oxidoreductase 2og-feii, both of which were reported to be upregulated in response to oxidative stress in arabidopsis. the protein phosphorylation category included several protein kinases, which may participate in signal transduction during many cellular processes under –p conditions, including metabolism, transcription, and cell cycle progression and differentiation. the lipid metabolic process category had many instances of lipase, which can alter membrane lipid composition to maintain internal pi homeostasis under −p. analysis of ips expression by rna-seq and qrt-pcr because rpkm values obtained from rna-seq indicated significant upregulation of taips the wheat ips homologue, in root and shoot samples after d of –p conditions, its expression levels were further analysed. ips a non-protein coding gene that includes highly conserved sequences near mir complementary regions in different plant species, has been found to be strongly upregulated under −p in rice. the homologous sequences taips, taips, and taips have been isolated and deposited in genbank, but molecular studies, such as expression analyses, have not been reported. based on rna-seq, expression levels increased 341-fold in roots and 13-fold in shoots under –p conditions. in this study, we did not distinguish homeologous forms, and we present the results for each gene as an integration of counts from the three genomes of hexaploid wheat as reported by pellny et al.. the primary goal of our study to identify –p-responsive wheat genes, but we noted that using an integration of counts from the three homeologous forms for each gene could cause problems with transcript abundance when reads were aligned back to the wheat transcript set using bowtie. a quantitative real-time pcr analysis was used to confirm these results under identical conditions. expression levels increased 368-fold in roots and 17-fold in shoots according to qrt-pcr, confirming the results obtained using rna-seq analysis. these data demonstrated that these transcriptomes accurately reflected the response of wheat to pi starvation and suggest that the ips1-mediated signalling cascade may also function in wheat. a highly conserved phr1-ips1-mir399-ubc24/pho signalling cascade members of the ips gene family from different plant species, including rice, have a highly-conserved 23-nt-long motif that exhibits complementarity with the mir involved in pi response. because ips is strongly upregulated in rice, we investigated the conserved responsive genes by comparing the wheat transcripts with rice transcripts used in our previous rna-seq study. to assign rice annotation project ids to responsive transcripts in wheat under –p conditions, we searched for homology to rap proteins using blastp. of the wheat transcripts in our constructed dataset, open reading frames were predicted for. among these, transcripts with reciprocal top hits were identified. a total of transcripts were found to be orthologous to their corresponding rice genes based on a threshold of ≥ 60% identity and ≥ 60% coverage. in addition, blastn was used to identify rice orthologues of wheat non-protein coding genes. orthologues of the responsive wheat transcripts were verified for their responses by comparing their data with those of the responsive rice transcripts. because most well-characterised transcripts have been found to be upregulated under −p conditions in plants, we focused on the upregulated transcripts. in wheat, of the upregulated transcripts in roots and of the upregulated transcripts in shoots possessed orthologues, which were then compared with those found in the responsive rice transcripts. thirty-nine transcripts in roots and transcripts in shoots were upregulated in both wheat and rice, suggesting their importance in adaptation to –p conditions in poaceae species. the transcripts upregulated in both roots and shoots are of particular interest. the pi regulatory mechanism has been elucidated through arabidopsis mutant analyses and involves phr, a myb-type transcription factor, acting as a key factor in regulating downstream –p-responsive gene expression, including that of ips through the p1bs cis-acting element of the phr binding site. this mechanism is conserved in vascular plants and unicellular algae. as in conserved genes in wheat and rice, the p1bs cis-element exists in the promoter of wheat-orthologous rice genes; ips is upregulated primarily in the roots of wheat and rice, suggesting conservation of the ips1-mediated signalling cascade under the control of phr. in arabidopsis, siz, mir and ubc24/pho are also involved in the cascade. according to genome sequencing information, siz is conserved in rice; mir expression has been confirmed in rice, and a potential orthologue of ubc24/pho has been identified by assembly using genomic dna and expressed sequence tags. transcribed sequences encoding mir mir399-binding sites, or protein sequences homologous to n- or c-terminal extensions of ubc24/pho have been observed in wheat. we also detected rice-orthologous wheat transcripts of phr ubc24/pho and siz although these transcripts were not highly responsive to −p in wheat or rice. mature mirnas shorter than the sequencing length of the short reads could be difficult to detect in our rna-seq libraries, however. based on the identification of hvmir399s and the upregulation of hvips, the ips1-mir399-ubc24/pho signalling cascade appears to be conserved in barley. consequently, the phr1-ips1-mir399-ubc24/pho signalling cascade is probably also conserved in wheat. in silico comparative analysis suggested the presence of conserved cis-acting elements, such as the p1bs-like motif, and trans-acting factors that are capable of regulating the sole putative wheat high-affinity phosphate transporter tapt which is expressed in a tissue-specific and pi-dependent fashion in both monocots and dicots. we also studied the upregulated expressions of other responsive rice-orthologous wheat transcripts. spx upregulated under –p conditions, acts as a negative regulator for the expression of phr and, consequently, pi transporters in rice. rns1-, pap, and pho may function in the remobilisation of pi. in arabidopsis, phr sumoylated by siz positively regulates rns expression. pht1; may function in the uptake of pi. atpht1; repression and atpht1; expression are attenuated in the pho mutant of arabidopsis. genes involved in lipid metabolism, such as those that function in the synthesis of galactolipids and sulfolipids, were strongly induced under −p conditions. these include sqd, mgd, and gdpd, which may work to re-route pi during lipid metabolism. udp-glucose pyrophosphorylase may act as a sugar-signalling network under –p conditions. these factors may be related to the phr1-ips1-mir399-ubc24/pho signalling cascade, as pi remobilisation, pi uptake, and pi-related metabolism are under the control of that cascade. upregulation of several transcripts in this cascade was also confirmed by qrt-pcr. these results suggested that the phr1-ips1-mir399-ubc24/pho signalling cascade functions as a potent adaptation to −p in wheat. transcriptional control of other −p responsive transcripts we identified several upregulated regulatory wheat transcripts involved in −p response from top hits to arabidopsis proteins using blastp. wrky and pho transcripts were upregulated in roots. in arabidopsis, enhanced wrky mutants are more sensitive to −p conditions, having lower pi contents in shoots compared with wild-type seedlings. wrky is involved in responses to –p by regulating phosphate expression. pho plays a role in pi translation from roots to shoots, aiding plant adaptation to –p. salt and mannitol stress-inducible wrky and skip, and wrky, which have pivotal roles in determining the balance between salicylic acid- and jasmonic acid-dependent defence pathways, were also upregulated in shoots. swi3c, swi3d, and gcn5-related n-acetyltransferase were upregulated in shoots. swi/snf chromatin-remodelling complexes mediate atp-dependent alterations of dna-histone contacts. histone h2a.z regulates the expression of several classes of phosphate starvation response genes. although mechanisms of transcription rate modulation entailing chromatin structure alteration have not been fully elucidated under –p in arabidopsis, swi/snf complexes and gcn histone acetyltransferase are necessary for full induction of several phosphatase genes and pho high-affinity phosphate transporter gene in yeast. finally, we investigated downregulated wheat transcripts under –p and found that psbq-like, psbq-like, psbp-, psbq-, lhcb, lhca, lhcb, and photosystem i reaction centre subunit n were downregulated in shoots. these genes are involved in photosynthesis. in brassica nigra leaf petiole suspension cells, the rate of photosynthesis and photosynthetic product partitioning were altered under –p. glyceraldehyde-3-phosphate dehydrogenase and sedoheptulose- -bisphosphatase were also downregulated in shoots under –p. these are key regulatory enzymes in co reduction and the regeneration phase of the calvin cycle for carbon fixation pathways, respectively. sbpase, which catalyses the dephosphorylation of sedoheptulose- -bisphosphate into sedoheptulose-7-phosphate and pi, is specific to the eukaryotic calvin cycle and plays vital roles in regulating the pathways in the cycle and improving photosynthetic capacity. downregulation under –p of transcripts associated with photosynthesis and alteration of the balance of carbon metabolism are well documented in arabidopsis.  in this study, we demonstrated the use of short-read sequence data to rapidly characterise a wheat transcriptome and have contributed significantly to the corpus of wheat transcript data. differentially-expressed transcripts under –p in wheat included, on average, % of the contigs not covered by cdnas in trifldb. the induction of these responses requires a sophisticated regulatory system, however, and details of this regulation in wheat have only recently begun to be elucidated. comparison of upregulated wheat and rice transcripts revealed that the phr1-ips1-mir399-ubc24/pho signalling cascade is conserved in both crop species. data from our previous study of rice transcripts and other supportive studies in other species confirmed that our analysis captured the transcriptome for –p response in wheat. this study thus represents a genomic approach to discovering wheat transcripts when genome sequences are unavailable. this contribution is significant to the development of genomic resources for wheat and other species and should accelerate the progress of functional genomic studies and breeding programs. 
the average expression level of the contigs was lower than that of the known cdnas, implying that these contigs included transcripts that were rarely represented in the full-length cdna library. a comparative analysis of the wheat- and rice-responsive transcripts for orthologous genes under pi-starvation revealed commonly upregulated transcripts, most of which appeared to be involved in a general response to pi starvation, namely, an ips1-mediated signalling cascade and its downstream functions such as pi remobilisation, pi uptake, and changes in pi metabolism. phosphorus is an essential macronutrient for plant growth and development. we carried out de novo assembly on million high-quality reads generated from rna-seq libraries. to elucidate these mechanisms, we used rna-seq methods to generate transcriptome profiles of the wheat variety ‘chinese spring’ responding to days of pi starvation. within the non-redundant transcript set, we identified 892– responsive transcripts in roots and shoots, corresponding on average to % of the contigs not covered by cdnas in trifldb under pi starvation. our transcriptome profiles demonstrated the impact of pi starvation on global gene expression in wheat. this study revealed that enhancement of the pi-mediated signalling cascade using ips is a potent adaptation mechanism to pi starvation that is conserved in both wheat and rice and validated the effectiveness of using short-read next-generation sequencing data for wheat transcriptome analysis in the absence of reference genome information. when compared with barley full-length cdnas, of the contigs were unalignable, suggesting that many might be distinct from barley transcripts.
 phosphorus is an essential macronutrient for plant growth and development. to modulate their p homeostasis, plants must balance p uptake, mobilisation, and partitioning to various organs. despite the worldwide importance of wheat as a cultivated food crop, molecular mechanisms associated with phosphate starvation in wheat remain unclear. to elucidate these mechanisms, we used rna-seq methods to generate transcriptome profiles of the wheat variety ‘chinese spring’ responding to days of pi starvation.  we carried out de novo assembly on million high-quality reads generated from rna-seq libraries. we then constructed a transcript dataset containing non-redundant wheat transcripts, comprising contigs and non-redundant full-length cdnas from the trifldb database. when compared with barley full-length cdnas, of the contigs were unalignable, suggesting that many might be distinct from barley transcripts. the average expression level of the contigs was lower than that of the known cdnas, implying that these contigs included transcripts that were rarely represented in the full-length cdna library. within the non-redundant transcript set, we identified 892– responsive transcripts in roots and shoots, corresponding on average to % of the contigs not covered by cdnas in trifldb under pi starvation. the relative expression level of the wheat ips homologue, taips was 341-fold higher in roots and 13-fold higher in shoots; this finding was further confirmed by qrt-pcr analysis. a comparative analysis of the wheat- and rice-responsive transcripts for orthologous genes under pi-starvation revealed commonly upregulated transcripts, most of which appeared to be involved in a general response to pi starvation, namely, an ips1-mediated signalling cascade and its downstream functions such as pi remobilisation, pi uptake, and changes in pi metabolism.  our transcriptome profiles demonstrated the impact of pi starvation on global gene expression in wheat. this study revealed that enhancement of the pi-mediated signalling cascade using ips is a potent adaptation mechanism to pi starvation that is conserved in both wheat and rice and validated the effectiveness of using short-read next-generation sequencing data for wheat transcriptome analysis in the absence of reference genome information. de novo assemblyrna-seqtranscriptomewheatphosphorusphosphate starvation 