Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

mRNA-Seq Reveals a Comprehensive Transcriptome Profile of Rice under Phosphate Stress

mRNA-Seq Reveals a Comprehensive Transcriptome Profile of Rice under Phosphate Stress Rice (2011) 4:50–65 DOI 10.1007/s12284-011-9064-0 mRNA-Seq Reveals a Comprehensive Transcriptome Profile of Rice under Phosphate Stress Youko Oono & Yoshihiro Kawahara & Hiroyuki Kanamori & Hiroshi Mizuno & Harumi Yamagata & Mayu Yamamoto & Satomi Hosokawa & Hiroshi Ikawa & Ikuko Akahane & Zuofeng Zhu & Jianzhong Wu & Takeshi Itoh & Takashi Matsumoto Received: 5 July 2011 /Accepted: 17 October 2011 /Published online: 17 November 2011 The Author(s) 2011. This article is published with open access at Springerlink.com Abstract Plants have developed several morphological and detected uniquely 10,388 transcripts with no match to any physiological strategies to adapt to phosphate stress. We RAP transcript. These transcripts that showed specific analyzed the inducible transcripts associated with phosphate response to Pi stress include those without ORFs that may starvation and over-abundant phosphate supply to characterize act as non-protein coding transcripts. With an accompanying the transcriptome in rice seedlings using the mRNA-Seq browser of the transcriptome under Pi stress, a deeper strategy. Fifty-three million reads obtained from 16 libraries understanding of the structural and functional features of both under various phosphate stress and recovery treatments were annotated and unannotated Pi stress-responsive transcripts can uniquely mapped to the rice genome. Transcripts identified provide useful information in improving Pi acquisition and specifically tagged to 40,574 (root) and 39,748 (shoot) Rice utilization in rice and other cereal crops. Annotation Project (RAP) transcripts. Additionally, we . . Keywords Phosphorus Phosphate starvation Over- . . . abundant phosphate supply Rice mRNA-Seq Electronic supplementary material The online version of this article Transcriptome (doi:10.1007/s12284-011-9064-0) contains supplementary material, which is available to authorized users. : : : : Y. Oono Y. Kawahara H. Kanamori H. Mizuno Introduction : : : : : H. Yamagata M. Yamamoto S. Hosokawa Z. Zhu J. Wu T. Itoh T. Matsumoto (*) Phosphorus (P) is a component of key cellular molecules Plant Genome Research Unit, Agrogenomics Research Center, National Institute of Agrobiological Sciences (NIAS), such as nucleic acids, proteins, phospholipids, phytic acid, 2-1-2 Kannondai, and ATP in plants. Plants absorb P almost exclusively in an Tsukuba, Ibaraki 305-8602, Japan inorganic phosphate form (Pi), primarily as H PO from 2 4 e-mail: mat@nias.affrc.go.jp the soil. Among the major nutrients necessary for plant H. Ikawa growth, P is the most dilute and the least mobile in soil; Forestry and Fisheries, Institute of the Society hence it is often a limiting factor for crop yield. Pi for Techno-innovation of Agriculture, starvation during farming is alleviated by the massive 446-1 Ippaizuka, Kamiyokoba, application of fertilizers. However, continuous usage of Tsukuba, Ibaraki 305-0854, Japan phosphate fertilizers may have a negative impact on the I. Akahane environment as the rock phosphate in the world is now in National Institute for Agro-Environmental Sciences (NIAES), short supply and maybe depleted within the next century 3-1-3 Kannondai, (Vance et al. 2003). On the other hand, massive Pi Tsukuba, Ibaraki 305-8604, Japan fertilization is also known to inhibit plant growth and Present Address: productivity. A thorough understanding of plant response to Z. Zhu stress due to Pi starvation (−P) as well as an over-abundant Department of Plant Genetics and Breeding, Pi supply (++P) is therefore indispensable in improving China Agricultural University, acquisition and utilization of this nutrient. Beijing 100093, China Rice (2011) 4:50–65 51 The molecular mechanisms by which plants respond and comprehensive overview of the primary molecular events acclimate to changes in the nutritional Pi concentration are resulting from phosphate stress. complex but of great importance and could be useful in developing strategies for elucidating the gene networks involved in plant response to various kinds of abiotic stress. Results The regulation system for the—P signaling pathway in plants has been proposed using studies in Arabidopsis and Morphological changes triggered by −P and ++P stress involved the sumoylation of the MYB transcription factor PHR1 by a small ubiquitin-like modifier SIZ1 in a process After the first 10 days in −P and ++P stress treatments, the dependent on E3 ligase (Miura et al. 2005; Gojon et al. rice seedlings began to show morphological changes in 2009). Located downstream of PHR1, miR399 and IPS1 roots and shoots that may be attributed to starvation or (induced by phosphate starvation1) are specifically induced over-abundant supply of phosphate. These changes gradu- by −P stress. The miR399 reciprocally regulates the ally became prominent so that after 30 days of –P stress, expression of PHO2/UBC24 and the resulting protein shoot growth was retarded with thin and fewer leaves functions as an ubiquitin-conjugating E2 enzyme (Bari et whereas the root became bristlier in texture and more al. 2006). A loss of function of PHO2/UBC24 could lead to brownish in color as compared with the control (Fig. 1a). excessive Pi accumulation in the shoot (Fujii et al. 2005). Under ++P stress however, the shoot gradually showed IPS1, a member of the Mt4/TPS1 family, was shown to signs of wilting with some yellowing of the leaf and mimic the target miR399 and offset its function (Franco- relatively inhibited root growth. The total phosphorus Zorrilla et al. 2007). The sequences near the miR399 content in −P stress-treated plants was lower than the complementary regions in different plant species, and all control for both shoots and roots (Fig. 1b). On the other known members of the IPS1 family in Arabidopsis, are hand, in ++P stress-treated plants, the total P content was highly conserved. Of note, if Pi is supplied to plants, IPS1 1.3-fold higher in the shoots but slightly lower in the roots transcripts rapidly disappear suggesting that they are than the control. The dry weight of −P and ++P stress- component riburegulators of the Pi signaling system. treated plants decreased gradually to 38.4% and 45.3% of PHR1 also upregulates many responsive genes that are the control plants after 30 days, respectively (Fig. 1c). The expressed under −P stress, including those encoding high- shoot/root weight ratio differed in the two stress treatments affinity Pi transporters, RNases, and acid phosphatases, after 10–15 days with a decrease in the −P stress-treated by binding to an imperfect palindromic sequence (P1BS, plants and an increase in ++P stress-treated plants as GNATATNC) in the promoter regions of these genes compared with the control (Fig. 1d). Although the roots (Rubio et al. 2001). However, although Pi stress due to weight did not decrease much under −P stress, there was a starvation has been characterized in detail, not much is significant decrease in the dry weight of the roots and known on the signaling pathway involved in Pi stress shoots under ++P stress. These results indicate that rice brought about by an over-abundant supply of phosphate in seedlings respond to stress due Pi starvation or an over- the soil. abundant Pi supply with visible changes in root and shoot Recently, the mRNA-Seq strategy using next generation morphology after 10 days of treatment. sequencers has become a useful tool for analyzing genome- wide gene expression and cataloguing of all transcripts Time course expression of OsIPS1 under P stress including mRNAs, non-coding RNAs and small RNAs. With a high resolution and sensitivity, the mRNA-Seq can We next analyzed the expression of OsIPS1, a Pi starvation provide detailed information on transcriptional structure of responsive non-protein coding regulatory transcript, to genes such as the precise location of transcription bound- confirm the validity of −P stress treatment. OsIPS1 was aries to a single-base resolution, reveal rare transcripts or gradually upregulated from days 1 to 10 in both root and variants, and identify splicing isoforms of known genes shoot under −P stress, after which the expression remained (Wang et al. 2009a). More importantly, it could accurately stable for up to 30 days (Fig. 2a). After transferring the quantify gene expression levels over a broad dynamic range plants to growth conditions with normal Pi concentration and detect transcripts expressed at either very low or very (+P) to allow recovery from stress, the expression was high levels including subtle changes which could not be immediately downregulated within 1 day. Under ++P stress characterized by microarray-based approaches (Li et al. however, OsIPS1 was not upregulated in both root and 2009; He et al. 2010; Lu et al. 2010; Mizuno et al. 2010). shoot (Fig. 2b). This gene was also not upregulated under We therefore used the massive parallel sequencing technol- starvation with other essential nutrients such as nitrogen, ogy by mRNA-Seq to elucidate the rice transcriptome calcium and magnesium (data not shown). Based on these under stress due to –P and ++P in order to provide a results, we focused our gene expression profiling analysis 52 Rice (2011) 4:50–65 Fig. 1 Effect of Pi stress on rice growth. a Phenotypic changes in rice represent the mean±SE for three replicates for each treatment. c, d plants after 30 days of growth in culture medium with normal Pi Changes in dry weight and shoot/root weight ratio (in mg) of rice concentration (+P, control), Pi starvation (−P), and over-abundant Pi seedlings +P, −P, and ++P treatment conditions. The values represent supply (++P). Both the shoot and root showed growth retardation under the mean±SE for three replicates for each treatment. Changes in dry Pi stress. b Total phosphorus content in roots and shoots after weight and shoot/root weight ratio became prominent at around 30 days in +P, −P, and ++P treatment conditions. The values 10 days after Pi stress treatment. by mRNA-Seq during the early growth stages particularly ces were produced and used for mapping onto the reference at 1 (early), 5 (middle), and 10 days (late) after −P and ++P Nipponbare genome sequence (Table 1). Of the 6 to 30 stress treatments. Additionally, the expression profiles at million quality-evaluated reads (Passed Filtered reads) from recovery after 10 days of −P and ++P stress treatments were each library, 12.3% to 39.8% were mapped to single also investigated (Supplementary Fig. S1). locations (unique) in the genome whereas 23.7% to 45.7% were mapped to multiple locations in the genome. The Sequencing and short-read mapping expression level of all unique transcripts mapped onto the genome were quantified as reads per kilobase of exon A total of 16 libraries were used for mRNA-Seq analysis model per million mapped (RPKM) values (Mortazavi et al. using the Illumina Genome Analyzer IIx (Illumina Inc., San 2008). On the other hand, 19.0% to 51.1% of the sequence Diego, CA, USA). Overall, 195 million short-read sequen- reads from each library had no match in the genome. These Rice (2011) 4:50–65 53 Fig. 2 OsIPS1 expression under Pi stress treatments. The expression of 30 days after stress and at 1 or 2 days after recovery from stress (Re1 and OsIPS1 in root and shoot under Pi starvation (a) and over-abundant Pi Re2). The amplification of β-tubulin in RT-PCR was used as control. supply (b) was analyzed by RT-PCR at various time points from 0 to unmapped reads may include low-quality reads, sequencing that 95% of the transcripts in the 44 K microarray platform errors, or sequences derived from adaptors and contami- could be validated by mRNA-Seq. A total of 11,888 nating organisms (Mizuno et al. 2010). transcripts from the root and 11,098 transcripts from the shoot were not represented in the 44 K array but showed mRNA-Seq vs. microarray analysis complete match with the annotated RAP transcripts. These transcripts correspond to gene models without EST and/or We correlated the mRNA-Seq data with the gene expres- full-length cDNA support as well as those predicted based sion profiles derived from analysis using the rice 44 K on gene prediction programs. The mRNA-Seq approach oligomicroarray platform (Agilent Technologies, Palo Alto, could therefore provide a more comprehensive analysis of CA, USA). First, we re-mapped the microarray probe the transcriptome including about half of RAP transcripts sequences into the updated rice genome assembly (Build not represented in the microarray platform. We used the 5.0 pseudomolecules) to make a valid comparison of the Cufflinks program for comparative assembly and estimation microarray data and the mRNA-Seq data. As a result, we of abundance of transcripts in each sample (Trapnell et al. confirmed that 30,026 RAP transcripts (23,113 loci) were 2010). By excluding the annotated transcripts in RAP and represented as probes in the 44 K array (Fig. 3). Sequence the MSU rice gene models (http://rice.plantbiology.msu. comparison with the mRNA-Seq data showed that 28,680 edu/), we were able to detect a total of 8,590 transcripts transcripts from the root and 28,650 transcripts from the from the root and 8,193 transcripts from the shoot which shoot matched with the microarray probes. This indicates have not yet been previously annotated. Table 1 Mapping of mRNA-Seq reads obtained from each sample into the reference rice genome sequence mRNA-Seq library Total PF reads Unique % Multiple % Unmapped % Root, +P_0d (control) 13,922,030 3,738,255 26.9 4,344,200 31.2 5,839,575 41.9 Root, −P_1d 8,411,989 2,314,238 27.5 2,827,342 33.6 3,270,409 38.9 Root, −P_5d 30,403,973 8,057,787 26.5 7,213,838 23.7 15,132,348 49.8 Root, −P_10d 6,696,463 2,495,317 37.3 2,459,102 36.7 1,742,044 26.0 Root, −P_10d_Re_1d 8,528,677 3,311,738 38.8 2,857,712 33.5 2,359,227 27.7 Root, ++P,_1d 15,694,862 4,738,300 30.2 5,013,389 31.9 5,943,173 37.9 Root, ++P_5d 13,586,681 2,642,094 19.4 4,000,708 29.4 6,943,879 51.1 Root, ++P_10d 11,027,653 1,353,233 12.3 4,597,218 41.7 5,077,202 46.0 Root, ++P_10d_Re_1d 16,771,413 3,827,077 22.8 5,270,793 31.4 7,673,543 45.8 Shoot, +P_0d (control) 8,515,740 2,205,725 25.9 3,181,292 37.4 3,128,723 36.7 Shoot, −P_1d 8,787,606 2,380,054 27.1 3,319,901 37.8 3,087,651 35.1 Shoot, −P_5d 9,281,244 2,278,998 24.6 4,240,468 45.7 2,761,778 29.8 Shoot, −P_10d 8,258,753 2,943,631 35.6 3,744,753 45.3 1,570,369 19.0 Shoot, ++P_1d 12,145,854 3,090,635 25.4 4,167,569 34.3 4,887,650 40.2 Shoot, ++P_5d 14,586,873 4,126,432 28.3 6,302,429 43.2 4,158,012 28.5 Shoot, ++P_10d 9,086,827 3,612,159 39.8 3,601,872 39.6 1,872,796 20.6 54 Rice (2011) 4:50–65 Fig. 3 Transcripts identified in RAP annotations and 44 K microarray probes. Unique tran- scripts mapped to the genome were classified as follows: (1) transcripts with match in RAD- DB annotation and microarray, (2) transcripts with match in RAP-DB annotation but no match in microarray, and (3) unannotated transcripts identi- fied by the Cufflinks program. In both root and shoot samples, about 95% of transcripts were supported by the microarray probes. About 50% of tran- scripts not identified by the array matched with the RAP-DB transcripts. A total of 10,388 unique unannotated transcripts were identified by Cufflinks program. Gene expression profiling the 44 K microarray platform the expression profiles for known Pi-related genes can be was performed in root samples at 1, 5, and 10 days after –P validated by mRNA-Seq. and ++P stress treatments. The signal intensities of 28,686 transcripts represented in the microarray were plotted Characterization of Pi stress-responsive transcripts against the RPKM values of corresponding transcripts (Supplementary Fig. S2). For overall gene expression, we We used the G test (FDR, <0.01) on the RPKM-derived observed an average correlation coefficient of >0.8, read counts to determine the differences in gene expression suggesting a clear validation of the microarray-based gene at various time points during −P and ++P treatments. A expression profiling data with the mRNA-Seq data. transcript is considered responsive if the expression level after recovery from stress (−P_Re_1d and ++P_Re_1d) Transcripts corresponding to known Pi-related genes reverted to the same level before treatment or the control condition. As a result, we were able to classify the We were able to identify transcripts with match to expression patterns into 24 expression patterns based on previously reported Pi-related genes including those in- the statistical significance between treatment conditions (−P volved in sumoylation (SIZ1) and signaling (OsIPS1 and and ++P) and the expression levels in root and shoot at each OsIPS2), transcription factors (phosphate starvation re- stage of stress designated as early for 1 day, middle for sponse (PHR)1 and PHR2), phosphatases (PAP2 and 5 days, and late for 10 days after stress treatment (Table 2). PAP10), etc. (Supplementary Table S1). As expected, Under stress recovery conditions, the responsive transcripts analysis using Cufflinks also revealed the abundance and which were upregulated in roots under −P stress showed differential expression of these transcripts. The expression similar downregulation as OsIPS1 after recovery from levels of transcripts from root and shoot samples at various stress starvation. The same pattern of fast and marked Pi starvation (−P_1d, −P_5d, and −P_10d) as well as over- recovery of expression after re-addition of Pi was also abundant Pi supply (++P_1d, ++P_5d, and ++P_10d) shown in Arabidopsis with >20-fold increase in AtPht1;4 treatments based on RPKM values indicate significant and AtPAP2 signal intensity as compared with the level expression of these Pi-related genes (Supplementary during P deprivation (Morcuende et al. 2007). Similarly, Table S2). Moreover, comparison of the fold change in upregulated transcripts under ++P stress conditions were RPKM values from control and treatment conditions with also downregulated after 1-day recovery from stress due to the normalized signal intensity clearly showed a similar over-abundant Pi supply. pattern of expression profile obtained by mRNA-Seq and The highest number of responsive (upregulated and microarray-based analysis in all treatments of –P and ++P downregulated) transcripts was observed in root at late stress (Supplementary Table S3). These results indicate that stage, which was almost twice the number in root at early Rice (2011) 4:50–65 55 Table 2 Number of Pi stress-responsive transcripts identified by G test Expression patterns Number of responsive transcripts Tissue Treatment Response Time point Total RAP/array RAP not Unannotated Novel Transcripts identified (%) with P1BS by array Root −P Upregulated Early (1 d) 230 201 27 2 12.6 43 Middle (5 d) 354 286 68 0 19.2 71 Late (10 d) 503 434 68 1 13.7 170 Dowregulated Early (1 d) 148 112 35 1 24.3 21 Middle (5 d) 171 156 15 0 8.8 31 Late (10 d) 246 221 24 1 10.2 54 ++P Upregulated Early (1 d) 233 192 40 1 17.6 58 Middle (5 d) 472 399 69 4 15.5 105 Late (10 d) 265 203 57 5 23.4 54 Dowregulated Early (1 d) 209 186 21 2 11.0 39 Middle (5 d) 203 181 22 0 10.8 39 Late (10 d) 208 186 22 0 10.6 52 Shoot −P Upregulated Early (1 d) 224 199 22 3 11.2 47 Middle (5 d) 117 103 14 0 12.0 25 Late (10 d) 216 182 32 2 15.7 74 Dowregulated Early (1 d) 454 400 53 1 11.9 110 Middle (5 d) 83 71 11 1 14.5 20 Late (10 d) 205 183 22 0 10.7 39 ++P Upregulated Early (1 d) 588 498 88 2 15.3 145 Middle (5 d) 466 397 68 1 14.8 94 Late (10 d) 166 139 27 0 16.3 40 Dowregulated Early (1 d) 1,136 1,014 118 4 10.7 252 Middle (5 d) 168 152 15 1 9.5 37 Late (10 d) 82 77 5 0 6.1 15 P1BS is significantly overrepresented in the 1-kb upstream region at 19.6% background level (responsive transcripts/non-responsive transcripts= 7,164/36,588) and FDR, <0.001 in Fisher’s exact test stage of –P stress (Table 2; Fig. 4). The highest number of 13.6%) of all responsive transcripts were classified as novel responsive transcripts was found in shoot at early stage and transcripts in each group. in root at middle stage of ++P stress. Downregulation of The Pi stress-responsive genes were also analyzed based transcription in shoot at early ++P treatment was most on the expression of PHR1. This gene is known to common among the 24 expression patterns. The expression upregulate −P stress-related genes such as high-affinity Pi patterns of representative transcripts such as WRKY transporters, RNases, and acid phosphatases by binding to (Os03t0321700), Pht1;2 (Os03t0150800), and OsIPS1 an imperfect palindromic sequence (P1BS, GNATATNC) in (Os03t0146800) clearly show upregulation at early, middle, the promoter regions (Rubio et al. 2001). The presence of and late stage of −P stress, respectively (Supplementary P1BS cis-acting element in the 1-kb upstream region of Pi Fig. S3). stress-responsive transcripts was therefore used as a The responsive transcripts identified based on G test also measure of response under stress. Using the signal scan included approximately 900 RAP transcripts which were search in the PLACE database (http://www.dna.affrc.go.jp/ not represented in the microarray as well as 28 unannotated PLACE/), we investigated the 1-kb upstream region of Pi transcripts revealed by the Cufflinks program (Table 2). stress-responsive representative RAP transcripts. The P1BS These unsupported RAP transcripts and unannotated tran- cis-acting elements were significantly overrepresented in scripts were classified as novel transcripts induced during – the 1-kb upstream region of transcripts upregulated in root P and ++P stress. Approximately 6.1% to 24.3% (average and shoot at late –P stress (Table 2). 56 Rice (2011) 4:50–65 and late) based on the different expression patterns (Supple- mentary Fig. S4). The GO categories showed specific trends for many Pi-related genes in response to –Pstress (Supple- mentary Fig. S5) and ++P stress (Supplementary Fig. S6). 1. Regulatory proteins/riboregulator: during –P stress, we observed specific pattern of response among genes for regulatory proteins/riboregulators such as OsIPS1, WRKY, SPX (SYG/PHO81/XPR1) domain-containing genes, histone, etc. OsIPS1 and OsIPS2 were strongly upregulated in root and shoot but PHO2/UBC24 was downregulated in shoot during late period of –P stress. The expression of OsIPS2 and PHO2/UBC24, which were not supported by the microarray, were confirmed here. The expression patterns of several transcription factors and signaling molecules under −P stress were clarified. In particular, WRKY (Os03t0321700) was upregulated in root during early −P stress. In Arabi- dopsis, AtWRKY6 and AtWRKY42 are known to modulate PHO1 transcription (Chen et al. 2009) whereas AtWRKY75 modulates Pi acquisition and root development (Devaiah et al. 2007). NAM (NAC) transcription factor (Os03t0624600) was downregu- lated in root at early stage of −P stress. The NAC gene family shows diverse functions in both plant develop- ment and stress responses as in the AtNAC1-mediated auxin signaling to promote lateral root development (Xie et al. 2000). Several SPX (SYG/PHO81/XPR1) domain-containing genes involved in response to environmental cues or internal regulation of nutrition homeostasis such as SPX1, SPX2, and SPX6 were upregulated only in shoot during late –P stress. In rice, SPX1 is partly involved in regulating the induction of some Pi transporters (Wang et al. 2009b). Histone H4 (Os04t0583600) and histone H2B (Os05t0574300) were upregulated in shoot during early –Pstress whereas histone H1 (Os03t0799000) and histone H2B Fig. 4 Distribution of upregulated and downregulated transcripts in (Os01t0152900) were downregulated in shoot during response to Pi stress. The total number of upregulated (upper)or late –P stress thereby supporting the chromatin-level downregulated (lower) transcripts identified by mRNA-Seq were regulation of response to Pi starvation (Smith et al. determined at various time points of stress (early, middle, and late) 2010). Under ++P stress treatments, OsIPS1 and during Pi starvation (−P) or phosphate over-abundant Pi supply (++P) treatments. Each bar shows the distribution of transcripts with match OsIPS2 did not show any specific pattern of expression in RAP-DB annotation and microarray (gray), transcripts with match in either root or shoot. Downregulation of auxin- in RAP-DB annotation but no match in microarray (white), and regulated transcripts such as Os03t0742900 in root at unannotated transcripts (black). early stage, Os05t0230700 in shoot at middle stage, and Os04t0519700 (auxin response factor) in root at Functional characterization of transcripts generated late stage may suggest specific functions of these genes in root and shoot growth under ++P stress. We observed by mRNA-Seq the downregulation of PTF1 in shoot at early ++P stress. Overexpression of this gene however, could enhance We used the Gene Ontology (GO) classification to estimate the functional categories of upregulated and downregulated tolerance to –P in transgenic rice (Yi et al. 2005). Among histone genes, ++P stress induced the upregula- genes in the root and shoot under various stress treatments (−P and ++P) and at different stages of stress (early, middle, tion of histone-like transcription factor (Os03t0413000) Rice (2011) 4:50–65 57 in shoot at early stage and histone H3 (Os06t0130900) 4. Metabolism (lipid): we observed specific pattern of in shoot at middle stage, and downregulation of histone response to Pi stress among transcripts related to lipid H2B (Os01t0152900) in root at late stage. and fatty-acid metabolism. During the late stage of –P 2. Metabolism (energy): many transcripts categorized as stress, there was prominent upregulation of lipid tricarboxylic acid (TCA) cycle components were induced metabolism-related genes such as UDP-sulfoquinovose during –P stress in both root and shoot, and most were synthase 1 (SQD1), SQD2, monogalactosyldiacylglycerol strongly upregulated from the middle to the late stage of – (MGDG) synthases, DGD1 (digalactosyl diacylglycerol P stress. The gene encoding isocitrate dehydrogenase 1), DGD2, and glycerophosphoryl diester phosphodies- (ICDH) was downregulated in root at the late stage of –P terase (GDPD), and transcripts encoding lipid transfer stress. As a gene involved in the degradation of citrate in proteins (Os07t0174400 and Os03t0794000). In particu- the TCA cycle, the repression of ICDH could induce an lar, SQD1, SQD2,and MGDG were upregulated in both increase in the internal citrate concentration (Delhaize et root and shoot during −P stress. Lipid and fatty-acid al. 2003). The gene encoding malate dehydrogenase metabolism-related genes involved in the synthesis of (MDH) (Os08t0434300) was upregulated in root at late galacto- and sulfolipids were strongly induced by −P stage of –P stress. In alfalfa, overexpression of MDH treatment. Lipase encoding transcripts such as caused an increase in the exudation of various organic Os01t0215000 and Os01t0710700 were upregulated in acids and subsequently resulted in increased P accumu- root at early ++P stress. Induction of these lipid lation in acid soils (Tesfaye et al. 2001). Among the OsA metabolism-related transcripts maybe associated with genes, the plasma membrane H -ATPase encoding OsA7 changes in the lipid composition and fluidity of the was strongly expressed under normal Pi concentration membranes that control metabolic and cellular processes. but the expression gradually decreased during –Pstress 5. Pi transport: among genes encoding high-affinity Pi (Supplementary Table S2), a response that maybe transporters, Pht1;2 and Pht1;8 were upregulatedinthe associated to nutrient uptake and translation during Pi roots in the middle stage whereas Pht1;3, Pht1;4, Pht1;6, starvation. During ++P stress, OsA2 was upregulated in Pht1;9,and Pht1;10 were upregulatedinthe rootsinthe root and shoot at early stage, OsA3 was upregulated in late stage of –P stress. Although high Pi concentration is shoot at middle stage, and OsA7 was downregulated in toxic to plants, the Pi transporter encoding genes were root at middle stage. Genes encoding glycolysis enzymes mostly upregulated during ++P stress including Pht1;4 were mostly upregulated as in the case of fructose- and Pht2;1 in shoot at early stage, Pht1;8 in shoot at late bisphosphate aldolase (Os11t0171300) in shoot at early stage, and Pht1;1 in root at middle stage. Only Pht1;8 stage, phosphofructokinase (Os05t0524400) in shoot at was downregulated in the root at middle stage of ++P middle stage, Os10t0405600 in shoot at late stage, and stress. These expression patterns may suggest specific enolase (Os06t0136600) in shoot at middle stage of ++P functions of these genes in Pi uptake and homeostasis stress. Enzymes involved in TCA cycle and glycolysis under conditions of over-abundant Pi supply. Among the may function to produce organic acids, which would transcripts for Pi transporters, only Pht1;11 was not permit the recycling of P from phosphorylated inter- detected (Supplementary Table S2) probably because it is mediates. Organic acids are also known to help release Pi specifically activated during mycorrhizal symbiosis from organic or insoluble inorganic Pi compounds (Paszkowski et al. 2002). outside the plant. 6. Pi remobilization: many OsRNS genes showed specific 3. Metabolism (carbohydrate): many genes encoding su- expression patterns during –P and ++P stress. It has crose metabolism-related proteins showed specific re- been suggested that RNS may have redundant functions sponse in root and shoot during –P and ++P stress. This or specific functions in different biological processes may be associated with the tightly controlled mechanisms and tissues under different biotic and abiotic stresses that allow the coordination of Pi homeostasis with carbon (MacIntosh et al. 2010). Purple acid phosphatase 10 status and photosynthesis (Wissuwa et al. 2005). Sugars (PAP10) and acid phosphatase 1 (ACP1) were upregu- generated in the shoots and transported through the lated in both root and shoot at late –P stress, which may phloem are involved in establishing a physiological, suggest a possible role in Pi acquisition and metabolism biochemical, and molecular response to −Pstress in during –P stress. plants (Hammond and White 2008). Sucrose phospha- 7. Morphological changes: From the middle to late stage tase encoding transcript (Os01t0376700) was upregu- of –P stress, genes encoding cell wall-related proteins lated in root at late –P stress whereas sucrose synthase such as pectin methylesterase 8 (Os01t0311800) and encoding transcripts (Os06t0194900 and Os03t0401300) cellulose synthase (Os07t0208500) were upregulated were downregulated in shoot and root at late–Pstress suggesting possible roles in the control of cell wall synthesis and extensibility during stress. Upregulation and root at middle ++P stress. 58 Rice (2011) 4:50–65 during the initial response to stress suggests that plant sequence (P1BS) in the 1-kb upstream or downstream growth was affected slightly during the first 10 days of – region (Supplementary Table S5). Among10,388 unanno- P stress. An auxin response factor (Os04t0519700) and tated transcripts, approximately one third (3,006 tran- an Aux/IAA binding protein (Os03t0742900) may be scripts) have amino acid sequence homology in the involved in regulating the root architecture in response to RefSeq and UniProt databases according to BLASTX −P stress by promoting lateral root elongation and search (identity, >30%; coverage, >30%) (Supplemental changing auxin distribution within root cells (Nacry et Table S6). Two unannotated Pi stress-responsive tran- al. 2005). Upregulation of anthocyanin metabolism- scripts showed homology; NP_CUFF.161528.1 to an related gene such as phenylalanine ammonia-lyase unknown protein and NP_CUFF.2572.1 to a leucine-rich (Os12t0520200) in root at late stage of –P, may be repeat family protein/protein kinase family protein. associated with its function in protecting cells from Validation experiments by qRT-PCR analyses confirmed the the stress. Many genes encoding auxin-related expression of these unannotated transcripts under −P or ++P regulatory proteins, cell wall-related proteins and stress (Fig. 5). NP_CUFF.26799.1 was upregulated in both cytoskeleton proteins were downregulated from the root and shoot at late −P stress and immediately down- early to late stage without any morphological changes for regulated after recovery from stress. NP_CUFF.141685.1 was the first 10 days during ++P stress. The genes encoding upregulated in root at late ++P stress and immediately protein synthesis-related proteins such as aminoacyl- downregulated after recovery from stress. Some of these tRNA synthetase (Os03t0749300, Os05t0150900) trans- transcripts may be regulated through the P1BS in their lation initiation factor 2 (Os03t0296400) and translation flanking regions, although they did not appear to have high elongation factor EF1B (Os06t0571400) were down- homology to one another. regulated, suggesting a repression of protein turnover and recycling of amino acids in roots and shoots. These results also show that ++P stress is a more severe form of Discussion stress for rice growth than −Pstress. 8. Homeostasis of other ions: Genes encoding transport Identification of responsive non-protein coding transcripts proteins such as potassium transporter (Os09t0448200), a under Pi stress heavy metal transporter (Os02t0530100) and an iron– sulfur cluster assembly related protein such as NifU-like We used the mRNA-seq approach to characterize the protein (Os01t0662600) showed specific expression transcriptome of rice under Pi stress at high resolution. patterns that may be associated to their functions in This is the first report of a whole transcriptome analysis involving stress due to Pi starvation and over-abundant Pi adjusting the nutrient balance under –P stress. As Pi is important in numerous energy-requiring metabolic and supply. In addition to the confirmation of almost 77% of transport processes, genes encoding transport proteins RAP transcripts, the results also provide evidence for such as heavy metal transporter (Os07t0258400, expression of transcripts not represented in the microarray Os01t0125600, and Os01t0976300) were upregulated including transcripts without FL-cDNA support, transcripts during ++P stress. based on gene prediction programs, or those identified by homology to other plants (Fig. 3). Furthermore, the Pi stress-responsive transcripts can be functionally character- Characterization of unannotated transcripts induced by −P ized based on the expression pattern in root and shoot at and ++P stress early, middle, or late stage of −P or ++P stress. A total of 10,388 transcripts identified using the Cufflinks program The unannotated Pi-responsive transcripts validated by G test have not yet been previously annotated. Some of these a ranged in length from 135 to 1,750 bp with an average of unannotated transcripts may have been induced by alterna- 600 bp (Supplementary Table S4). These transcripts were tive promoters as they have been identified in the promoter distributed in the 12 chromosomes. The nearest adjacent region of RAP transcripts, within the 5′ and 3′ exons transcripts were located from a distance of 1 to more than supported by ESTs or FL-cDNAs, or within the extended 3,000 bp. Most transcripts showed specific expression regions of the exons. Although unannotated transcripts may pattern in root and shoot at early, middle or late stage also be found in the introns, we could not differentiate of −P and ++P stress. Two unannotated transcripts, intron-derived transcripts from newly formed exon tran- namely, NP_CUFF.172640.1 and NP_CUFF.136416.1, were script and the mRNA-seq approach could not define found to contain an intron which showed their direction based whether the RNA is sense or antisense. In order to on the junction sequence. Analysis of flanking regions also characterize intron-derived transcripts, it maybe necessary revealed that most of these transcripts contain a PHR1-binding to obtain longer reads and pair-end sequence data. Rice (2011) 4:50–65 59 NP_CUFF.111084.1 1200 ++P ++P ++P * * 600 20 2 10 0 0 0 15 10 1015 10 10 15 10 1015 10 10 15 10 1015 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day NP_CUFF.88762.1: Root_++P_Early_Down NP_CUFF.141685.1: Root_++P_Late_Up 2.5 ++P ++P ++P 0.8 0.6 1.5 0.4 0.2 0.5 * 0 0 0 15 10 1015 10 10 15 10 15 10 15 10 15 10 10 10 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day ++P ++P ++P 25 * * 20 10 * 0 0 15 10 1015 10 10 15 10 1015 10 10 15 10 1015 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day NP_CUFF.2572.1: Shoot_++P_Early_Down NP_CUFF.76959.1: Shoot_++P_Early_Down 1.4 0.8 1.2 ++P ++P ++P 0.6 0.8 0.4 0.4 0.2 0.2 15 10 15 10 15 10 1015 10 10 10 10 15 10 1015 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day Fig. 5 Validation of responsive unannotated transcripts. The expression expression on day 10 during –P/++P and at 1 day after recovery from of OsIPS1 and several unannotated transcripts in response to –Pand ++P stress (Re1). The bars represent the mean±SE level of transcripts from was validated by qRT-PCR analysis. OsIPS1 is a –P upregulation index three experiments. The asterisk indicates the responsive point based on G gene with statistical significance in the roots (upper)and shoots (lower) test (FDR, <0.01). The primers used in qRT-PCR analysis for as determined by G test. Transcript expression levels were normalized unannotated transcripts are described in Supplementary Table S5. using an internal control (Ubiqutin1) and plotted relative to the We have identified a total of 10,388 unique unannotated various Pi stress treatments (Supplementary Table S4). transcripts from the root and shoot with 3,006 transcripts These transcripts comprised less than 1% of the responsive showing homology to various proteins (Supplementary transcripts and had much lower expression levels than Table S6). Among them, 28 transcripts have been found RAP-representative transcripts for all conditions. For to be rice-specific and 303 transcripts showed homology to example, the average RPKM values (±SD) for unannotated other known plant genes. The remaining unannotated and RAP-representative transcripts in the root control transcripts (7,382) without homology to any protein may sample were 0.8 (±2.6) and 23.8 (±94.5), respectively. As include non-protein coding transcripts, novel protein tran- these transcripts also showed small differential expression scripts, rare transcripts that are expressed at low copies, between conditions that could not be detected by statistical transcripts with very low expression levels, or even tests, it maybe difficult to characterize all these transcripts. transcripts which may have lethal functions in Escherichia Two unannotated transcripts showed protein homology to coli. We particularly focused on the 28 unannotated known proteins, one to leucine family gene and another to transcripts, which were differentially expressed under an unknown protein. Furthermore, the presence of P1BS Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression 60 Rice (2011) 4:50–65 cis-acting elements was confirmed in the promoter regions gene expression between the root and shoot (Supplemen- of ten unannotated transcripts. These findings establish the tary Fig. S7). These morphological and physiological utility of mRNA-Seq in identifying uncharacterized tran- changes are associated with the changes in gene expression. scripts that are expressed in response to Pi stress. The differences in gene expression between −P and ++P The expression of approximately 70% of the responsive stress in root and shoot (Supplementary Figure S7a, b)as unannotated transcripts was confirmed by qRT-PCR anal- well as the differences in gene expression between root and ysis. Several transcripts (NP_CUFF.26799.1 and shoot under −Pand ++Pstress (Supplementary Figure S7c, NP_CUFF.141685.1) showed a sharp and reversible re- d) suggest that Pi stress treatments have a significant effect sponse to −P or ++P stress (Fig. 5). Most of these on the overall expression profile of many genes involved in transcripts do not contain any ORFs and may function as plant response to stress. We have confirmed the response of non-coding transcripts similar to OsIPS1. Some responsive many previously reported phosphate-related genes. Although unannotated transcripts could not be validated because of some of the genes that showed significant changes in very low expression levels that could not be detected by expression may also be involved in the developmental stage, qRT-PCR. Overall, these validation experiments further there is no doubt that most of the responsive genes were suggest that mRNA-Seq can provide strong evidence of expressed due to Pi stress as shown by the statistical transcript expression. It has also been efficiently used in confirmation of the differences in gene expression over time confirming the function of miRNAs such as miR399 in during −P and ++P treatments. After recovery from stress, response to Pi stress. In Arabidopsis, the interaction we also found that many genes upregulated or down- between the Pi- and nitrate-limited signaling pathways regulated genes during Pi stress, including OsIPS1,reverted affecting E3 ligase gene NLA and anthocyanin synthesis to the same level of expression before stress treatments. was found to be mediated by miR827, which targets the Among the –P responsive transcripts, 316 transcripts were transcripts of specific proteins that contain the SPX domain commonly expressed in both root and shoot. These transcripts (Pant et al. 2009). The rice homolog OsmiR827 exhibits likely function in basic –P signaling cascade that impacts Pi different preferences in tissue expression and regulates transport and Pi remobilization in plants, and included different classes of genes sharing a common SPX domain upregulated transcripts such as OsIPS1, OsIPS2, SQD1, with AtmiR827 (Lin et al. 2010). The existence of miRNA- OsRNS3,and PAP10. On the other hand, PHO2/UBC24, dependent pathway in distinct aspects of Pi homeostasis in which regulates Pi allocation between roots and shoots, was these two species may suggest a species-specific pathway downregulated in the shoots. A total of 1,336 –Presponsive where many unannotated transcripts may play significant transcripts in the roots were not responsive to in the shoots roles. It has been reported that miR2111 may target the E3 (Supplementary Fig. S7c). These transcripts likely function ligase gene and a calcineurin-like phosphoesterase gene in the cell-rescue system that impacts Pi uptake, carbon- under −P stress (Hsieh et al. 2009). Additionally, miRNAs assimilation reduction and root-system morphology by such as miR169, miR395, and miR398, which respond to enhancing root growth, and included upregulated transcripts various nutritional stresses, were also downregulated in that encode many phosphate transporters (Pht1;2, Pht1;3; response to Pi stress (Hsieh et al. 2009). Therefore, some of Pht1;6, Pht1;8, Pht1;9, and Pht1;10), an enzyme in the TCA the unannotated transcripts detected here including non- cycle (PEPC/Os08t0366000) and an auxin response factor protein coding transcripts and small RNAs may have (Os04t0519700). It has been reported that Pht1;8 and Pht1;9 important functions in the signaling of plant root responses were located downstream of the PHO2/UBC24-dependent to changing Pi concentrations and the coordination of signaling pathway (Bari et al. 2006). homeostatic pathways of Pi and other nutrients. A comparison of the ++P responsive transcripts in root and shoot showed that 358 transcripts were commonly expressed. Responsive genes related to morphological However, many transcripts with different functions responded and physiological changes in root and shoot under ++P stress. In particular, PHO2/ UBC24 was upregulated in root at early ++P stress. This Here, we found that –P stress inriceisassociatedwithseveral expression may be sufficient to repress the –P upregulated morphological changes such as brown coloration, roughening transcripts such as Pht1;8, which is downstream of the of the roots and a decreased ratio of shoot-to-root weight PHO2/UBC24-dependent signaling pathway. Pht1;1 and (Fig. 1). Under ++P stress, we also observed weak roots and Pht1;4 were also upregulated in root during ++P, which as shoots, and an increased ratio of shoot-to-root weight with suggested by Bari et al. (2006) may not be a part of the increased Pi uptake. Interestingly, although the ratio of PHO2/UBC24-dependent signaling pathway. There were shoot-to-root weight clearly reversed between 10 and 15 days 2,248 ++P responsive transcripts in shoot that were not in the –P and ++P stress-treated plants, the weights of the responsive in root. These included upregulated genes such as plants decreased concordantly as reflected in differences in Pht1;8 and Pht2;1 which may contribute to the establish- Rice (2011) 4:50–65 61 ment of high Pi in shoot. Pht1;8 may also be upregulated by starvation-induced genes (Bustos et al. 2010). Therefore, other factors without downregulating PHO2/UBC24 tran- PHR1 is maybe upregulated in late stages of −P stress via a scription. Protein synthesis (Os01t0678600) and fatty-acid P1BS cis-acting element. Since the rest of the genes metabolism (Os02t0205500) transcripts were downregulated. upregulated in late stage of −P stress do not have a P1BS A large number of responsive transcripts detected in shoot at cis-acting element in the promoter regions, other unknown early ++P stress may have diverse functions in Pi metabo- cis-acting elements may also be involved in upregulation in lism, photosynthesis and ATP biosynthesis, and may be late stage of −P stress as well as in other expressions transiently activated for plant growth. At the same time, high patterns in response to Pi stress. Some transcripts may be Pi concentrations may inhibit plant growth because the genes regulated through the release from repression of PHO2/ encoding many nuclear proteins, such as poly(A) polymerase UBC24 or regulated by other factors such as PHL1. Both (e.g., Os06t0558700), many protein synthesis-related pro- PHR1 and PHL1 were found to be partially redundant and teins, such as translation initiation factor 2 (e.g., have a central role in the control of physiological and Os03t0296400) and aminoacyl-tRNA synthetases (e.g., molecular responses to –P stress (Bustos et al. 2010). Os06t0645400) were downregulated in shoot. The plant We searched gene loci with two or more alternative may then acclimate and use less energy because the number splicing isoforms supported by FL-cDNAs to determine of responsive transcripts decreases from the middle to the whether alternative splicing mechanisms affect gene ex- late stage of stress in shoot. Interestingly, Pht1;4 was pression during Pi stress. Among 34,780 representative loci upregulated in root and shoot under both –Pand ++P stress. on the rice genome (IRGSP Build 5 Pseudomolecule), This finding suggests that Pht1;4 can easilyadjusttothe alternative splicing isoforms were confirmed in 5,625 cellular Pi condition. In Arabidopsis, it has been reported (16.2%) representative loci. To examine changes in the that AtPHT1;4 was induced by sucrose and is thought to be usage of alternative splicing isoforms against Pi stress, the upstream of the hexokinase sugar-sensing pathways (Lejay et isoforms were compared in terms of Pi stress-responsive al. 2008). Thus, morphological changes for adaptation in rice expression patterns. As a result, the alternative splicing due to changes in Pi concentration may be controlled by isoforms of 194 loci (0.6%) showed different expression local and systemic gene expression. patterns. This suggests that the Pi stress-responsive expres- The differences in total phosphorus content of rice plants sion is also regulated by alternative splicing mechanisms grown under normal Pi concentration (+P), Pi starvation (−P), with respect to the timing of the response or tissue and over-abundant Pi supply (++P) are direct manifestations of specificity. Further analysis of the upstream regions of Pi transporter activities. We identified several responsive responsive transcripts and accompanying changes in alter- transporter genes such as Pht1;4, Pht1;8, and Pht2;1. Genetic native splicing may reveal expression networks underlying modification of Pi transport proteins are maybe the key in Pi homeostasis in plants. developing tolerance to both –P and ++P stress. Britto and Kronzucker (2008) reported that the high-affinity transport A genome-wide transcriptional analysis of Pi stress system is strongly downregulated under K -replete conditions and is conversely strongly upregulated under K -starvation This study has shown that RNA-Seq can be used to determine conditions and that the LATS (low-affinity transport system)- transcript expression levels more accurately and capture the + + range K influx may be increased generally by increased K transcriptome dynamics associated with Pi stress across (Britto and Kronzucker. 2008). We also observed upregula- different time points and tissues. Interestingly, SIZ1, PHR1, tion of several low- and high-affinity transporter genes such and PHR2, which are located upstream of OsIPS1,were as Pht1;1, Pht1;8, and Pht2;1 under ++P stress. Therefore, a slightly downregulated after 10 days under −Pstress in roots portion of the Pi influx may be regulated by a transport and shoots and slightly upregulated at 1 d after recovery in system similartothat ofK influx. roots, an expression pattern which was in contrast with that of OsIPS1. These genes were previously reported to be Promoter analysis and alternative splicing isoforms upregulated slightly under −P stress based on qRT-PCR of responsive genes analysis and northern analysis (Rubio et al. 2001;Miura et al. 2005). This apparent discrepancy indicates that a more The flanking regions of most Pi stress-responsive tran- detailed analysis of gene expression during different stages scripts contain a P1BS cis-acting element which were of stress via mRNA-Seq could be useful in detecting subtle significantly over represented in the 1-kb upstream regions changes in gene expression. of 170 transcripts (33.8%) from the roots and 74 transcripts Using mRNA-Seq, we were able to identify many novel (34.3%) from the shoots that were upregulated at the late transcripts which could not be revealed by microarray- stage of −P stress. In Arabidopsis, the direct targets of based technology. We were also able to characterize even subtle changes in their expression at different stages of PHR1 are greatly enriched in P1BS-containing Pi 62 Rice (2011) 4:50–65 stress under –P and ++P treatments. These results will be (Life Technologies, Rockville, MD, USA). The RNA was useful in a detailed functional analysis of an accurate gene treated with DNase I (Takara, Shiga, Japan), and the first- set provided in the high-quality map-based genome strand cDNA was synthesized using the SuperScript® III sequence of the model rice cultivar Nipponbare. Towards first-strand synthesis system (Invitrogen, Carlsbad, CA, this goal, we also constructed transcriptome viewer in USA) according to the manufacturer’s protocol. The GBrowse format to facilitate visualization of all uniquely cDNAs produced by reverse transcription were amplified mapped mRNA-Seq reads in the rice genome under Pi using a pair of gene-specific primers for each gene. stress (Fig. 6; http://rnaseq-pi.dna.affrc.go.jp/cgi-bin/gb2/ gbrowse/RNAseq_p/). A comprehensive rice transcriptome Sequencing and mapping of short reads onto the rice under Pi stress could be an indispensable tool in co- genome expression analysis of related genes, deciphering gene networks involved in Pi stress, and identifying genes that Extraction of total RNA, construction of cDNA, and maybe used to enhance P efficiency and improving P sequencing with the Illumina Genome Analyzer IIx were acquisition, which are all critical to crop productivity. performed as described previously (Mizuno et al. 2010). Trimming of Illumina adaptor sequences and low-quality bases (Q<20) at the 5′ and 3′ ends of each read was Materials and methods performed by the fastx_clipper program in the FASTX- Toolkit and an in-house C language program. Using the Plant materials and stress treatment pre-processed Illumina reads, the transcript structures were reconstructed by a series of programs, namely, the Bowtie Rice (Oryza sativa ssp. japonica cv. Nipponbare) seeds were version 0.12.4 for short-read mapping (Langmead et al. germinated and grown by hydroponic culture in nutrient 2009), TopHat ver version 1.0.13 for defining exon–intron media (Yoshida et al. 1976) in a growth chamber. After junctions (Trapnell et al. 2009), and Cufflinks version 0.8.2 14 days, the seedlings were subjected to Pi stress treatment by for gene structure predictions (Trapnell et al. 2010). transferring in a similar media with the following Pi Unannotated transcribed regions with no overlap to any concentrations: 0 mg Pi/L (18.7 mg NaCl/L; alternatively known loci was analyzed by comparing known transcript NaH PO )for –P stress, and 100 mg Pi/L (439 mg structures annotated in RAP-DB (http://rapdb.dna.affrc.go. 2 4 KH PO /L; alternatively NaH PO ) for ++P stress. The jp/) and the MSU Rice Genome Annotation Project 2 4 2 4 plants were maintained under Pi stress conditions for 30 days database (http://rice.plantbiology.msu.edu/) with the at 28°C and 70–80% humidity in a 16-h light/8-h dark cycle Cufflinks-predicted structures. Those regions that over- with the light period from 6:00 AM to 10:00 pm. The plants lapped with rRNA (1 transcript), tRNA (12 transcripts), were moved to different positions in the growth chamber and repetitive regions with similarity of ≥90% (3,888 several times during the light period so that all plants in each transcripts) in RAP-DB were discarded, and the remaining treatment were exposed to the same amount of light. Root unannotated transcripts were used for further analysis. To and shoot samples were collected at approximately 9:00 AM, estimate the expression levels of each transcript, all pre- frozen in liquid nitrogen, and stored at −80°C until processed reads were mapped into the IRGSP Build 5 subsequent analyses. pseudomolecules (http://rgp.dna.affrc.go.jp/E/IRGSP/ Build5/build5.html) by BWA (version 0.5.7) with default Measurement of dry weight and P uptake parameters (Li and Durbin. 2009). The expression level for each known and unannotated transcript was calculated as Shoots and roots of seedlings under +P (control), −P, and ++P RPKM (Reads Per Kilobase exon Model per Million conditions were separately sampled and oven-dried at 70°C mapped reads) based on the number of uniquely mapped for 2 days. The dry weight of each sample was determined and reads that completely overlap with the exonic regions. The the phosphorus content was analyzed by the phosphomolyb- resulting mRNA-Seq data were deposited in the DDBJ denum blue reaction using the U-0080D spectrophotometer Sequence Read Archive (accession no. DRA000314). (Hitachi, Tokyo, Japan) and Spectroquant phosphate test kit (Merck, Darmstadt, Germany). Rice 44 K microarray analysis Time course analysis of OsIPS1 expression Root RNA samples from plants under under+P (control), −P, and ++P conditions were used for gene expression profiling The expression of OsIPS1 was analyzed by RT-PCR. Total using the rice 44 K oligomicroarray (Agilent Technologies). RNA was extracted from root and shoot samples of 1–10- Labeling, hybridization and scanning were performed accord- ing to the manufacturer’s protocol as described previously day-old seedlings under –P condition using Trizol® reagent Rice (2011) 4:50–65 63 Fig. 6 Browser of the rice tran- scriptome under –P and ++P stress. The mRNA-Seq data mapped on the IRGSP Build 5 pseudomolecules can be accessed in a GBrowse format. The position of the transcript and corresponding microarray probe (if available) is indicated. For each transcript, the expres- sion profiles of mapped reads including RAP transcripts (e.g., OsIPS1) and unannotated tran- script under normal Pi concen- tration (control), –P and ++P stress, and recovery condition are shown as graphs. (Mizuno et al. 2010). Three biological replicates for each Analysis of responsive transcripts under Pi stress treatment were used for analysis. All the analysis data from microarray experiments are deposited in the NCBI Gene To detect stress response transcripts from all RAP transcripts Expression Omnibus (accession no. GSE25647). and unannotated transcripts, statistical analysis by G test was conducted between two stages of stress treatments for each Validation of gene expression by qRT-PCR transcript. The numbers of mapped reads on a given transcript and those on other regions for two stages were used as Expression of unannotated transcripts was analyzed by variables in 2×2 contingency tables for each test. All p values quantitative RT-PCR using the LightCycler® 480 system were corrected for a FDR (Benjamini and Hochberg. 1995) (Roche, Basel, Switzerland). Total RNA from root and shoot of 1% using the R package (ver 2.7.1) and in-house Perl was used for reverse transcription using the Transcriptor First- scripts (Benjamini and Hochberg. 1995). The GO terms Strand synthesis kit. The reaction mixture and reaction assigned to each transcript were obtained from RAP-DB for condition was described previously (Mizuno et al. 2010). each GO category for biological process, molecular function, The detection threshold cycle for each reaction was and cellular component. We also searched for the P1BS cis- normalized using Ubiquitin1 with 5′-CCAGGACAAGAT acting element in the 1-kb upstream regions of Pi stress- GATCTGCC-3′ and 5′-AAGAAGCTGAAGCATCCAGC-3′ responsive RAP-representative transcripts using the PLACE as primers. Three technical replicates for each treatment were database (http://www.dna.affrc.go.jp/PLACE/) and signal used for analysis. scan search program (Prestridge 1991). 64 Rice (2011) 4:50–65 Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory- Construction of transcriptome viewer efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25. We constructed a browser to show the distribution of Lejay L, Wirth J, Pervent M, Cross JMF, Tillard P, Gojon A. mapped mRNA-Seq reads for each condition along with the Oxidative pentose phosphate pathway-dependent sugar sensing as a mechanism for regulation of root ion transporters by reference genome sequence and RAP annotation in a photosynthesis. Plant Physiol. 2008;146:2036–53. GBrowse format Li H, Durbin R. Fast and accurate short read alignment with Burrows– Wheeler transform. Bioinformatics. 2009;25:1754–60. Li LH, Qiu XH, Li XH, Wang SP, Lian XM. The expression profile of Acknowledgments We thank Ms. F. Aota, Ms. K. Ohtsu, and Ms. K. genes in rice roots under low phosphorus stress. Sci China Ser C. Yamada for technical assistance, and Dr. Y. Nagamura and Ms. R. 2009;52:1055–64. Motoyama for assistance with the microarray experiments. This work Lin SI, Chiou TJ, Santi C, Jobet E, Lacut E, El Kholti N, et al. was supported by the Ministry of Agriculture, Forestry and Fisheries of Complex regulation of two target genes encoding SPX-MFS Japan (Genomics for Agricultural Innovation, RTR-0001). proteins by rice miR827 in response to phosphate starvation. Plant Cell Physiol. 2010;51:2119–31. Open Access This article is distributed under the terms of the Creative Lu TT, Lu GJ, Fan DL, Zhu CR, Li W, Zhao QA, et al. Function Commons Attribution Noncommercial License which permits any annotation of the rice transcriptome at single-nucleotide resolu- noncommercial use, distribution, and reproduction in any medium, tion by RNA-seq. Genome Res. 2010;20:1238–49. provided the original author(s) and source are credited. MacIntosh GC, Hillwig MS, Meyer A, Flagel L. RNase T2 genes from rice and the evolution of secretory ribonucleases in plants. Mol Genet Genomics. 2010;283:381–96. Miura K, Rus A, Sharkhuu A, Yokoi S, Karthikeyan AS, Raghothama KG, et al. The Arabidopsis SUMO E3 ligase SIZ1 controls phosphate References deficiency responses. Proc Natl Acad Sci USA. 2005;102:7760–5. Mizuno H, Kawahara Y, Sakai H, Kanamori H, Wakimoto H, Bari R, Datt Pant B, Stitt M, Scheible WR. PHO2, microRNA399, and Yamagata H, et al. Massive parallel sequencing of mRNA in PHR1 define a phosphate-signaling pathway in plants. Plant identification of unannotated salinity stress-inducible transcripts Physiol. 2006;141:988–99. in rice (Oryza sativa L.). BMC Genomics. 2010;11:683. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a Morcuende R, Bari R, Gibon Y, Zheng W, Pant BD, Blasing O, et al. practical and powerful approach to multiple testing. R Statist Soc Genome-wide reprogramming of metabolism and regulatory B. 1995;57:289–300. networks of Arabidopsis in response to phosphorus. Plant Cell Britto DT, Kronzucker HJ. Cellular mechanisms of potassium Environ. 2007;30:85–112. transport in plants. Physiol Plantarum. 2008;133:637–50. Mortazavi A, Williams BA, Mccue K, Schaeffer L, Wold B. Mapping Bustos R, Castrillo G, Linhares F, Puga MI, Rubio V, Perez-Perez J, et and quantifying mammalian transcriptomes by RNA-Seq. Nat al. A Central regulatory system largely controls transcriptional Methods. 2008;5:621–8. activation and repression responses to phosphate starvation in Nacry P, Canivenc G, Muller B, Azmi A, Van Onckelen H, Rossignol Arabidopsis. Plos Genet. 2010;6(9):e1001102. M, et al. A role for auxin redistribution in the responses of the Chen YF, Li LQ, Xu Q, Kong YH, Wang H, Wu WH. The WRKY6 root system architecture to phosphate starvation in Arabidopsis. transcription factor modulates PHOSPHATE1 expression in response Plant Physiol. 2005;138:2061–74. to low Pi stress in Arabidopsis. Plant Cell. 2009;21:3554–66. Pant BD, Musialak-Lange M, Nuc P, May P, Buhtz A, Kehr J, et al. Delhaize E, Ryan PR, Hocking PJ, Richardson AE. Effects of altered Identification of nutrient-responsive Arabidopsis and rapeseed citrate synthase and isocitrate dehydrogenase expression on microRNAs by comprehensive real-time polymerase chain internal citrate concentrations and citrate efflux from tobacco reaction profiling and small RNA sequencing. Plant Physiol. (Nicotiana tabacum L.) roots. Plant Soil. 2003;248:137–44. 2009;150:1541–55. Devaiah BN, Karthikeyan AS, Raghothama KG. WRKY75 transcription Paszkowski U, Kroken S, Roux C, Briggs SP. Rice phosphate factor is a modulator of phosphate acquisition and root development transporters include an evolutionarily divergent gene specifically in Arabidopsis. Plant Physiol. 2007;143:1789–801. activated in arbuscular mycorrhizal symbiosis. Proc Natl Acad Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio- Sci USA. 2002;99:13324–9. Somoza I, et al. Target mimicry provides a new mechanism for Prestridge DS. SIGNAL SCAN: a computer program that scans DNA regulation of microRNA activity. Nat Genet. 2007;39:1033–37. sequences for eukaryotic transcriptional elements. Comput Appl Fujii H, Chiou TJ, Lin SI, Aung K, Zhu JK. A miRNA involved in Biosci. 1991;7:203–6. phosphatestarvation response in Arabidopsis.CurrBiol. Rubio V, Linhares F, Solano R, Martin AC, Iglesias J, Leyva A, et al. 2005;15:2038–43. A conserved MYB transcription factor involved in phosphate Gojon A, Nacry P, Davidian JC. Root uptake regulation: a central starvation signaling both in vascular plants and in unicellular process for NPS homeostasis in plants. Curr Opin Plant Biol. algae. Gene Dev. 2001;15:2122–33. 2009;12:328–38. Smith AP, Jain A, Deal RB, Nagarajan VK, Poling MD, Raghothama Hammond JP, White PJ. Sucrose transport in the phloem: integrating root KG, et al. Histone H2A.Z regulates the expression of several responses to phosphorus starvation. J Exp Bot. 2008;59:93–109. classes of phosphate starvation response genes but not as a He GM, Zhu XP, Elling AA, Chen LB, Wang XF, Guo L, et al. Global transcriptional activator. Plant Physiol. 2010;152:217–25. epigenetic and transcriptional trends among two rice subspecies Tesfaye M, Temple SJ, Allan DL, Vance CP, Samac DA. Overexpression and their reciprocal hybrids. Plant Cell. 2010;22:17–33. of malate dehydrogenase in transgenic alfalfa enhances organic acid Hsieh LC, Lin SI, Shih ACC, Chen JW, Lin WY, Tseng CY, et al. synthesis and confers tolerance to aluminum. Plant Physiol. Uncovering small RNA-mediated responses to phosphate defi- 2001;127:1836–44. ciency in Arabidopsis by deep sequencing. Plant Physiol. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice 2009;151:2120–32. junctions with RNA-Seq. Bioinformatics. 2009;25:1105–11. Rice (2011) 4:50–65 65 Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren Wissuwa M, Gamat G, Ismail AM. Is root growth under phosphorus MJ, et al. Transcript assembly and quantification by RNA-Seq deficiency affected by source or sink limitations? J Exp Bot. reveals unannotated transcripts and isoform switching during cell 2005;56:1943–50. differentiation. Nat Biotechnol. 2010;28:511–U174. Xie Q, Frugis G, Colgan D, Chua NH. Arabidopsis NAC1 transduces Vance CP, Uhde-Stone C, Allan DL. Phosphorus acquisition and use: auxin signal downstream of TIR1 to promote lateral root critical adaptations by plants for securing a nonrenewable development. Gene Dev. 2000;14:3024–36. resource. New Phytol. 2003;157:423–47. Yi K, Wu Z, Zhou J, Du L, Guo L, Wu Y, et al. OsPTF1, a novel Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcription factor involved in tolerance to phosphate starvation transcriptomics. Nat Rev Genet. 2009a;10:57–63. in rice. Plant Physiol. 2005;138:2087–96. Wang C, Ying S, Huang H, Li K, Wu P, Shou H. Involvement of Yoshida S, Forno AD, Cock HJ, Gomez AK. Laboratory manual for OsSPX1 in phosphate homeostasis in rice. Plant J. 2009b;57:895– physiological studies of rice. 3rd ed. Manila: International Rice 904. Research Institute; 1976. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Rice Springer Journals

mRNA-Seq Reveals a Comprehensive Transcriptome Profile of Rice under Phosphate Stress

Download PDF
16 pages

Loading next page...
 
/lp/springer-journals/mrna-seq-reveals-a-comprehensive-transcriptome-profile-of-rice-under-u5MONvHUuE

References (43)

Publisher
Springer Journals
Copyright
Copyright © 2011 by The Author(s)
Subject
Life Sciences; Plant Sciences; Plant Genetics & Genomics; Plant Breeding/Biotechnology; Agriculture; Plant Ecology
ISSN
1939-8425
eISSN
1939-8433
DOI
10.1007/s12284-011-9064-0
Publisher site
See Article on Publisher Site

Abstract

Rice (2011) 4:50–65 DOI 10.1007/s12284-011-9064-0 mRNA-Seq Reveals a Comprehensive Transcriptome Profile of Rice under Phosphate Stress Youko Oono & Yoshihiro Kawahara & Hiroyuki Kanamori & Hiroshi Mizuno & Harumi Yamagata & Mayu Yamamoto & Satomi Hosokawa & Hiroshi Ikawa & Ikuko Akahane & Zuofeng Zhu & Jianzhong Wu & Takeshi Itoh & Takashi Matsumoto Received: 5 July 2011 /Accepted: 17 October 2011 /Published online: 17 November 2011 The Author(s) 2011. This article is published with open access at Springerlink.com Abstract Plants have developed several morphological and detected uniquely 10,388 transcripts with no match to any physiological strategies to adapt to phosphate stress. We RAP transcript. These transcripts that showed specific analyzed the inducible transcripts associated with phosphate response to Pi stress include those without ORFs that may starvation and over-abundant phosphate supply to characterize act as non-protein coding transcripts. With an accompanying the transcriptome in rice seedlings using the mRNA-Seq browser of the transcriptome under Pi stress, a deeper strategy. Fifty-three million reads obtained from 16 libraries understanding of the structural and functional features of both under various phosphate stress and recovery treatments were annotated and unannotated Pi stress-responsive transcripts can uniquely mapped to the rice genome. Transcripts identified provide useful information in improving Pi acquisition and specifically tagged to 40,574 (root) and 39,748 (shoot) Rice utilization in rice and other cereal crops. Annotation Project (RAP) transcripts. Additionally, we . . Keywords Phosphorus Phosphate starvation Over- . . . abundant phosphate supply Rice mRNA-Seq Electronic supplementary material The online version of this article Transcriptome (doi:10.1007/s12284-011-9064-0) contains supplementary material, which is available to authorized users. : : : : Y. Oono Y. Kawahara H. Kanamori H. Mizuno Introduction : : : : : H. Yamagata M. Yamamoto S. Hosokawa Z. Zhu J. Wu T. Itoh T. Matsumoto (*) Phosphorus (P) is a component of key cellular molecules Plant Genome Research Unit, Agrogenomics Research Center, National Institute of Agrobiological Sciences (NIAS), such as nucleic acids, proteins, phospholipids, phytic acid, 2-1-2 Kannondai, and ATP in plants. Plants absorb P almost exclusively in an Tsukuba, Ibaraki 305-8602, Japan inorganic phosphate form (Pi), primarily as H PO from 2 4 e-mail: mat@nias.affrc.go.jp the soil. Among the major nutrients necessary for plant H. Ikawa growth, P is the most dilute and the least mobile in soil; Forestry and Fisheries, Institute of the Society hence it is often a limiting factor for crop yield. Pi for Techno-innovation of Agriculture, starvation during farming is alleviated by the massive 446-1 Ippaizuka, Kamiyokoba, application of fertilizers. However, continuous usage of Tsukuba, Ibaraki 305-0854, Japan phosphate fertilizers may have a negative impact on the I. Akahane environment as the rock phosphate in the world is now in National Institute for Agro-Environmental Sciences (NIAES), short supply and maybe depleted within the next century 3-1-3 Kannondai, (Vance et al. 2003). On the other hand, massive Pi Tsukuba, Ibaraki 305-8604, Japan fertilization is also known to inhibit plant growth and Present Address: productivity. A thorough understanding of plant response to Z. Zhu stress due to Pi starvation (−P) as well as an over-abundant Department of Plant Genetics and Breeding, Pi supply (++P) is therefore indispensable in improving China Agricultural University, acquisition and utilization of this nutrient. Beijing 100093, China Rice (2011) 4:50–65 51 The molecular mechanisms by which plants respond and comprehensive overview of the primary molecular events acclimate to changes in the nutritional Pi concentration are resulting from phosphate stress. complex but of great importance and could be useful in developing strategies for elucidating the gene networks involved in plant response to various kinds of abiotic stress. Results The regulation system for the—P signaling pathway in plants has been proposed using studies in Arabidopsis and Morphological changes triggered by −P and ++P stress involved the sumoylation of the MYB transcription factor PHR1 by a small ubiquitin-like modifier SIZ1 in a process After the first 10 days in −P and ++P stress treatments, the dependent on E3 ligase (Miura et al. 2005; Gojon et al. rice seedlings began to show morphological changes in 2009). Located downstream of PHR1, miR399 and IPS1 roots and shoots that may be attributed to starvation or (induced by phosphate starvation1) are specifically induced over-abundant supply of phosphate. These changes gradu- by −P stress. The miR399 reciprocally regulates the ally became prominent so that after 30 days of –P stress, expression of PHO2/UBC24 and the resulting protein shoot growth was retarded with thin and fewer leaves functions as an ubiquitin-conjugating E2 enzyme (Bari et whereas the root became bristlier in texture and more al. 2006). A loss of function of PHO2/UBC24 could lead to brownish in color as compared with the control (Fig. 1a). excessive Pi accumulation in the shoot (Fujii et al. 2005). Under ++P stress however, the shoot gradually showed IPS1, a member of the Mt4/TPS1 family, was shown to signs of wilting with some yellowing of the leaf and mimic the target miR399 and offset its function (Franco- relatively inhibited root growth. The total phosphorus Zorrilla et al. 2007). The sequences near the miR399 content in −P stress-treated plants was lower than the complementary regions in different plant species, and all control for both shoots and roots (Fig. 1b). On the other known members of the IPS1 family in Arabidopsis, are hand, in ++P stress-treated plants, the total P content was highly conserved. Of note, if Pi is supplied to plants, IPS1 1.3-fold higher in the shoots but slightly lower in the roots transcripts rapidly disappear suggesting that they are than the control. The dry weight of −P and ++P stress- component riburegulators of the Pi signaling system. treated plants decreased gradually to 38.4% and 45.3% of PHR1 also upregulates many responsive genes that are the control plants after 30 days, respectively (Fig. 1c). The expressed under −P stress, including those encoding high- shoot/root weight ratio differed in the two stress treatments affinity Pi transporters, RNases, and acid phosphatases, after 10–15 days with a decrease in the −P stress-treated by binding to an imperfect palindromic sequence (P1BS, plants and an increase in ++P stress-treated plants as GNATATNC) in the promoter regions of these genes compared with the control (Fig. 1d). Although the roots (Rubio et al. 2001). However, although Pi stress due to weight did not decrease much under −P stress, there was a starvation has been characterized in detail, not much is significant decrease in the dry weight of the roots and known on the signaling pathway involved in Pi stress shoots under ++P stress. These results indicate that rice brought about by an over-abundant supply of phosphate in seedlings respond to stress due Pi starvation or an over- the soil. abundant Pi supply with visible changes in root and shoot Recently, the mRNA-Seq strategy using next generation morphology after 10 days of treatment. sequencers has become a useful tool for analyzing genome- wide gene expression and cataloguing of all transcripts Time course expression of OsIPS1 under P stress including mRNAs, non-coding RNAs and small RNAs. With a high resolution and sensitivity, the mRNA-Seq can We next analyzed the expression of OsIPS1, a Pi starvation provide detailed information on transcriptional structure of responsive non-protein coding regulatory transcript, to genes such as the precise location of transcription bound- confirm the validity of −P stress treatment. OsIPS1 was aries to a single-base resolution, reveal rare transcripts or gradually upregulated from days 1 to 10 in both root and variants, and identify splicing isoforms of known genes shoot under −P stress, after which the expression remained (Wang et al. 2009a). More importantly, it could accurately stable for up to 30 days (Fig. 2a). After transferring the quantify gene expression levels over a broad dynamic range plants to growth conditions with normal Pi concentration and detect transcripts expressed at either very low or very (+P) to allow recovery from stress, the expression was high levels including subtle changes which could not be immediately downregulated within 1 day. Under ++P stress characterized by microarray-based approaches (Li et al. however, OsIPS1 was not upregulated in both root and 2009; He et al. 2010; Lu et al. 2010; Mizuno et al. 2010). shoot (Fig. 2b). This gene was also not upregulated under We therefore used the massive parallel sequencing technol- starvation with other essential nutrients such as nitrogen, ogy by mRNA-Seq to elucidate the rice transcriptome calcium and magnesium (data not shown). Based on these under stress due to –P and ++P in order to provide a results, we focused our gene expression profiling analysis 52 Rice (2011) 4:50–65 Fig. 1 Effect of Pi stress on rice growth. a Phenotypic changes in rice represent the mean±SE for three replicates for each treatment. c, d plants after 30 days of growth in culture medium with normal Pi Changes in dry weight and shoot/root weight ratio (in mg) of rice concentration (+P, control), Pi starvation (−P), and over-abundant Pi seedlings +P, −P, and ++P treatment conditions. The values represent supply (++P). Both the shoot and root showed growth retardation under the mean±SE for three replicates for each treatment. Changes in dry Pi stress. b Total phosphorus content in roots and shoots after weight and shoot/root weight ratio became prominent at around 30 days in +P, −P, and ++P treatment conditions. The values 10 days after Pi stress treatment. by mRNA-Seq during the early growth stages particularly ces were produced and used for mapping onto the reference at 1 (early), 5 (middle), and 10 days (late) after −P and ++P Nipponbare genome sequence (Table 1). Of the 6 to 30 stress treatments. Additionally, the expression profiles at million quality-evaluated reads (Passed Filtered reads) from recovery after 10 days of −P and ++P stress treatments were each library, 12.3% to 39.8% were mapped to single also investigated (Supplementary Fig. S1). locations (unique) in the genome whereas 23.7% to 45.7% were mapped to multiple locations in the genome. The Sequencing and short-read mapping expression level of all unique transcripts mapped onto the genome were quantified as reads per kilobase of exon A total of 16 libraries were used for mRNA-Seq analysis model per million mapped (RPKM) values (Mortazavi et al. using the Illumina Genome Analyzer IIx (Illumina Inc., San 2008). On the other hand, 19.0% to 51.1% of the sequence Diego, CA, USA). Overall, 195 million short-read sequen- reads from each library had no match in the genome. These Rice (2011) 4:50–65 53 Fig. 2 OsIPS1 expression under Pi stress treatments. The expression of 30 days after stress and at 1 or 2 days after recovery from stress (Re1 and OsIPS1 in root and shoot under Pi starvation (a) and over-abundant Pi Re2). The amplification of β-tubulin in RT-PCR was used as control. supply (b) was analyzed by RT-PCR at various time points from 0 to unmapped reads may include low-quality reads, sequencing that 95% of the transcripts in the 44 K microarray platform errors, or sequences derived from adaptors and contami- could be validated by mRNA-Seq. A total of 11,888 nating organisms (Mizuno et al. 2010). transcripts from the root and 11,098 transcripts from the shoot were not represented in the 44 K array but showed mRNA-Seq vs. microarray analysis complete match with the annotated RAP transcripts. These transcripts correspond to gene models without EST and/or We correlated the mRNA-Seq data with the gene expres- full-length cDNA support as well as those predicted based sion profiles derived from analysis using the rice 44 K on gene prediction programs. The mRNA-Seq approach oligomicroarray platform (Agilent Technologies, Palo Alto, could therefore provide a more comprehensive analysis of CA, USA). First, we re-mapped the microarray probe the transcriptome including about half of RAP transcripts sequences into the updated rice genome assembly (Build not represented in the microarray platform. We used the 5.0 pseudomolecules) to make a valid comparison of the Cufflinks program for comparative assembly and estimation microarray data and the mRNA-Seq data. As a result, we of abundance of transcripts in each sample (Trapnell et al. confirmed that 30,026 RAP transcripts (23,113 loci) were 2010). By excluding the annotated transcripts in RAP and represented as probes in the 44 K array (Fig. 3). Sequence the MSU rice gene models (http://rice.plantbiology.msu. comparison with the mRNA-Seq data showed that 28,680 edu/), we were able to detect a total of 8,590 transcripts transcripts from the root and 28,650 transcripts from the from the root and 8,193 transcripts from the shoot which shoot matched with the microarray probes. This indicates have not yet been previously annotated. Table 1 Mapping of mRNA-Seq reads obtained from each sample into the reference rice genome sequence mRNA-Seq library Total PF reads Unique % Multiple % Unmapped % Root, +P_0d (control) 13,922,030 3,738,255 26.9 4,344,200 31.2 5,839,575 41.9 Root, −P_1d 8,411,989 2,314,238 27.5 2,827,342 33.6 3,270,409 38.9 Root, −P_5d 30,403,973 8,057,787 26.5 7,213,838 23.7 15,132,348 49.8 Root, −P_10d 6,696,463 2,495,317 37.3 2,459,102 36.7 1,742,044 26.0 Root, −P_10d_Re_1d 8,528,677 3,311,738 38.8 2,857,712 33.5 2,359,227 27.7 Root, ++P,_1d 15,694,862 4,738,300 30.2 5,013,389 31.9 5,943,173 37.9 Root, ++P_5d 13,586,681 2,642,094 19.4 4,000,708 29.4 6,943,879 51.1 Root, ++P_10d 11,027,653 1,353,233 12.3 4,597,218 41.7 5,077,202 46.0 Root, ++P_10d_Re_1d 16,771,413 3,827,077 22.8 5,270,793 31.4 7,673,543 45.8 Shoot, +P_0d (control) 8,515,740 2,205,725 25.9 3,181,292 37.4 3,128,723 36.7 Shoot, −P_1d 8,787,606 2,380,054 27.1 3,319,901 37.8 3,087,651 35.1 Shoot, −P_5d 9,281,244 2,278,998 24.6 4,240,468 45.7 2,761,778 29.8 Shoot, −P_10d 8,258,753 2,943,631 35.6 3,744,753 45.3 1,570,369 19.0 Shoot, ++P_1d 12,145,854 3,090,635 25.4 4,167,569 34.3 4,887,650 40.2 Shoot, ++P_5d 14,586,873 4,126,432 28.3 6,302,429 43.2 4,158,012 28.5 Shoot, ++P_10d 9,086,827 3,612,159 39.8 3,601,872 39.6 1,872,796 20.6 54 Rice (2011) 4:50–65 Fig. 3 Transcripts identified in RAP annotations and 44 K microarray probes. Unique tran- scripts mapped to the genome were classified as follows: (1) transcripts with match in RAD- DB annotation and microarray, (2) transcripts with match in RAP-DB annotation but no match in microarray, and (3) unannotated transcripts identi- fied by the Cufflinks program. In both root and shoot samples, about 95% of transcripts were supported by the microarray probes. About 50% of tran- scripts not identified by the array matched with the RAP-DB transcripts. A total of 10,388 unique unannotated transcripts were identified by Cufflinks program. Gene expression profiling the 44 K microarray platform the expression profiles for known Pi-related genes can be was performed in root samples at 1, 5, and 10 days after –P validated by mRNA-Seq. and ++P stress treatments. The signal intensities of 28,686 transcripts represented in the microarray were plotted Characterization of Pi stress-responsive transcripts against the RPKM values of corresponding transcripts (Supplementary Fig. S2). For overall gene expression, we We used the G test (FDR, <0.01) on the RPKM-derived observed an average correlation coefficient of >0.8, read counts to determine the differences in gene expression suggesting a clear validation of the microarray-based gene at various time points during −P and ++P treatments. A expression profiling data with the mRNA-Seq data. transcript is considered responsive if the expression level after recovery from stress (−P_Re_1d and ++P_Re_1d) Transcripts corresponding to known Pi-related genes reverted to the same level before treatment or the control condition. As a result, we were able to classify the We were able to identify transcripts with match to expression patterns into 24 expression patterns based on previously reported Pi-related genes including those in- the statistical significance between treatment conditions (−P volved in sumoylation (SIZ1) and signaling (OsIPS1 and and ++P) and the expression levels in root and shoot at each OsIPS2), transcription factors (phosphate starvation re- stage of stress designated as early for 1 day, middle for sponse (PHR)1 and PHR2), phosphatases (PAP2 and 5 days, and late for 10 days after stress treatment (Table 2). PAP10), etc. (Supplementary Table S1). As expected, Under stress recovery conditions, the responsive transcripts analysis using Cufflinks also revealed the abundance and which were upregulated in roots under −P stress showed differential expression of these transcripts. The expression similar downregulation as OsIPS1 after recovery from levels of transcripts from root and shoot samples at various stress starvation. The same pattern of fast and marked Pi starvation (−P_1d, −P_5d, and −P_10d) as well as over- recovery of expression after re-addition of Pi was also abundant Pi supply (++P_1d, ++P_5d, and ++P_10d) shown in Arabidopsis with >20-fold increase in AtPht1;4 treatments based on RPKM values indicate significant and AtPAP2 signal intensity as compared with the level expression of these Pi-related genes (Supplementary during P deprivation (Morcuende et al. 2007). Similarly, Table S2). Moreover, comparison of the fold change in upregulated transcripts under ++P stress conditions were RPKM values from control and treatment conditions with also downregulated after 1-day recovery from stress due to the normalized signal intensity clearly showed a similar over-abundant Pi supply. pattern of expression profile obtained by mRNA-Seq and The highest number of responsive (upregulated and microarray-based analysis in all treatments of –P and ++P downregulated) transcripts was observed in root at late stress (Supplementary Table S3). These results indicate that stage, which was almost twice the number in root at early Rice (2011) 4:50–65 55 Table 2 Number of Pi stress-responsive transcripts identified by G test Expression patterns Number of responsive transcripts Tissue Treatment Response Time point Total RAP/array RAP not Unannotated Novel Transcripts identified (%) with P1BS by array Root −P Upregulated Early (1 d) 230 201 27 2 12.6 43 Middle (5 d) 354 286 68 0 19.2 71 Late (10 d) 503 434 68 1 13.7 170 Dowregulated Early (1 d) 148 112 35 1 24.3 21 Middle (5 d) 171 156 15 0 8.8 31 Late (10 d) 246 221 24 1 10.2 54 ++P Upregulated Early (1 d) 233 192 40 1 17.6 58 Middle (5 d) 472 399 69 4 15.5 105 Late (10 d) 265 203 57 5 23.4 54 Dowregulated Early (1 d) 209 186 21 2 11.0 39 Middle (5 d) 203 181 22 0 10.8 39 Late (10 d) 208 186 22 0 10.6 52 Shoot −P Upregulated Early (1 d) 224 199 22 3 11.2 47 Middle (5 d) 117 103 14 0 12.0 25 Late (10 d) 216 182 32 2 15.7 74 Dowregulated Early (1 d) 454 400 53 1 11.9 110 Middle (5 d) 83 71 11 1 14.5 20 Late (10 d) 205 183 22 0 10.7 39 ++P Upregulated Early (1 d) 588 498 88 2 15.3 145 Middle (5 d) 466 397 68 1 14.8 94 Late (10 d) 166 139 27 0 16.3 40 Dowregulated Early (1 d) 1,136 1,014 118 4 10.7 252 Middle (5 d) 168 152 15 1 9.5 37 Late (10 d) 82 77 5 0 6.1 15 P1BS is significantly overrepresented in the 1-kb upstream region at 19.6% background level (responsive transcripts/non-responsive transcripts= 7,164/36,588) and FDR, <0.001 in Fisher’s exact test stage of –P stress (Table 2; Fig. 4). The highest number of 13.6%) of all responsive transcripts were classified as novel responsive transcripts was found in shoot at early stage and transcripts in each group. in root at middle stage of ++P stress. Downregulation of The Pi stress-responsive genes were also analyzed based transcription in shoot at early ++P treatment was most on the expression of PHR1. This gene is known to common among the 24 expression patterns. The expression upregulate −P stress-related genes such as high-affinity Pi patterns of representative transcripts such as WRKY transporters, RNases, and acid phosphatases by binding to (Os03t0321700), Pht1;2 (Os03t0150800), and OsIPS1 an imperfect palindromic sequence (P1BS, GNATATNC) in (Os03t0146800) clearly show upregulation at early, middle, the promoter regions (Rubio et al. 2001). The presence of and late stage of −P stress, respectively (Supplementary P1BS cis-acting element in the 1-kb upstream region of Pi Fig. S3). stress-responsive transcripts was therefore used as a The responsive transcripts identified based on G test also measure of response under stress. Using the signal scan included approximately 900 RAP transcripts which were search in the PLACE database (http://www.dna.affrc.go.jp/ not represented in the microarray as well as 28 unannotated PLACE/), we investigated the 1-kb upstream region of Pi transcripts revealed by the Cufflinks program (Table 2). stress-responsive representative RAP transcripts. The P1BS These unsupported RAP transcripts and unannotated tran- cis-acting elements were significantly overrepresented in scripts were classified as novel transcripts induced during – the 1-kb upstream region of transcripts upregulated in root P and ++P stress. Approximately 6.1% to 24.3% (average and shoot at late –P stress (Table 2). 56 Rice (2011) 4:50–65 and late) based on the different expression patterns (Supple- mentary Fig. S4). The GO categories showed specific trends for many Pi-related genes in response to –Pstress (Supple- mentary Fig. S5) and ++P stress (Supplementary Fig. S6). 1. Regulatory proteins/riboregulator: during –P stress, we observed specific pattern of response among genes for regulatory proteins/riboregulators such as OsIPS1, WRKY, SPX (SYG/PHO81/XPR1) domain-containing genes, histone, etc. OsIPS1 and OsIPS2 were strongly upregulated in root and shoot but PHO2/UBC24 was downregulated in shoot during late period of –P stress. The expression of OsIPS2 and PHO2/UBC24, which were not supported by the microarray, were confirmed here. The expression patterns of several transcription factors and signaling molecules under −P stress were clarified. In particular, WRKY (Os03t0321700) was upregulated in root during early −P stress. In Arabi- dopsis, AtWRKY6 and AtWRKY42 are known to modulate PHO1 transcription (Chen et al. 2009) whereas AtWRKY75 modulates Pi acquisition and root development (Devaiah et al. 2007). NAM (NAC) transcription factor (Os03t0624600) was downregu- lated in root at early stage of −P stress. The NAC gene family shows diverse functions in both plant develop- ment and stress responses as in the AtNAC1-mediated auxin signaling to promote lateral root development (Xie et al. 2000). Several SPX (SYG/PHO81/XPR1) domain-containing genes involved in response to environmental cues or internal regulation of nutrition homeostasis such as SPX1, SPX2, and SPX6 were upregulated only in shoot during late –P stress. In rice, SPX1 is partly involved in regulating the induction of some Pi transporters (Wang et al. 2009b). Histone H4 (Os04t0583600) and histone H2B (Os05t0574300) were upregulated in shoot during early –Pstress whereas histone H1 (Os03t0799000) and histone H2B Fig. 4 Distribution of upregulated and downregulated transcripts in (Os01t0152900) were downregulated in shoot during response to Pi stress. The total number of upregulated (upper)or late –P stress thereby supporting the chromatin-level downregulated (lower) transcripts identified by mRNA-Seq were regulation of response to Pi starvation (Smith et al. determined at various time points of stress (early, middle, and late) 2010). Under ++P stress treatments, OsIPS1 and during Pi starvation (−P) or phosphate over-abundant Pi supply (++P) treatments. Each bar shows the distribution of transcripts with match OsIPS2 did not show any specific pattern of expression in RAP-DB annotation and microarray (gray), transcripts with match in either root or shoot. Downregulation of auxin- in RAP-DB annotation but no match in microarray (white), and regulated transcripts such as Os03t0742900 in root at unannotated transcripts (black). early stage, Os05t0230700 in shoot at middle stage, and Os04t0519700 (auxin response factor) in root at Functional characterization of transcripts generated late stage may suggest specific functions of these genes in root and shoot growth under ++P stress. We observed by mRNA-Seq the downregulation of PTF1 in shoot at early ++P stress. Overexpression of this gene however, could enhance We used the Gene Ontology (GO) classification to estimate the functional categories of upregulated and downregulated tolerance to –P in transgenic rice (Yi et al. 2005). Among histone genes, ++P stress induced the upregula- genes in the root and shoot under various stress treatments (−P and ++P) and at different stages of stress (early, middle, tion of histone-like transcription factor (Os03t0413000) Rice (2011) 4:50–65 57 in shoot at early stage and histone H3 (Os06t0130900) 4. Metabolism (lipid): we observed specific pattern of in shoot at middle stage, and downregulation of histone response to Pi stress among transcripts related to lipid H2B (Os01t0152900) in root at late stage. and fatty-acid metabolism. During the late stage of –P 2. Metabolism (energy): many transcripts categorized as stress, there was prominent upregulation of lipid tricarboxylic acid (TCA) cycle components were induced metabolism-related genes such as UDP-sulfoquinovose during –P stress in both root and shoot, and most were synthase 1 (SQD1), SQD2, monogalactosyldiacylglycerol strongly upregulated from the middle to the late stage of – (MGDG) synthases, DGD1 (digalactosyl diacylglycerol P stress. The gene encoding isocitrate dehydrogenase 1), DGD2, and glycerophosphoryl diester phosphodies- (ICDH) was downregulated in root at the late stage of –P terase (GDPD), and transcripts encoding lipid transfer stress. As a gene involved in the degradation of citrate in proteins (Os07t0174400 and Os03t0794000). In particu- the TCA cycle, the repression of ICDH could induce an lar, SQD1, SQD2,and MGDG were upregulated in both increase in the internal citrate concentration (Delhaize et root and shoot during −P stress. Lipid and fatty-acid al. 2003). The gene encoding malate dehydrogenase metabolism-related genes involved in the synthesis of (MDH) (Os08t0434300) was upregulated in root at late galacto- and sulfolipids were strongly induced by −P stage of –P stress. In alfalfa, overexpression of MDH treatment. Lipase encoding transcripts such as caused an increase in the exudation of various organic Os01t0215000 and Os01t0710700 were upregulated in acids and subsequently resulted in increased P accumu- root at early ++P stress. Induction of these lipid lation in acid soils (Tesfaye et al. 2001). Among the OsA metabolism-related transcripts maybe associated with genes, the plasma membrane H -ATPase encoding OsA7 changes in the lipid composition and fluidity of the was strongly expressed under normal Pi concentration membranes that control metabolic and cellular processes. but the expression gradually decreased during –Pstress 5. Pi transport: among genes encoding high-affinity Pi (Supplementary Table S2), a response that maybe transporters, Pht1;2 and Pht1;8 were upregulatedinthe associated to nutrient uptake and translation during Pi roots in the middle stage whereas Pht1;3, Pht1;4, Pht1;6, starvation. During ++P stress, OsA2 was upregulated in Pht1;9,and Pht1;10 were upregulatedinthe rootsinthe root and shoot at early stage, OsA3 was upregulated in late stage of –P stress. Although high Pi concentration is shoot at middle stage, and OsA7 was downregulated in toxic to plants, the Pi transporter encoding genes were root at middle stage. Genes encoding glycolysis enzymes mostly upregulated during ++P stress including Pht1;4 were mostly upregulated as in the case of fructose- and Pht2;1 in shoot at early stage, Pht1;8 in shoot at late bisphosphate aldolase (Os11t0171300) in shoot at early stage, and Pht1;1 in root at middle stage. Only Pht1;8 stage, phosphofructokinase (Os05t0524400) in shoot at was downregulated in the root at middle stage of ++P middle stage, Os10t0405600 in shoot at late stage, and stress. These expression patterns may suggest specific enolase (Os06t0136600) in shoot at middle stage of ++P functions of these genes in Pi uptake and homeostasis stress. Enzymes involved in TCA cycle and glycolysis under conditions of over-abundant Pi supply. Among the may function to produce organic acids, which would transcripts for Pi transporters, only Pht1;11 was not permit the recycling of P from phosphorylated inter- detected (Supplementary Table S2) probably because it is mediates. Organic acids are also known to help release Pi specifically activated during mycorrhizal symbiosis from organic or insoluble inorganic Pi compounds (Paszkowski et al. 2002). outside the plant. 6. Pi remobilization: many OsRNS genes showed specific 3. Metabolism (carbohydrate): many genes encoding su- expression patterns during –P and ++P stress. It has crose metabolism-related proteins showed specific re- been suggested that RNS may have redundant functions sponse in root and shoot during –P and ++P stress. This or specific functions in different biological processes may be associated with the tightly controlled mechanisms and tissues under different biotic and abiotic stresses that allow the coordination of Pi homeostasis with carbon (MacIntosh et al. 2010). Purple acid phosphatase 10 status and photosynthesis (Wissuwa et al. 2005). Sugars (PAP10) and acid phosphatase 1 (ACP1) were upregu- generated in the shoots and transported through the lated in both root and shoot at late –P stress, which may phloem are involved in establishing a physiological, suggest a possible role in Pi acquisition and metabolism biochemical, and molecular response to −Pstress in during –P stress. plants (Hammond and White 2008). Sucrose phospha- 7. Morphological changes: From the middle to late stage tase encoding transcript (Os01t0376700) was upregu- of –P stress, genes encoding cell wall-related proteins lated in root at late –P stress whereas sucrose synthase such as pectin methylesterase 8 (Os01t0311800) and encoding transcripts (Os06t0194900 and Os03t0401300) cellulose synthase (Os07t0208500) were upregulated were downregulated in shoot and root at late–Pstress suggesting possible roles in the control of cell wall synthesis and extensibility during stress. Upregulation and root at middle ++P stress. 58 Rice (2011) 4:50–65 during the initial response to stress suggests that plant sequence (P1BS) in the 1-kb upstream or downstream growth was affected slightly during the first 10 days of – region (Supplementary Table S5). Among10,388 unanno- P stress. An auxin response factor (Os04t0519700) and tated transcripts, approximately one third (3,006 tran- an Aux/IAA binding protein (Os03t0742900) may be scripts) have amino acid sequence homology in the involved in regulating the root architecture in response to RefSeq and UniProt databases according to BLASTX −P stress by promoting lateral root elongation and search (identity, >30%; coverage, >30%) (Supplemental changing auxin distribution within root cells (Nacry et Table S6). Two unannotated Pi stress-responsive tran- al. 2005). Upregulation of anthocyanin metabolism- scripts showed homology; NP_CUFF.161528.1 to an related gene such as phenylalanine ammonia-lyase unknown protein and NP_CUFF.2572.1 to a leucine-rich (Os12t0520200) in root at late stage of –P, may be repeat family protein/protein kinase family protein. associated with its function in protecting cells from Validation experiments by qRT-PCR analyses confirmed the the stress. Many genes encoding auxin-related expression of these unannotated transcripts under −P or ++P regulatory proteins, cell wall-related proteins and stress (Fig. 5). NP_CUFF.26799.1 was upregulated in both cytoskeleton proteins were downregulated from the root and shoot at late −P stress and immediately down- early to late stage without any morphological changes for regulated after recovery from stress. NP_CUFF.141685.1 was the first 10 days during ++P stress. The genes encoding upregulated in root at late ++P stress and immediately protein synthesis-related proteins such as aminoacyl- downregulated after recovery from stress. Some of these tRNA synthetase (Os03t0749300, Os05t0150900) trans- transcripts may be regulated through the P1BS in their lation initiation factor 2 (Os03t0296400) and translation flanking regions, although they did not appear to have high elongation factor EF1B (Os06t0571400) were down- homology to one another. regulated, suggesting a repression of protein turnover and recycling of amino acids in roots and shoots. These results also show that ++P stress is a more severe form of Discussion stress for rice growth than −Pstress. 8. Homeostasis of other ions: Genes encoding transport Identification of responsive non-protein coding transcripts proteins such as potassium transporter (Os09t0448200), a under Pi stress heavy metal transporter (Os02t0530100) and an iron– sulfur cluster assembly related protein such as NifU-like We used the mRNA-seq approach to characterize the protein (Os01t0662600) showed specific expression transcriptome of rice under Pi stress at high resolution. patterns that may be associated to their functions in This is the first report of a whole transcriptome analysis involving stress due to Pi starvation and over-abundant Pi adjusting the nutrient balance under –P stress. As Pi is important in numerous energy-requiring metabolic and supply. In addition to the confirmation of almost 77% of transport processes, genes encoding transport proteins RAP transcripts, the results also provide evidence for such as heavy metal transporter (Os07t0258400, expression of transcripts not represented in the microarray Os01t0125600, and Os01t0976300) were upregulated including transcripts without FL-cDNA support, transcripts during ++P stress. based on gene prediction programs, or those identified by homology to other plants (Fig. 3). Furthermore, the Pi stress-responsive transcripts can be functionally character- Characterization of unannotated transcripts induced by −P ized based on the expression pattern in root and shoot at and ++P stress early, middle, or late stage of −P or ++P stress. A total of 10,388 transcripts identified using the Cufflinks program The unannotated Pi-responsive transcripts validated by G test have not yet been previously annotated. Some of these a ranged in length from 135 to 1,750 bp with an average of unannotated transcripts may have been induced by alterna- 600 bp (Supplementary Table S4). These transcripts were tive promoters as they have been identified in the promoter distributed in the 12 chromosomes. The nearest adjacent region of RAP transcripts, within the 5′ and 3′ exons transcripts were located from a distance of 1 to more than supported by ESTs or FL-cDNAs, or within the extended 3,000 bp. Most transcripts showed specific expression regions of the exons. Although unannotated transcripts may pattern in root and shoot at early, middle or late stage also be found in the introns, we could not differentiate of −P and ++P stress. Two unannotated transcripts, intron-derived transcripts from newly formed exon tran- namely, NP_CUFF.172640.1 and NP_CUFF.136416.1, were script and the mRNA-seq approach could not define found to contain an intron which showed their direction based whether the RNA is sense or antisense. In order to on the junction sequence. Analysis of flanking regions also characterize intron-derived transcripts, it maybe necessary revealed that most of these transcripts contain a PHR1-binding to obtain longer reads and pair-end sequence data. Rice (2011) 4:50–65 59 NP_CUFF.111084.1 1200 ++P ++P ++P * * 600 20 2 10 0 0 0 15 10 1015 10 10 15 10 1015 10 10 15 10 1015 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day NP_CUFF.88762.1: Root_++P_Early_Down NP_CUFF.141685.1: Root_++P_Late_Up 2.5 ++P ++P ++P 0.8 0.6 1.5 0.4 0.2 0.5 * 0 0 0 15 10 1015 10 10 15 10 15 10 15 10 15 10 10 10 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day ++P ++P ++P 25 * * 20 10 * 0 0 15 10 1015 10 10 15 10 1015 10 10 15 10 1015 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day NP_CUFF.2572.1: Shoot_++P_Early_Down NP_CUFF.76959.1: Shoot_++P_Early_Down 1.4 0.8 1.2 ++P ++P ++P 0.6 0.8 0.4 0.4 0.2 0.2 15 10 15 10 15 10 1015 10 10 10 10 15 10 1015 10 10 Re1 Re1 Re1 Re1 Re1 Re1 day day day Fig. 5 Validation of responsive unannotated transcripts. The expression expression on day 10 during –P/++P and at 1 day after recovery from of OsIPS1 and several unannotated transcripts in response to –Pand ++P stress (Re1). The bars represent the mean±SE level of transcripts from was validated by qRT-PCR analysis. OsIPS1 is a –P upregulation index three experiments. The asterisk indicates the responsive point based on G gene with statistical significance in the roots (upper)and shoots (lower) test (FDR, <0.01). The primers used in qRT-PCR analysis for as determined by G test. Transcript expression levels were normalized unannotated transcripts are described in Supplementary Table S5. using an internal control (Ubiqutin1) and plotted relative to the We have identified a total of 10,388 unique unannotated various Pi stress treatments (Supplementary Table S4). transcripts from the root and shoot with 3,006 transcripts These transcripts comprised less than 1% of the responsive showing homology to various proteins (Supplementary transcripts and had much lower expression levels than Table S6). Among them, 28 transcripts have been found RAP-representative transcripts for all conditions. For to be rice-specific and 303 transcripts showed homology to example, the average RPKM values (±SD) for unannotated other known plant genes. The remaining unannotated and RAP-representative transcripts in the root control transcripts (7,382) without homology to any protein may sample were 0.8 (±2.6) and 23.8 (±94.5), respectively. As include non-protein coding transcripts, novel protein tran- these transcripts also showed small differential expression scripts, rare transcripts that are expressed at low copies, between conditions that could not be detected by statistical transcripts with very low expression levels, or even tests, it maybe difficult to characterize all these transcripts. transcripts which may have lethal functions in Escherichia Two unannotated transcripts showed protein homology to coli. We particularly focused on the 28 unannotated known proteins, one to leucine family gene and another to transcripts, which were differentially expressed under an unknown protein. Furthermore, the presence of P1BS Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression Relative Expression 60 Rice (2011) 4:50–65 cis-acting elements was confirmed in the promoter regions gene expression between the root and shoot (Supplemen- of ten unannotated transcripts. These findings establish the tary Fig. S7). These morphological and physiological utility of mRNA-Seq in identifying uncharacterized tran- changes are associated with the changes in gene expression. scripts that are expressed in response to Pi stress. The differences in gene expression between −P and ++P The expression of approximately 70% of the responsive stress in root and shoot (Supplementary Figure S7a, b)as unannotated transcripts was confirmed by qRT-PCR anal- well as the differences in gene expression between root and ysis. Several transcripts (NP_CUFF.26799.1 and shoot under −Pand ++Pstress (Supplementary Figure S7c, NP_CUFF.141685.1) showed a sharp and reversible re- d) suggest that Pi stress treatments have a significant effect sponse to −P or ++P stress (Fig. 5). Most of these on the overall expression profile of many genes involved in transcripts do not contain any ORFs and may function as plant response to stress. We have confirmed the response of non-coding transcripts similar to OsIPS1. Some responsive many previously reported phosphate-related genes. Although unannotated transcripts could not be validated because of some of the genes that showed significant changes in very low expression levels that could not be detected by expression may also be involved in the developmental stage, qRT-PCR. Overall, these validation experiments further there is no doubt that most of the responsive genes were suggest that mRNA-Seq can provide strong evidence of expressed due to Pi stress as shown by the statistical transcript expression. It has also been efficiently used in confirmation of the differences in gene expression over time confirming the function of miRNAs such as miR399 in during −P and ++P treatments. After recovery from stress, response to Pi stress. In Arabidopsis, the interaction we also found that many genes upregulated or down- between the Pi- and nitrate-limited signaling pathways regulated genes during Pi stress, including OsIPS1,reverted affecting E3 ligase gene NLA and anthocyanin synthesis to the same level of expression before stress treatments. was found to be mediated by miR827, which targets the Among the –P responsive transcripts, 316 transcripts were transcripts of specific proteins that contain the SPX domain commonly expressed in both root and shoot. These transcripts (Pant et al. 2009). The rice homolog OsmiR827 exhibits likely function in basic –P signaling cascade that impacts Pi different preferences in tissue expression and regulates transport and Pi remobilization in plants, and included different classes of genes sharing a common SPX domain upregulated transcripts such as OsIPS1, OsIPS2, SQD1, with AtmiR827 (Lin et al. 2010). The existence of miRNA- OsRNS3,and PAP10. On the other hand, PHO2/UBC24, dependent pathway in distinct aspects of Pi homeostasis in which regulates Pi allocation between roots and shoots, was these two species may suggest a species-specific pathway downregulated in the shoots. A total of 1,336 –Presponsive where many unannotated transcripts may play significant transcripts in the roots were not responsive to in the shoots roles. It has been reported that miR2111 may target the E3 (Supplementary Fig. S7c). These transcripts likely function ligase gene and a calcineurin-like phosphoesterase gene in the cell-rescue system that impacts Pi uptake, carbon- under −P stress (Hsieh et al. 2009). Additionally, miRNAs assimilation reduction and root-system morphology by such as miR169, miR395, and miR398, which respond to enhancing root growth, and included upregulated transcripts various nutritional stresses, were also downregulated in that encode many phosphate transporters (Pht1;2, Pht1;3; response to Pi stress (Hsieh et al. 2009). Therefore, some of Pht1;6, Pht1;8, Pht1;9, and Pht1;10), an enzyme in the TCA the unannotated transcripts detected here including non- cycle (PEPC/Os08t0366000) and an auxin response factor protein coding transcripts and small RNAs may have (Os04t0519700). It has been reported that Pht1;8 and Pht1;9 important functions in the signaling of plant root responses were located downstream of the PHO2/UBC24-dependent to changing Pi concentrations and the coordination of signaling pathway (Bari et al. 2006). homeostatic pathways of Pi and other nutrients. A comparison of the ++P responsive transcripts in root and shoot showed that 358 transcripts were commonly expressed. Responsive genes related to morphological However, many transcripts with different functions responded and physiological changes in root and shoot under ++P stress. In particular, PHO2/ UBC24 was upregulated in root at early ++P stress. This Here, we found that –P stress inriceisassociatedwithseveral expression may be sufficient to repress the –P upregulated morphological changes such as brown coloration, roughening transcripts such as Pht1;8, which is downstream of the of the roots and a decreased ratio of shoot-to-root weight PHO2/UBC24-dependent signaling pathway. Pht1;1 and (Fig. 1). Under ++P stress, we also observed weak roots and Pht1;4 were also upregulated in root during ++P, which as shoots, and an increased ratio of shoot-to-root weight with suggested by Bari et al. (2006) may not be a part of the increased Pi uptake. Interestingly, although the ratio of PHO2/UBC24-dependent signaling pathway. There were shoot-to-root weight clearly reversed between 10 and 15 days 2,248 ++P responsive transcripts in shoot that were not in the –P and ++P stress-treated plants, the weights of the responsive in root. These included upregulated genes such as plants decreased concordantly as reflected in differences in Pht1;8 and Pht2;1 which may contribute to the establish- Rice (2011) 4:50–65 61 ment of high Pi in shoot. Pht1;8 may also be upregulated by starvation-induced genes (Bustos et al. 2010). Therefore, other factors without downregulating PHO2/UBC24 tran- PHR1 is maybe upregulated in late stages of −P stress via a scription. Protein synthesis (Os01t0678600) and fatty-acid P1BS cis-acting element. Since the rest of the genes metabolism (Os02t0205500) transcripts were downregulated. upregulated in late stage of −P stress do not have a P1BS A large number of responsive transcripts detected in shoot at cis-acting element in the promoter regions, other unknown early ++P stress may have diverse functions in Pi metabo- cis-acting elements may also be involved in upregulation in lism, photosynthesis and ATP biosynthesis, and may be late stage of −P stress as well as in other expressions transiently activated for plant growth. At the same time, high patterns in response to Pi stress. Some transcripts may be Pi concentrations may inhibit plant growth because the genes regulated through the release from repression of PHO2/ encoding many nuclear proteins, such as poly(A) polymerase UBC24 or regulated by other factors such as PHL1. Both (e.g., Os06t0558700), many protein synthesis-related pro- PHR1 and PHL1 were found to be partially redundant and teins, such as translation initiation factor 2 (e.g., have a central role in the control of physiological and Os03t0296400) and aminoacyl-tRNA synthetases (e.g., molecular responses to –P stress (Bustos et al. 2010). Os06t0645400) were downregulated in shoot. The plant We searched gene loci with two or more alternative may then acclimate and use less energy because the number splicing isoforms supported by FL-cDNAs to determine of responsive transcripts decreases from the middle to the whether alternative splicing mechanisms affect gene ex- late stage of stress in shoot. Interestingly, Pht1;4 was pression during Pi stress. Among 34,780 representative loci upregulated in root and shoot under both –Pand ++P stress. on the rice genome (IRGSP Build 5 Pseudomolecule), This finding suggests that Pht1;4 can easilyadjusttothe alternative splicing isoforms were confirmed in 5,625 cellular Pi condition. In Arabidopsis, it has been reported (16.2%) representative loci. To examine changes in the that AtPHT1;4 was induced by sucrose and is thought to be usage of alternative splicing isoforms against Pi stress, the upstream of the hexokinase sugar-sensing pathways (Lejay et isoforms were compared in terms of Pi stress-responsive al. 2008). Thus, morphological changes for adaptation in rice expression patterns. As a result, the alternative splicing due to changes in Pi concentration may be controlled by isoforms of 194 loci (0.6%) showed different expression local and systemic gene expression. patterns. This suggests that the Pi stress-responsive expres- The differences in total phosphorus content of rice plants sion is also regulated by alternative splicing mechanisms grown under normal Pi concentration (+P), Pi starvation (−P), with respect to the timing of the response or tissue and over-abundant Pi supply (++P) are direct manifestations of specificity. Further analysis of the upstream regions of Pi transporter activities. We identified several responsive responsive transcripts and accompanying changes in alter- transporter genes such as Pht1;4, Pht1;8, and Pht2;1. Genetic native splicing may reveal expression networks underlying modification of Pi transport proteins are maybe the key in Pi homeostasis in plants. developing tolerance to both –P and ++P stress. Britto and Kronzucker (2008) reported that the high-affinity transport A genome-wide transcriptional analysis of Pi stress system is strongly downregulated under K -replete conditions and is conversely strongly upregulated under K -starvation This study has shown that RNA-Seq can be used to determine conditions and that the LATS (low-affinity transport system)- transcript expression levels more accurately and capture the + + range K influx may be increased generally by increased K transcriptome dynamics associated with Pi stress across (Britto and Kronzucker. 2008). We also observed upregula- different time points and tissues. Interestingly, SIZ1, PHR1, tion of several low- and high-affinity transporter genes such and PHR2, which are located upstream of OsIPS1,were as Pht1;1, Pht1;8, and Pht2;1 under ++P stress. Therefore, a slightly downregulated after 10 days under −Pstress in roots portion of the Pi influx may be regulated by a transport and shoots and slightly upregulated at 1 d after recovery in system similartothat ofK influx. roots, an expression pattern which was in contrast with that of OsIPS1. These genes were previously reported to be Promoter analysis and alternative splicing isoforms upregulated slightly under −P stress based on qRT-PCR of responsive genes analysis and northern analysis (Rubio et al. 2001;Miura et al. 2005). This apparent discrepancy indicates that a more The flanking regions of most Pi stress-responsive tran- detailed analysis of gene expression during different stages scripts contain a P1BS cis-acting element which were of stress via mRNA-Seq could be useful in detecting subtle significantly over represented in the 1-kb upstream regions changes in gene expression. of 170 transcripts (33.8%) from the roots and 74 transcripts Using mRNA-Seq, we were able to identify many novel (34.3%) from the shoots that were upregulated at the late transcripts which could not be revealed by microarray- stage of −P stress. In Arabidopsis, the direct targets of based technology. We were also able to characterize even subtle changes in their expression at different stages of PHR1 are greatly enriched in P1BS-containing Pi 62 Rice (2011) 4:50–65 stress under –P and ++P treatments. These results will be (Life Technologies, Rockville, MD, USA). The RNA was useful in a detailed functional analysis of an accurate gene treated with DNase I (Takara, Shiga, Japan), and the first- set provided in the high-quality map-based genome strand cDNA was synthesized using the SuperScript® III sequence of the model rice cultivar Nipponbare. Towards first-strand synthesis system (Invitrogen, Carlsbad, CA, this goal, we also constructed transcriptome viewer in USA) according to the manufacturer’s protocol. The GBrowse format to facilitate visualization of all uniquely cDNAs produced by reverse transcription were amplified mapped mRNA-Seq reads in the rice genome under Pi using a pair of gene-specific primers for each gene. stress (Fig. 6; http://rnaseq-pi.dna.affrc.go.jp/cgi-bin/gb2/ gbrowse/RNAseq_p/). A comprehensive rice transcriptome Sequencing and mapping of short reads onto the rice under Pi stress could be an indispensable tool in co- genome expression analysis of related genes, deciphering gene networks involved in Pi stress, and identifying genes that Extraction of total RNA, construction of cDNA, and maybe used to enhance P efficiency and improving P sequencing with the Illumina Genome Analyzer IIx were acquisition, which are all critical to crop productivity. performed as described previously (Mizuno et al. 2010). Trimming of Illumina adaptor sequences and low-quality bases (Q<20) at the 5′ and 3′ ends of each read was Materials and methods performed by the fastx_clipper program in the FASTX- Toolkit and an in-house C language program. Using the Plant materials and stress treatment pre-processed Illumina reads, the transcript structures were reconstructed by a series of programs, namely, the Bowtie Rice (Oryza sativa ssp. japonica cv. Nipponbare) seeds were version 0.12.4 for short-read mapping (Langmead et al. germinated and grown by hydroponic culture in nutrient 2009), TopHat ver version 1.0.13 for defining exon–intron media (Yoshida et al. 1976) in a growth chamber. After junctions (Trapnell et al. 2009), and Cufflinks version 0.8.2 14 days, the seedlings were subjected to Pi stress treatment by for gene structure predictions (Trapnell et al. 2010). transferring in a similar media with the following Pi Unannotated transcribed regions with no overlap to any concentrations: 0 mg Pi/L (18.7 mg NaCl/L; alternatively known loci was analyzed by comparing known transcript NaH PO )for –P stress, and 100 mg Pi/L (439 mg structures annotated in RAP-DB (http://rapdb.dna.affrc.go. 2 4 KH PO /L; alternatively NaH PO ) for ++P stress. The jp/) and the MSU Rice Genome Annotation Project 2 4 2 4 plants were maintained under Pi stress conditions for 30 days database (http://rice.plantbiology.msu.edu/) with the at 28°C and 70–80% humidity in a 16-h light/8-h dark cycle Cufflinks-predicted structures. Those regions that over- with the light period from 6:00 AM to 10:00 pm. The plants lapped with rRNA (1 transcript), tRNA (12 transcripts), were moved to different positions in the growth chamber and repetitive regions with similarity of ≥90% (3,888 several times during the light period so that all plants in each transcripts) in RAP-DB were discarded, and the remaining treatment were exposed to the same amount of light. Root unannotated transcripts were used for further analysis. To and shoot samples were collected at approximately 9:00 AM, estimate the expression levels of each transcript, all pre- frozen in liquid nitrogen, and stored at −80°C until processed reads were mapped into the IRGSP Build 5 subsequent analyses. pseudomolecules (http://rgp.dna.affrc.go.jp/E/IRGSP/ Build5/build5.html) by BWA (version 0.5.7) with default Measurement of dry weight and P uptake parameters (Li and Durbin. 2009). The expression level for each known and unannotated transcript was calculated as Shoots and roots of seedlings under +P (control), −P, and ++P RPKM (Reads Per Kilobase exon Model per Million conditions were separately sampled and oven-dried at 70°C mapped reads) based on the number of uniquely mapped for 2 days. The dry weight of each sample was determined and reads that completely overlap with the exonic regions. The the phosphorus content was analyzed by the phosphomolyb- resulting mRNA-Seq data were deposited in the DDBJ denum blue reaction using the U-0080D spectrophotometer Sequence Read Archive (accession no. DRA000314). (Hitachi, Tokyo, Japan) and Spectroquant phosphate test kit (Merck, Darmstadt, Germany). Rice 44 K microarray analysis Time course analysis of OsIPS1 expression Root RNA samples from plants under under+P (control), −P, and ++P conditions were used for gene expression profiling The expression of OsIPS1 was analyzed by RT-PCR. Total using the rice 44 K oligomicroarray (Agilent Technologies). RNA was extracted from root and shoot samples of 1–10- Labeling, hybridization and scanning were performed accord- ing to the manufacturer’s protocol as described previously day-old seedlings under –P condition using Trizol® reagent Rice (2011) 4:50–65 63 Fig. 6 Browser of the rice tran- scriptome under –P and ++P stress. The mRNA-Seq data mapped on the IRGSP Build 5 pseudomolecules can be accessed in a GBrowse format. The position of the transcript and corresponding microarray probe (if available) is indicated. For each transcript, the expres- sion profiles of mapped reads including RAP transcripts (e.g., OsIPS1) and unannotated tran- script under normal Pi concen- tration (control), –P and ++P stress, and recovery condition are shown as graphs. (Mizuno et al. 2010). Three biological replicates for each Analysis of responsive transcripts under Pi stress treatment were used for analysis. All the analysis data from microarray experiments are deposited in the NCBI Gene To detect stress response transcripts from all RAP transcripts Expression Omnibus (accession no. GSE25647). and unannotated transcripts, statistical analysis by G test was conducted between two stages of stress treatments for each Validation of gene expression by qRT-PCR transcript. The numbers of mapped reads on a given transcript and those on other regions for two stages were used as Expression of unannotated transcripts was analyzed by variables in 2×2 contingency tables for each test. All p values quantitative RT-PCR using the LightCycler® 480 system were corrected for a FDR (Benjamini and Hochberg. 1995) (Roche, Basel, Switzerland). Total RNA from root and shoot of 1% using the R package (ver 2.7.1) and in-house Perl was used for reverse transcription using the Transcriptor First- scripts (Benjamini and Hochberg. 1995). The GO terms Strand synthesis kit. The reaction mixture and reaction assigned to each transcript were obtained from RAP-DB for condition was described previously (Mizuno et al. 2010). each GO category for biological process, molecular function, The detection threshold cycle for each reaction was and cellular component. We also searched for the P1BS cis- normalized using Ubiquitin1 with 5′-CCAGGACAAGAT acting element in the 1-kb upstream regions of Pi stress- GATCTGCC-3′ and 5′-AAGAAGCTGAAGCATCCAGC-3′ responsive RAP-representative transcripts using the PLACE as primers. Three technical replicates for each treatment were database (http://www.dna.affrc.go.jp/PLACE/) and signal used for analysis. scan search program (Prestridge 1991). 64 Rice (2011) 4:50–65 Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory- Construction of transcriptome viewer efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25. We constructed a browser to show the distribution of Lejay L, Wirth J, Pervent M, Cross JMF, Tillard P, Gojon A. mapped mRNA-Seq reads for each condition along with the Oxidative pentose phosphate pathway-dependent sugar sensing as a mechanism for regulation of root ion transporters by reference genome sequence and RAP annotation in a photosynthesis. Plant Physiol. 2008;146:2036–53. GBrowse format Li H, Durbin R. Fast and accurate short read alignment with Burrows– Wheeler transform. Bioinformatics. 2009;25:1754–60. Li LH, Qiu XH, Li XH, Wang SP, Lian XM. The expression profile of Acknowledgments We thank Ms. F. Aota, Ms. K. Ohtsu, and Ms. K. genes in rice roots under low phosphorus stress. Sci China Ser C. Yamada for technical assistance, and Dr. Y. Nagamura and Ms. R. 2009;52:1055–64. Motoyama for assistance with the microarray experiments. This work Lin SI, Chiou TJ, Santi C, Jobet E, Lacut E, El Kholti N, et al. was supported by the Ministry of Agriculture, Forestry and Fisheries of Complex regulation of two target genes encoding SPX-MFS Japan (Genomics for Agricultural Innovation, RTR-0001). proteins by rice miR827 in response to phosphate starvation. Plant Cell Physiol. 2010;51:2119–31. Open Access This article is distributed under the terms of the Creative Lu TT, Lu GJ, Fan DL, Zhu CR, Li W, Zhao QA, et al. Function Commons Attribution Noncommercial License which permits any annotation of the rice transcriptome at single-nucleotide resolu- noncommercial use, distribution, and reproduction in any medium, tion by RNA-seq. Genome Res. 2010;20:1238–49. provided the original author(s) and source are credited. MacIntosh GC, Hillwig MS, Meyer A, Flagel L. RNase T2 genes from rice and the evolution of secretory ribonucleases in plants. Mol Genet Genomics. 2010;283:381–96. Miura K, Rus A, Sharkhuu A, Yokoi S, Karthikeyan AS, Raghothama KG, et al. The Arabidopsis SUMO E3 ligase SIZ1 controls phosphate References deficiency responses. Proc Natl Acad Sci USA. 2005;102:7760–5. Mizuno H, Kawahara Y, Sakai H, Kanamori H, Wakimoto H, Bari R, Datt Pant B, Stitt M, Scheible WR. PHO2, microRNA399, and Yamagata H, et al. Massive parallel sequencing of mRNA in PHR1 define a phosphate-signaling pathway in plants. Plant identification of unannotated salinity stress-inducible transcripts Physiol. 2006;141:988–99. in rice (Oryza sativa L.). BMC Genomics. 2010;11:683. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a Morcuende R, Bari R, Gibon Y, Zheng W, Pant BD, Blasing O, et al. practical and powerful approach to multiple testing. R Statist Soc Genome-wide reprogramming of metabolism and regulatory B. 1995;57:289–300. networks of Arabidopsis in response to phosphorus. Plant Cell Britto DT, Kronzucker HJ. Cellular mechanisms of potassium Environ. 2007;30:85–112. transport in plants. Physiol Plantarum. 2008;133:637–50. Mortazavi A, Williams BA, Mccue K, Schaeffer L, Wold B. Mapping Bustos R, Castrillo G, Linhares F, Puga MI, Rubio V, Perez-Perez J, et and quantifying mammalian transcriptomes by RNA-Seq. Nat al. A Central regulatory system largely controls transcriptional Methods. 2008;5:621–8. activation and repression responses to phosphate starvation in Nacry P, Canivenc G, Muller B, Azmi A, Van Onckelen H, Rossignol Arabidopsis. Plos Genet. 2010;6(9):e1001102. M, et al. A role for auxin redistribution in the responses of the Chen YF, Li LQ, Xu Q, Kong YH, Wang H, Wu WH. The WRKY6 root system architecture to phosphate starvation in Arabidopsis. transcription factor modulates PHOSPHATE1 expression in response Plant Physiol. 2005;138:2061–74. to low Pi stress in Arabidopsis. Plant Cell. 2009;21:3554–66. Pant BD, Musialak-Lange M, Nuc P, May P, Buhtz A, Kehr J, et al. Delhaize E, Ryan PR, Hocking PJ, Richardson AE. Effects of altered Identification of nutrient-responsive Arabidopsis and rapeseed citrate synthase and isocitrate dehydrogenase expression on microRNAs by comprehensive real-time polymerase chain internal citrate concentrations and citrate efflux from tobacco reaction profiling and small RNA sequencing. Plant Physiol. (Nicotiana tabacum L.) roots. Plant Soil. 2003;248:137–44. 2009;150:1541–55. Devaiah BN, Karthikeyan AS, Raghothama KG. WRKY75 transcription Paszkowski U, Kroken S, Roux C, Briggs SP. Rice phosphate factor is a modulator of phosphate acquisition and root development transporters include an evolutionarily divergent gene specifically in Arabidopsis. Plant Physiol. 2007;143:1789–801. activated in arbuscular mycorrhizal symbiosis. Proc Natl Acad Franco-Zorrilla JM, Valli A, Todesco M, Mateos I, Puga MI, Rubio- Sci USA. 2002;99:13324–9. Somoza I, et al. Target mimicry provides a new mechanism for Prestridge DS. SIGNAL SCAN: a computer program that scans DNA regulation of microRNA activity. Nat Genet. 2007;39:1033–37. sequences for eukaryotic transcriptional elements. Comput Appl Fujii H, Chiou TJ, Lin SI, Aung K, Zhu JK. A miRNA involved in Biosci. 1991;7:203–6. phosphatestarvation response in Arabidopsis.CurrBiol. Rubio V, Linhares F, Solano R, Martin AC, Iglesias J, Leyva A, et al. 2005;15:2038–43. A conserved MYB transcription factor involved in phosphate Gojon A, Nacry P, Davidian JC. Root uptake regulation: a central starvation signaling both in vascular plants and in unicellular process for NPS homeostasis in plants. Curr Opin Plant Biol. algae. Gene Dev. 2001;15:2122–33. 2009;12:328–38. Smith AP, Jain A, Deal RB, Nagarajan VK, Poling MD, Raghothama Hammond JP, White PJ. Sucrose transport in the phloem: integrating root KG, et al. Histone H2A.Z regulates the expression of several responses to phosphorus starvation. J Exp Bot. 2008;59:93–109. classes of phosphate starvation response genes but not as a He GM, Zhu XP, Elling AA, Chen LB, Wang XF, Guo L, et al. Global transcriptional activator. Plant Physiol. 2010;152:217–25. epigenetic and transcriptional trends among two rice subspecies Tesfaye M, Temple SJ, Allan DL, Vance CP, Samac DA. Overexpression and their reciprocal hybrids. Plant Cell. 2010;22:17–33. of malate dehydrogenase in transgenic alfalfa enhances organic acid Hsieh LC, Lin SI, Shih ACC, Chen JW, Lin WY, Tseng CY, et al. synthesis and confers tolerance to aluminum. Plant Physiol. Uncovering small RNA-mediated responses to phosphate defi- 2001;127:1836–44. ciency in Arabidopsis by deep sequencing. Plant Physiol. Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice 2009;151:2120–32. junctions with RNA-Seq. Bioinformatics. 2009;25:1105–11. Rice (2011) 4:50–65 65 Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren Wissuwa M, Gamat G, Ismail AM. Is root growth under phosphorus MJ, et al. Transcript assembly and quantification by RNA-Seq deficiency affected by source or sink limitations? J Exp Bot. reveals unannotated transcripts and isoform switching during cell 2005;56:1943–50. differentiation. Nat Biotechnol. 2010;28:511–U174. Xie Q, Frugis G, Colgan D, Chua NH. Arabidopsis NAC1 transduces Vance CP, Uhde-Stone C, Allan DL. Phosphorus acquisition and use: auxin signal downstream of TIR1 to promote lateral root critical adaptations by plants for securing a nonrenewable development. Gene Dev. 2000;14:3024–36. resource. New Phytol. 2003;157:423–47. Yi K, Wu Z, Zhou J, Du L, Guo L, Wu Y, et al. OsPTF1, a novel Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcription factor involved in tolerance to phosphate starvation transcriptomics. Nat Rev Genet. 2009a;10:57–63. in rice. Plant Physiol. 2005;138:2087–96. Wang C, Ying S, Huang H, Li K, Wu P, Shou H. Involvement of Yoshida S, Forno AD, Cock HJ, Gomez AK. Laboratory manual for OsSPX1 in phosphate homeostasis in rice. Plant J. 2009b;57:895– physiological studies of rice. 3rd ed. Manila: International Rice 904. Research Institute; 1976.

Journal

Rice – Springer Journals

Published: Nov 17, 2011

There are no references for this article.