Spatio-temporal transcript profiling of rice roots and shoots in response to phosphate starvation and recovery.

Using rice (Oryza sativa) as a model crop species, we performed an in-depth temporal transcriptome analysis, covering the early and late stages of Pi deprivation as well as Pi recovery in roots and shoots, using next-generation sequencing. Analyses of 126 paired-end RNA sequencing libraries, spanning nine time points, provided a comprehensive overview of the dynamic responses of rice to Pi stress. Differentially expressed genes were grouped into eight sets based on their responses to Pi starvation and recovery, enabling the complex signaling pathways involved in Pi homeostasis to be untangled. A reference annotation-based transcript assembly was also generated, identifying 438 unannotated loci that were differentially expressed under Pi starvation. Several genes also showed induction of unannotated splice isoforms under Pi starvation. Among these, PHOSPHATE2 (PHO2), a key regulator of Pi homeostasis, displayed a Pi starvation-induced isoform, which was associated with increased translation activity. In addition, microRNA (miRNA) expression profiles after long-term Pi starvation in roots and shoots were assessed, identifying 20 miRNA families that were not previously associated with Pi starvation, such as miR6250. In this article, we present a comprehensive spatio-temporal transcriptome analysis of plant responses to Pi stress, revealing a large number of potential key regulators of Pi homeostasis in plants.


INTRODUCTION
P, in addition to being one of the most important macronutrients for all living organisms, is also one of the most limiting nutrients for plant growth and development. Although phosphorus (P) is often present in nonlimiting concentrations in the soil, the availability of Pi, the main source of P used by plants, is often hindered due to its low solubility in soil and high sorption capacity (Poirier and Bucher, 2002). As a consequence, application of Pi fertilizers has been the main alternative to overcome low soil Pi content to maintain and increase crop productivity. Yet, this process is not economically or environmentally sustainable as the price of phosphate rock, a nonrenewable resource, increases and run-off of excess Pi fertilizers damages the environment. Thus, studies to better understand the mechanisms involved in regulating Pi homeostasis in plants have been initiated to aid the development of crops with increased phosphorus use efficiency (Veneklaas et al., 2012;Wu et al., 2013).
Pi-limiting conditions trigger an array of adaptive responses at the morphological, physiological, biochemical, and molecular levels, such as reduced plant growth, altered root system architecture, and secretion of organic acids, phosphatases, and nucleases to increase Pi acquisition (Poirier and Bucher, 2002;Rouached et al., 2010;Chiou and Lin, 2011;Plaxton and Tran, 2011). With the emergence of microarray technology, several global transcriptome analyses have been performed in Arabidopsis thaliana (Misson et al., 2005;Morcuende et al., 2007;Müller et al., 2007;Bustos et al., 2010;Nilsson et al., 2010;Thibaud et al., 2010;Woo et al., 2012), and other plants, including rice (Oryza sativa; Wasaki et al., 2003Wasaki et al., , 2006Zheng et al., 2009), greatly improving our knowledge of the complex mechanisms regulating Pi homeostasis (Chiou and Lin, 2011;Péret et al., 2011;Plaxton and Tran, 2011;Jain et al., 2012;Secco et al., 2012). Briefly, it appears that the responses of plants to Pi starvation are highly regulated at the transcript level and are largely conserved between Arabidopsis and rice, with the Myb transcription factor At-PHR1 (PHOSPHATE STARVATION RESPONSE1), and its paralogs At-PHL1 (PHR1-LIKE) and Os-PHR2 in Arabidopsis and rice, respectively, being key regulators in the Pi starvation signaling pathway (Rubio et al., 2001;Zhou et al., 2008;Bustos et al., 2010). Under Pi-limited conditions, At-PHR1 is sumoylated by At-SIZ1 (small ubiquitin-like modifier E3 ligase) (Miura et al., 2005), inducing the expression level of numerous genes, including the microRNA miR399, which leads to the degradation of its target PHO2 mRNA that encodes an E2 ubiquitin conjugase enzyme, which in turn negatively regulates genes involved in Pi responses. In Arabidopsis, PHO2 degrades PHO1 (PHOSPHATE1; Liu et al., 2012), a protein involved in Pi (A) Schematic representation of the experimental design. Seeds were germinated in water for 2 d and then transferred to a Pi-sufficient hydroponic (0.320 mM) solution for 2 weeks before being either transferred to Pi-deficient media (0 mM) or maintained in Pi-sufficient solution for 21 d. After 21 d of treatment, half of the Pi-starved plants were supplied with Pi-sufficient media for up to a day. (B) Morphological appearance of rice seedling grown for 21 d on Pi-sufficient or -deficient media. Bar = 5 cm.
(C) to (G) Fresh weight (FW) ([C] and [D]) and Pi concentration ([F] and [G]) were assessed for shoots and roots. Root-to-shoot biomass ratio is represented (E). (H) Relative expression of Pi starvation-responsive genes MGD2, PT6, and SPX2, in the roots. Expression levels are relative to the +Pi condition. Errors bars are SD and n = 10, except for (H), where n = 3. Asterisk indicates statistically significant (P < 0.05) differences between samples grown under +Pi and -Pi conditions. transfer from the roots to the shoots, thus modulating Pi homeostasis . In addition, PHO2 indirectly affects the expression levels of some high-affinity Pi transporters, such as PHT1;8 and PHT1;9 Bari et al., 2006;Hu et al., 2011). The PHO2-miRNA399 pathway is also regulated by target mimicry, where the non-protein-coding transcript INDUCED BY PHOSPHATE STARVATION1 (IPS1) is induced by Pi starvation and has a complementary region to miR399, thus ultimately inhibiting the action of miR399-charged silencing complexes on PHO2 mRNA (Franco-Zorrilla et al., 2007).
Although these studies contributed to our understanding of the global responses to Pi starvation, they were limited to the genes present on microarrays and mainly focused on a few time points and specific tissues, thus failing to capture the dynamic and genome-wide range of responses triggered by such a stress. In addition to mRNA regulation, microRNAs (miRNAs) function as key regulators of nutrient stress signaling (Chiou, 2007), with miR399 and miR827 regulating Pi homeostasis (Fujii et al., 2005;Chiou et al., 2006;Hsieh et al., 2009;Pant et al., 2009;Lin et al., 2010;Kant et al., 2011;Kuo and Chiou, 2011;Wang et al., 2012). Yet, to date, only a single study was performed to assess miRNA regulation in rice under Pi stress, using plants starved for a short period of time (5 d) and focusing only on miR399 and miR827 (Jeong et al., 2011), thus revealing the need for the identification of Piresponsive miRNAs in rice.
The recent development of next-generation sequencing technology has enabled researchers to study whole transcriptomes at single-base resolution, in a multitude of samples/conditions, in a cost-effective manner. It is now possible to look at several levels of regulation, under a large number of conditions within a same study, thus providing an integrated understanding of the mechanisms involved in the regulation of Pi homeostasis. In addition, unlike microarrays, RNA sequencing (RNAseq) enables accurate quantification of mRNA and noncoding RNA expression, as well as revealing any previously unannotated transcripts and splicing isoforms (Trapnell et al., 2012). Using RNAseq to identify loci/isoforms associated with specific conditions, such as Pi deficiency, may lead to the discovery of key regulators of Pi homeostasis in rice (Trapnell et al., 2010).
In this study, we aimed to shed further light on the complex responses of rice plants to Pi stress using next-generation RNAseq coupled with a comprehensive time-course experiment, analyzing the early and late responses to Pi deprivation as well as Pi resupply. We generated a comprehensive and integrated data set that provides insight into the molecular responses of the rice transcriptome, and small RNAome, following Pi deprivation and resupply. This was accomplished with unprecedented resolution, thereby revealing key regulators of Pi homeostasis in plants.

Experimental Design and Physiological Responses
To better understand the complex mechanisms regulating Pi homeostasis in rice, we performed a time-course experiment, where pregerminated seedlings were grown hydroponically for 2 weeks on Pi-sufficient medium (0.32 mM Pi), before transferring half of the plants to Pi-deficient solution (0 mM Pi) for 21 d ( Figure 1A). After 3 weeks of Pi deprivation treatment, half of these plants where then resupplied with Pi-sufficient media for up to 24 h. In total, nine time points were selected to cover shortand long-term responses to Pi deprivation as well as the effects of Pi resupply on Pi-starved plants. To confirm the effectiveness of the Pi deprivation, physiological and molecular analyses were performed. Growth analysis showed that the shoot biomass was significantly affected by 21 d of Pi deficiency, being 2.3-fold lower than that of shoots grown under Pi-sufficient conditions, but we observed no differences in root biomass ( Figures 1B to 1D). Consequently, plants grown under Pi deprivation conditions showed a higher root-to-shoot biomass ratio compared with plants grown under Pi-sufficient conditions ( Figures 1B and 1E), which is in accordance with previous studies (Reymond et al., 2006;Jiang et al., 2007). We observed significant differences in the root-to-shoot ratio between the two conditions as early as 7 d after initiating the Pi treatment, reaching more than 2.2-fold difference after 21 d of Pi starvation.
Pi concentration analysis revealed that both shoots and roots showed a gradual decrease in Pi content as the length of the treatment increased ( Figures 1F and 1G). In the roots, Pi deprivation led to a 50% reduction in Pi concentration as early as 24 h after the initiation of the treatment, while it took 3 d to have similar effects in the leaves. Interestingly, even after 7 d of Pi treatment, the Pi concentration in both root and shoot was still twofold higher than that observed after 21 d of Pi starvation, suggesting that after 7 d of Pi deprivation, plants may not be completely starved (Figures 1F and 1G). The Pi concentration observed after 21 d of Pi starvation in the roots and shoots was similar to that reported in previous long-term Pi starvation studies (Wang et al., 2009a;Secco et al., 2010;Jia et al., 2011;Rouached et al., 2011). In addition, the shoot Pi concentration observed after 21 d of Pi starvation was similar to that of rice shoots when grown in P-poor soils (4.5 mmol/g fresh weight) (H. Shou, unpublished data). Pi resupply resulted in rapid increase of the internal Pi content in both roots and shoots. Indeed, within 6 h of Pi resupply, the shoot Pi concentration was already almost twice as high as that observed in plants grown under Pi starvation conditions and was similar to that of plants grown on Pi-sufficient media after 24 h of Pi resupply ( Figure 1F).
Finally, we used quantitative RT-PCR (qRT-PCR) to assess the molecular responses of rice plants to Pi deprivation and resupply ( Figure 1H). We measured expression of four known Pi starvationinduced marker genes, namely, the noncoding RNA IPS1 (Hou et al., 2005), the high-affinity Pi transporter PT6 (Ai et al., 2009), SPX2 (Wang et al., 2009b), and MONOGALACTOSYLDIACYL-GLYCEROL SYNTHASE2 (MGD2) involved in galactolipid synthesis (Kobayashi et al., 2009) in roots and shoots ( Figure 1H; see Supplemental Figure 1 online). Expression of these Pi starvationinduced marker genes increased with the duration of Pi deprivation, reaching a maximum after 21 d of Pi starvation, before returning toward basal level after Pi resupply ( Figure 1H). In the roots, IPS1 showed the strongest induction after long-term Pi starvation, being upregulated by more than 1000-fold, while MGD2 and PT6 were induced by more than 90 and 60 times, respectively. In the shoots, IPS1 and SPX2 displayed the highest induction upon long term Pi starvation, showing more than 60-fold the expression levels observed for plants grown under Pi-sufficient conditions (see Supplemental Figure 1 online).

Profiling the Short-and Medium-Term Responses to Pi Deprivation Using RNAseq
To analyze the effects of Pi status on the transcriptome of rice roots and shoots, we selected nine time points ( Figure 1A) and used three biological replicates per condition for RNAseq, representing a total of 126 libraries ( Figure 1A). Sequencing libraries were generated by first removing rRNAs and subsequently analyzing the remaining RNAs on an Illumina Hisequation 1000 using a 2X101 paired-end run. Overall, more than 3.6 billion paired-end reads passing filter were generated, and these were used for mapping onto the International Rice Genome Sequencing Project-1.0 (IRGSP-1.0) Nipponbare rice genome sequence. On average, each sample had 30 million paired-end reads, of which an average of 68 and 77% uniquely mapped to the rice genome (IRGSP-1.0) for root and shoot samples, respectively (see Supplemental Table   1 online). Thus, each sample had an average of 20 million uniquely mapped paired-end reads. The whole data set can be visualized at http://www.plantenergy.uwa.edu.au/annoj/Secco_2013.html using the AnnoJ genome browser as well as in the Rice Annotation Project Database (RAP-DB) genome browser.
To search for unannotated transcripts and isoforms, we made a reference annotation-based transcript (RABT) assembly using the Cufflinks software (Trapnell et al., 2012). To do so, and to potentially identify loci and isoforms that are either tissue or condition specific, we used mapped reads corresponding to the 21-d samples (+Pi and -Pi, root and shoot). Comparison of the newly generated assembly to the original Michigan State University (MSU) Rice Genome Annotation Project assembly (v7) assembly led to the identification of putative misannotated loci and previously unannotated loci. Here, we consider misannotated loci as loci assembled (A) Summary of the reference annotation-based transcript assembly and the corresponding mis-and unannotated loci identified. (B) Number of significantly differentially expressed genes after varying length of Pi deprivation in roots (R) and shoots (S) using Cufflinks (P < 0.05, FDR < 0.05). Number in parenthesis represents the number of previously unannotated loci among the differentially expressed genes. (C) Venn diagram representing the overlap of the early molecular responses in the roots upon Pi deprivation (1 h to 3 d). Number in parenthesis represents the number of genes upregulated by Pi deprivation. (D) AgriGO representation of the overrepresented GO terms in the 30 genes upregulated by 1 h of Pi deprivation in the roots was generated using singular analysis enrichment (Fisher's test, P < 0.05, FDR < 0.05). Number in parenthesis represents the FDR value. by Cufflinks but spanning two or more loci annotated in MSU v7. This is often due to the incorrect annotation of the initial gene model or annotated genes being too close. By contrast, previously unannotated loci are an assembled group of reads not overlapping with any known gene model annotation from MSU v7, thus including potentially true novel loci, unannotated portions of annotated loci, such as untranslated regions (UTRs) or exons, or genes not annotated in MSU v7, such as rRNAs. In total, the RABT assembly identified 948 putative misannotated loci and 2400 previously unannotated loci ( Figure 2A). Transcript abundance was then estimated by FPKM (for fragment per kilobase of exon per million fragments mapped) using Cuffdiff 2, a software known to robustly perform differential analysis using RNAseq data (Trapnell et al., 2012).
Differential expression analysis showed that only a few genes were significantly differentially expressed in either root or shoot within 24 h to 3 d of Pi deprivation, respectively (P < 0.05, false discovery rate [FDR] < 0.05) ( Figure 2B; see Supplemental Data Sets 1 to 3 online). In addition, these genes showed little overlap between the different early time points, suggesting either that the early steps of Pi deprivation sensing and signaling occur transiently or that these responses may be the consequences of nonspecific stresses due to the manipulation of the plants, and the Pi starvation sensing mechanism initiates later. Indeed, among the 30 genes upregulated by 1 h of Pi deprivation in the roots, only six of these were also induced after 6 h, 24 h, and 3 d of Pi deprivation ( Figure 2C). Interestingly, several genes involved in photosynthesis and carbon fixation were specifically induced after 1 h of Pi deprivation in the roots and not in the control plants (see Supplemental Data Sets 1 and 3 online). This is intriguing, as with the hydroponic system used in this study, the roots were kept in the dark during the entire experiment, except for less than a minute during the transfer of plants between nutrient solution replacements, and special care was taken to avoid contamination from green photosynthesizing tissues. Gene Ontology (GO) analysis of this subset of genes revealed a significant overrepresentation of genes involved in photosynthesis and the generation of precursor metabolites and energy, associated with an overrepresentation of ribulose-bisphosphate carboxylase activity ( Figure 2D). This set of genes was only upregulated at this very early time point (1 h) and was downregulated at the later time points. Among the 26 annotated genes downregulated by 1 h of Pi deprivation in the roots, several genes encoded functions involved in metal uptake, including two nicotianamine synthases, two metal transporters, and three thionins (see Supplemental Data Set 3 online). In the shoots, among the six significantly upregulated genes by 1 h of Pi deprivation, two genes encoded metallothioneins (LOC_Os12g38270 and LOC_Os12g38290). Metallothioneins affect metal tolerance and homeostasis and scavenge reactive oxygen species (Cobbett and Goldsbrough, 2002), which could be a mechanism to overcome the increase in certain ion concentration, such as iron, upon Pi starvation. Under Pi starvation, plants overaccumulated some ions, including iron (Misson et al., 2005;Hirsch et al., 2006;Lei et al., 2011), as a result of enhanced availability in the media.
After 6 h of Pi deprivation, several genes involved in iron homeostasis were upregulated in both roots and shoots. In roots, two putative vacuolar iron transporters (annotated as integral membrane proteins, LOC_Os04g45520 and LOC_Os09g23300) were among the six annotated upregulated genes, while in the shoots, two ferritins (LOC_Os11g01530 and LOC_Os12g01530) and one iron transporter (LOC_Os09g23300) were among the five induced annotated genes. The induction of iron storage genes and the repression of metal transporters (1 h in the roots) has been reported in Arabidopsis, and this is thought to be a response to iron overload due to low Pi (Misson et al., 2005;Thibaud et al., 2010). In Arabidopsis, a Pi starvation-responsive ferritin, Fer1, contains a P1BS element (Bournier et al., 2013), which is the binding site for PHR1, a major regulator of Pi starvation response. Further analysis confirmed that both PHR1 and its paralog, PHL1, could bind to Fer1, hence regulating the expression of Fer1 upon Pi starvation, demonstrating a direct link between iron and Pi homeostasis in Arabidopsis (Bournier et al., 2013). Promoter analysis of the genes involved in iron homeostasis that were induced by short-term Pi deprivation revealed that only LOC_OS12g01530, an ortholog of Arabidopsis Fer1, displayed a P1BS motif within the first two kilobases of its promoter, suggesting that the direct link between iron and Pi homeostasis observed in Arabidopsis is also conserved in rice. In addition, two genes (LOC_Os08g06010, a putative glycerol 3-phosphate permease, and LOC_Os03g40670, a putative glycerophosphoryl diester phosphodiesterase) are suggested to be involved in Pi remobilization (Cheng et al., 2011) and were also upregulated in the roots after 6 h of Pi starvation. The expression of these two genes gradually increased with the duration of Pi deprivation, reaching a maximum after 21 d and decreasing after Pi resupply, suggesting that this regulation is part of a continued response to Pi starvation. Within 24 h of Pi deprivation in the roots, 23 genes were upregulated, including two high-affinity Pi transporters, PT2 and PT4, which were significantly upregulated by more than 2.4-and 3.2-fold, respectively.
Medium-term Pi deprivation (3 d -Pi) led to a significant decrease in the internal Pi content by more than two-and threefold in both shoots and roots, respectively, but the number of differentially regulated genes remained moderately low, with only 28 and 179 genes being upregulated in the shoots and roots, respectively. Among the 179 significantly upregulated genes in the roots, several known Pi-responsive genes, including the high-affinity Pi transporters PT3 and PT6 as well as the noncoding RNA IPS1, were among the most induced genes, being upregulated by more than 11-fold. Surprisingly, in the shoots, medium-term Pi deprivation led to the upregulation of only three annotated genes, including IPS1 ( Figure 2B; see Supplemental Data Set 3 online).

Profiling the Long-Term Responses to Pi Starvation Using RNAseq
Long-term Pi deprivation (7 d or more) resulted in severe changes in gene expression, with 3761 and 3029 genes (2654 and 1950, using a twofold cutoff) being significantly differentially expressed after 21 d of Pi starvation in roots and shoots, respectively ( Figure  2B; see Supplemental Data Set 3 online). In the roots, genes such as the high-affinity Pi transporters PT3 and PT10, IPS1, MGD2, and a lipase (LOC_Os11g43760) were among the most induced genes, being induced more than 150-fold (Table 1; more than sevenfold [log 2 ]). It is also worth noting that several genes, including a hydrolase (LOC_Os05g46460), a putative low-temperature and salt-responsive protein (LOC_Os06g44220), and PT9, were only expressed under Pi-deficient conditions (FPKM > 30 in -Pi, while in +Pi FPKM = 0) ( Table 1). Among the genes most induced by longterm Pi starvation, several transcripts are not present on the rice Affymetrix arrays, including LOC_Os06g44220, encoding a putative low-temperature and salt-responsive protein, and thus constitutes potential key candidate genes involved in Pi homeostasis that have been missed by previous genome-wide studies. GO analysis of the subset of genes induced by 7 d of Pi deprivation in the roots revealed a significant overrepresentation (Fisher's test, P < 0.05, FDR < 0.05) of genes involved in lipid metabolism, response to oxidative stress, and inorganic anion transport (see Supplemental Figure 2 online). To assess the response to longterm Pi starvation in rice, when the internal Pi levels are at a minimum, the PAGEMAN tool (Usadel et al., 2006) was used to identify overrepresented (Fisher's test, P < 0.05) functions among the Only annotated genes with FPKM > 3 (P < 0.05, FDR < 0.05) are shown. Asterisks indicate genes with known function in Pi homeostasis. N/A indicates that the corresponding genes were not expressed in 1 Pi conditions but were induced in 2 Pi.indicates that no gene description is available for that locus, due to reads overlapping two loci. differentially expressed genes upon 21 d of Pi starvation in both roots and shoots ( Figure 3A). Interestingly, most of the genes upregulated by Pi starvation in roots and shoots were overrepresented in the same functional categories. Only a limited number of functional categories were overrepresented. We also observed that the genes upregulated under Pi starvation were enriched in lipid metabolism, phenylpropanoid metabolism, cytochrome P450 and transport functions, and, more specifically, Pi transport, as previously reported (Misson et al., 2005;Morcuende et al., 2007). By contrast, genes downregulated by long-term Pi starvation were enriched in genes involved in photosynthesis and N metabolism in both roots and shoots. Also, the genes downregulated by Pi starvation in the shoots were enriched in the protein synthesis functional category, mainly through the overrepresentation of rRNAs. Despite using rRNA-depleted RNA to generate the sequencing libraries, we still detected reads originating from rRNA sequences. This phenomenon is not surprising as rRNAs are highly expressed genes; thus, removal of 95 to 99% of rRNA still results in a significant number of reads. Additionally, N metabolism, including ammonia and Asp metabolism as well as tetrapyrrole synthesis, were all also enriched in the downregulated genes.
To compare the power of RNAseq to traditional microarray analyses, we assessed the proportion of differentially expressed genes identified by RNAseq in response to 21 d of Pi starvation in the roots having probe sets on the Affymetrix microarray chips (see Supplemental Figure 3 online). This analysis revealed that among the 3761 differentially expressed genes after 21 d of Pi starvation in the roots, 23% of these could not be detected by microarrays as a result of annotated genes lacking probe sets, genes being misannotated in MSU v7, or genes being previously unannotated; thus, RNAseq enabled the identification of a large number of key regulators of Pi homeostasis.

RNAseq Analysis of Pi Resupply to Pi-Starved Plants
To study the effects of Pi resupply, plants that were starved for 21 d were supplied with Pi sufficient media (0.32 mM Pi) for up to 24 h, allowing the classification of differentially expressed genes into eight groups according to their expression patterns in response to the Pi treatments ( Figure 3B). This classification revealed genes sets that responded to Pi starvation and recovered after Pi resupply, genes that responded to Pi starvation but persisted in their response during Pi resupply (classes 2 and 5), genes differentially expressed under Pi starvation, but not responsive to Pi resupply (classes 3 and 6), and genes that were not affected by Pi starvation but responded to Pi resupply (classes 7 and 8) (see Supplemental Data Sets 4 and 5 online).
Classes 1 and 4 genes respond to both Pi starvation and Pi resupply and thus likely function in Pi homeostasis. The proportion of genes returning to basal expression levels after Pi refeeding increased with the duration of the Pi resupply in both roots and shoots. Indeed, while 24% (270 genes) of the upregulated genes in the roots had returned toward basal level after 1 h of Pi refeeding, 61% (949 genes) of them had returned toward basal level after 24 h of Pi resupply ( Figure 3B). GO analyses (Fisher's test, P < 0.05, FDR < 0.05) of the three time points of class 1 revealed an overrepresentation of similar functional categories ( Figure 4A, green box), with the later time points having genes specifically enriched in lipid metabolism and catalytic activity ( Figure 4A). Thus, as expected, the genes belonging to this class were highly conserved among the three time points ( Figure 4B). Analysis of the 270 genes returning toward basal level after 1 h of Pi resupply revealed an overrepresentation of genes involved in the response to extracellular stimulus and cell communication ( Figure 4A). Among this set of genes, several key genes involved in Pi homeostasis were present, including some PTs, SPX genes, and genes involved in lipid metabolism, such as MGD2 ( Table 2). The presence of these genes, which are known to have key roles in Pi homeostasis, within this set (Table 2) suggests that the other genes in this group may have important roles in Pi homeostasis as well. A similar pattern was observed in the shoot, though to a lesser extent, with 1 h of Pi resupply resulting in 189 genes returning toward basal level. Among the 181 upregulated genes returning toward basal level, GO analyses (Fisher's test, P < 0.05, FDR < 0.05) revealed an overrepresentation of genes involved in phosphatase activity and included key Pi homeostasis regulators, such as some PTs, SPX genes, and genes involved in lipid metabolism, such as MGD2.
In parallel, the number of continuously positively regulated genes in the roots decreased from 24% (279 genes) of the total number of upregulated genes in the roots after 1 h of Pi refeeding to only 2% (34 genes) after 24 h of Pi resupply ( Figure 3B). In the shoots, a limited number of genes showed a continuous response, with 30 genes induced by 24 h of Pi resupply, compared with Pi-starved plants. Classes 3 and 6 are more likely to represent genes that are not directly linked to Pi homeostasis as they are induced or repressed during Pi starvation but do not respond to Pi resupply. While the number of genes belonging to the persistent negative class remained similar independently of the duration of the Pi resupply, the numbers of genes belonging to the persistent positive response class decreased with the duration of Pi resupply ( Figure 3B).
Finally, resupplying Pi to Pi-starved plants resulted in the differential expression of another class of genes, genes that were not affected by Pi starvation, but only responded to Pi resupply (classes 7 and 8) (see Supplemental Data Set 5 online). Indeed, 1 h of Pi resupply resulted in the differential expression of more than 1200 genes (789 when using a twofold cutoff) in the roots, with more than 900 (653 when using a twofold cutoff) being upregulated ( Figure 3B). The differential expression of such a large number of genes after 1 h of treatment suggests that sensing external Pi occurs rapidly and triggers a large set of responses, unlike the one observed in response to Pi starvation, where only 30 genes were upregulated in the roots after 1 h of Pi starvation. After 24 h of Pi resupply, the number of upregulated genes in the roots was reduced by threefold, with only 334 genes being induced ( Figure 3B). Interestingly, the overlap between the genes upregulated after 1 h of Pi resupply and the later time points was very low, with only 6% (65 genes) of the genes upregulated by 1 h of Pi resupply also being upregulated after 24 h of Pi resupply ( Figure 4C). GO analyses of the genes upregulated in the three time points confirmed an evolution in the response to Pi refeeding ( Figure 4A). Indeed, GO analysis of the 979 upregulated genes by 1 h Pi resupply revealed an overrepresentation (Fisher's test, P < 0.05, FDR < 0.05) of genes involved in transcription factor and kinase activity, representing 20% of the upregulated genes ( Figure 4A, black box). Among the various transcription factors induced by 1 h of Pi resupply, WRKYs and MYBs were the most represented with 14 and 10 members, respectively, and were significantly overrepresented (x 2 P < 0.05) ( Figure 4D). However, longer Pi resupply resulted in the upregulation of genes involved in protein synthesis, with the translation functional category being the most significantly overrepresented ( Figure 5A, yellow boxes).
While the number of genes specifically responding to Pi resupply decreased in the roots with the length of the Pi resupply, an opposite trend was observed in the shoots, with 33 to 978 genes upregulated by 1 and 24 h of Pi resupply, respectively, suggesting that the shoot responds more slowly to Pi resupply than the roots ( Figure 3B).

Identification of Conserved Pi Starvation-Induced Genes between Rice and Arabidopsis
To compare the responses of rice to long-term Pi starvation to those observed in other plants, we analyzed data from previous Arabidopsis genome-wide studies of long-term Pi starvation in roots or whole seedlings, including those of Misson et al. (2005), Morcuende et al. (2007), and Woo et al. (2012). From these data sets, we identified a set of conserved Arabidopsis genes induced in roots ( Figure 5A). Overall, 130 genes were found to be upregulated by Pi starvation in the three studies, thus defining a set of core Arabidopsis Pi starvation-induced genes. Rice orthologs were then identified using riceDB (Narsai et al., 2013), which integrates the InParanoid (Ostlund et al., 2010) and Gramene (Liang et al., 2008) orthology prediction methods, resulting in the  identification of 217 rice orthologs ( Figure 5B). Interestingly, of the 217 rice orthologs, 76 genes were represented in our set of 1543 upregulated genes in the root after 21d + 24 h -Pi, thus defining a set of genes that are induced by Pi starvation in both rice and Arabidopsis. Among these 76 genes, 95% (72 genes) were genes that were responding to both Pi starvation and recovery in rice (class 1) and included genes such as the high-affinity transporters, genes from the PHO1 and SPX families, and some phosphatases (see Supplemental Data Set 6 online).

Expression Patterns of the PT and SPX Genes
Since several PT and SPX genes were among the genes that were highly upregulated by Pi starvation, and returning toward basal level of expression (class 1) as early as 1 h after Pi resupply, we decided to examine their regulation in response to Pi stress in more detail. Surprisingly, while these two gene families have been shown to have key roles in Pi homeostasis, very little is known regarding their regulation upon Pi stress and only a few members have been characterized in plants. Among the 13 high-affinity Pi transporters encoded by the rice genome, only four have been characterized, namely PT1, PT2, PT6, and PT8 (Ai et al., 2009;Jia et al., 2011;Sun et al., 2012). Similar observations can be made with the SPX gene family where only SPX1 and SPX3 have been characterized (Wang et al., 2009a(Wang et al., , 2009b. While most of the PTs are preferentially expressed in the roots, PT1, PT4, and PT8 were also expressed in the shoots (FPKM > 2), though to a lesser extent ( Figure 6A). Changes in the expression of PTs in the shoots was largely unaffected by the first week of Pi deprivation, with only PT4 showing significant induction by more than 5.5-fold. Yet, after 21 d of Pi starvation, in addition to PT4, two others PTs were induced in the shoots, PT8 and PT5, with PT8 showing the highest expression level of all the PTs, with almost 50 FPKM observed. Pi resupply resulted in a return toward basal level as early as 1 h. In the roots, under Pi-sufficient conditions, PT1, PT2, PT4, and PT8 were the most expressed PTs (FPKM > 10), while the other PTs showed very little gene expression (FPKM < 0.5) ( Figures  6A and 6B). After 21 d of Pi starvation, eight of the PTs displayed a high level of expression (FPKM > 25), with PT3 and PT10 showing the greatest induction in gene expression, being upregulated by more than 1000-and 500-fold and reaching more than 100 and 25 FPKM, respectively. Interestingly, PT3 and PT10 were already induced by more than 15-fold as early as 24 h of Pi deprivation. In opposition, 1 h of Pi resupply resulted in a return toward basal level of expression for all of the PTs, except for PT1, which still showed a slight 1.5-fold increase during that period. However, after 6 h of Pi resupply, PT1 was back to prestarvation levels.
Among the SPX genes, SPX1 and SPX4 showed the highest level of expression in the roots under Pi-sufficient conditions (FPKM > 15) ( Figures 6A and 6B). Very little changes in gene expression could be observed during the first 3 d of Pi deprivation in the roots, with only SPX5 induced by 11-fold. After 21 d of Pi starvation, SPX1 displayed the highest expression level with more than 170 FPKM, corresponding to an induction of more than 10fold ( Figure 6B). However the greatest upregulation was observed for SPX5, which showed a 260-fold induction after 21 d of Pi starvation, and was associated with a FPKM of 125. In the shoots, all the SPX genes, with the exception of SPX4, were upregulated by Pi starvation, with SPX5 being the most upregulated gene showing more than a 1500-fold increase in expression. Interestingly, SPX2 was the most abundant SPX gene, with more than 140 FPKM after 21 d of Pi starvation.

Identification of Mis-and Unannotated Loci Involved in Pi Homeostasis
The generation of a reference annotation-based transcript assembly resulted in the identification of 948 misannotated loci and 2400 previously unannotated loci, when compared with the MSU v7 gene annotation, out of which 206 and 438 were significantly differentially expressed in at least one condition, respectively (Figure 2A; see Supplemental Data Set 7 online). While the aim of this study is not to generate an extensive list of loci or isoforms that require annotation, we provide a few examples of some previously unannotated loci/isoforms induced by Pi stress, highlighting the importance of improving the quality of the rice annotation to better understand the genome wide responses of rice to Pi stress.
Among the 206 putative misannotated loci that were differentially expressed under Pi starvation were loci such as XLOC_057362 spanning two annotated loci, namely, LOC_Os09g32840, encoding a phosphatase and LOC_Os09g32850, encoding a putative nucleotide pyrophosphatase, which was induced by more than 80-fold after 21 d of Pi starvation in the roots ( Figure 7B). The locus is supported by gapped reads, which arise from spliced transcripts, spanning the two original loci, and thus representing bona fide evidence of the presence of a single gene ( Figure 7A). While some putative misannotated loci may be correctly annotated in other gene annotations, such as in the new release of RAP-DB (Kawahara et al., 2013;Sakai et al., 2013), some remain misannotated, such as XLOC_057362.
Among the 276 previously unannotated loci differentially regulated after 21 d of Pi deprivation, only 35 were significantly upregulated by more than twofold (with a FPKM > 3, and size >400 bp). Interestingly, 14 of these can now be found in the recently updated RAP-DB, including genes such as IPS2, genes encoding for various hypothetical proteins, and non-protein-coding proteins, five corresponded to rRNAs, three were fragments of known genes, and 10 appeared to be previously unannotated transcribed loci. This result also validates the quality of the RABT assembly generated. Among the genes recently added in the latest version of RAP-DB, Os02g0609000 (XLOC_024511) is of particular interest ( Figure 7C). This transcript was one of the most highly induced genes in the roots with more than 1400-fold induction, as well as showing a 30-fold increase in the shoot. This locus also overlapped with the genomic location of O. sativa (osa)-miR827a, one of the key miRNAs involved in Pi homeostasis, which acts through the regulation of two target genes, SPX-MFS1 and SPX-MFS2 (Lin et al., 2010;Wang et al., 2012). This locus was considered to be the precursor of osa-miR827, but no cDNA or ESTs were reported for this locus so far (Lacombe et al., 2008), probably due to the fact that this gene is only expressed under Pi-deficient conditions. Furthermore, the transcript contains an open reading frame capable of coding for a protein of 10 kD (92 amino acids).
Among the previously unannotated loci identified, only three were differentially upregulated in the roots, namely, XLOC_054531, XLOC_030435, and XLOC_031078. XLOC_054531 showed the highest induction upon Pi starvation, with more than an 11-fold increase in gene expression observed after 21 d of Pi starvation. XLOC_030435, composed of two exons, appeared to be only expressed in the roots and was upregulated threefold after long term Pi starvation ( Figure 7D). Analysis of the coding potential of these three mRNAs showed that they could encode small proteins, yet no homology was observed with known proteins in the National Center for Biotechnology Information database. To determine if mRNAs corresponding to some of the loci were likely to be translated into proteins, total RNA from roots grown under Pi-deficient conditions for 30 d was fractionated into polysomal and nonpolysomal fractions. Among the four genes tested, Os02g0609000, XLOC_054531, XLOC_030435, and XLOC_031078, only XLOC_030435 showed enrichment in the polysomal fraction, thus suggesting active translation (see Supplemental Figure 4 online). In opposition, transcripts from XLOC_031078 and Os02g0609000 were enriched in the nonpolysomal fraction (see Supplemental Figure 4 online), suggesting that these transcripts are likely to be noncoding. This finding is consistent with the role of Os02g0609000 being the precursor of miR827. In addition, sequence analysis of XLOC_031078 showed high sequence similarity to Os01g0554400, a non-protein-coding transcript.
In the shoots, XLOC_045298 and XLOC_018325 showed the highest increase in gene expression after 21 d of Pi starvation with more than 16-fold and eightfold increases seen, respectively. The expression level of most of these genes returned to basal level after (A) Graphical representation of gapped reads, which are reads mapping across an intron. (B) Example of a putative misannotated locus, XLOC_057362, spanning two annotated genes, LOC_Os09g32840 and LOC_Os09g32850, and induced by Pi starvation. (C) XLOC_024511, a locus not present in MSU v7, referred to as Os02g0609000 in RAP-DB, is highly expressed in both root and shoot under Pi starvation and is the precursor of miR827a. (D) XLOC_030435, a loci neither annotated in MSU nor in RAP-DB, only expressed in the roots, and induced by Pi starvation. Top panels are screen captures from the AnnoJ genome browser (http://www.plantenergy.uwa.edu.au/annoj/Secco_2013.html), representing gene annotations from MSU v7, the RAP-DB IRGSP-1.00, and the RABT assembly as well as expression tracks. Bottom panels are screen captures from Integrative Genome Viewer representing individual RNAseq reads, where thick lines represents reads and thin lines represent gapped reads.
Pi resupply for 24 h, suggesting that these genes are most likely involved in Pi homeostasis and are not responding as a consequence of secondary effects.
We also identified previously unannotated loci that specifically responded to Pi resupply, with 61 induced upon Pi resupply (with a FPKM > 3, more than twofold change cutoff, and size >400 bp). One hour of Pi resupply led to the upregulation of 29 of these loci, with XLOC_006597 showing the highest upregulation (50-fold increase) in roots (see Supplemental Figure 5A online). This locus was still the most upregulated previously unannotated locus after 6 h of Pi resupply, being induced by >20-fold. This locus is likely to be a misannotated gene, as it appears to span multiple genes, namely, LOC_Os01g18360, LOC_Os01g18370, and LOC_Os01g18375, with LOC_Os01g18360, encoding an auxin-responsive element, thus forming a much longer gene than initially predicted by MSU or RAP-DB. Additionally, XLOC_018878, which was the second most highly upregulated locus among the previously unannotated loci after 1 h of Pi resupply, was not in RAP-DB, and represents a small 500-nucleotide-long transcript with a single exon (see Supplemental Figure 5B online).
We also identified genes with alternative transcription start sites use upon Pi stress, using the differential promoter use parameter in Cufflinks. In this way, PHO2 (LOC_Os05g48390), one of the key regulators of Pi homeostasis ( Figure 8A), was identified. PHO2 has an unusual gene structure, with a large 59 UTR containing several binding sites for the Pi starvation-responsive miRNA, miR399 ( Figure 8A). While PHO2 is not significantly regulated at the transcript level in the roots or shoots (see Supplemental Data Sets 1 and 2 online), an alternative isoform could be detected, specifically under Pi starvation in both root and shoot. We term these isoforms PHO2.1 and PHO2.2, with PHO2.2 being the alternative isoform we detected. PHO2.2 overlapped with LOC_Os05g48400, a gene annotated as an expressed protein. Yet, this locus is likely an artifact, as it is not present in the new release from RAP-DB, and gapped reads spanned the two annotated loci, suggesting that LOC_Os05g48400 is part of PHO2.2 ( Figure 8A). To confirm the existence of this isoform, the sequences encoding the N termini of both isoforms were cloned and sequenced (see Supplemental Figure 6 online), revealing that the only difference between the isoforms resided within their 59 UTRs, not altering the coding sequence of PHO2. Indeed, the first exon of the 59 UTR of PHO2.1 was replaced by two upstream fragments in PHO2.2 ( Figure 8A). Yet, both isoforms contained miR399 binding sites, suggesting that both isoforms are subject to miR399-mediated transcript cleavage ( Figure 8A). In addition, sequence analysis of the N-terminal region of the two isoforms revealed the presence of the P1BS element in the second exon of the 59 UTR of the PHO2.2 isoform. The P1BS element is the DNA motif to which the key transcription factor, PHR2, binds and is often found in Pi-responsive genes , potentially explaining the fact that PHO2.2 is only expressed under Pi starvation. To further confirm this expression pattern, we performed qRT-PCR analyses on shoot samples, showing that PHO2.2 was specifically expressed upon Pi starvation ( Figure 8B). In addition, to test whether the differences observed in the 59 UTRs of both PHO2 transcript isoforms could result in altered rate of translation, their association with polysomes was determined using RNA fractionated onto a Suc gradient. PHO2.1 was enriched in the nonpolysomal fraction, but PHO2.2 was preferentially associated with the polysomes, suggesting that PHO2.2 is more actively translated into protein than isoform 1 ( Figure 8C). The Pi starvation-induced mRNA isoform PHO2.2 is likely to constitute a new level of regulation of Pi homeostasis in rice, but further work will be required to identify the role of this isoform.
Among the other unannotated isoforms showing specific induction under Pi starvation were genes such as an ATP binding cassette transporter (LOC_Os06g37364), a nucleotide pyrophosphatase (LOC_Os09g32840), and a cytochrome P450 (LOC_Os01g52560). Interestingly, unlike with PHO2, some isoforms were only expressed in the root, including the cytochrome P450, LOC_OS01g52560 and LOC_Os09g32840, which encodes a nucleotide pyrophosphatase.

Genome-Wide Analysis of the MiRNAs
To increase our knowledge of the types of miRNAs involved in the long-term response to Pi starvation in rice, we generated and sequenced small RNAs libraries from roots and shoots of rice plants grown for 21 d under Pi-sufficient and -deficient conditions. On average, root and shoot samples had 13 and 60 million reads passing filter, respectively. Reads containing the 39 adapter sequence were then trimmed and size selected (between 15 to 30 nucleotides long) before being analyzed with miRanalyzer (Hackenberg et al., 2009). We found that 58 to 90% of the reads aligned to the rice genome (IRGSP-1.0), but only 0.7 to 11.5% of the reads could be mapped to known miRNAs (see Supplemental Table  2 online). Distribution of the mapped small RNA reads revealed that the majority of the shoot small RNAs had a size of 21 and 24 nucleotides, while most of the root small RNAs reads had a length of 22 and 24 nucleotides (see Supplemental Figure 7 online). In addition, out of the 547 rice mature miRNAs tested by miRanalyzer, an average of 170 to 370 could be detected in the roots and shoots, respectively. MiRNA analysis revealed that osa-miR166 was the most abundant miRNA in roots and shoots, accounting for 80 and 50% of all miRNA reads, respectively. Differential expression analysis revealed that osa-miR827 and osa-miR399 were the most induced families in both root and shoot, with miR827 being the most induced miRNA after 21 d of Pi starvation, with more than a 140-fold increase seen in roots ( Figures 9A and 9B; see Supplemental Data Set 8 online). These two miRNA families have been shown to be induced by Pi starvation in both Arabidopsis and rice Bari et al., 2006;Lin et al., 2010;Jeong et al., 2011;Wang et al., 2012), with Arabidopsis miR399 negatively regulating genes such as PHO2 and the two Pi transporters Pht1;8 Differentially expressed miRNAs in response to 21 d of Pi starvation in the roots (A) and shoots (B) were identified using the miRanalyzer software. Only significantly differentially expressed miRNAs with a fold change greater than 2 are shown (q value < 0.05). Asterisk indicates miRNAs that were not previously associated with Pi starvation. and Pht1;9 (Bari et al., 2006) and osa-miR827 regulating the expression of two putative Pi transporters, SPX-MFS1 and SPX-MFS2 (Lin et al., 2010;Wang et al., 2012), thus affecting Pi homeostasis. In addition to osa-miR399 and osa-miR827 families, 17 and 28 other miRNA families were differentially regulated upon 21 d of Pi starvation in root and shoot, respectively, corresponding to 33 unique Pi starvation-responsive miRNAs families (Figure 9). Comparison of these 33 rice miRNAs families to known Pi starvation-responsive miRNAs families identified in Arabidopsis and rapeseed (Brassica napus) (Hsieh et al., 2009;Pant et al., 2009;Kuo and Chiou, 2011) resulted in the identification of 20 miRNA families that were not previously associated with Pi stress (Figure 9). In addition, 80% of these Pi starvation-responsive miRNAs families were only found in rice (see Supplemental Table 3 online). Among these, miR6250 showed a seven-and two-fold induction upon Pi stress in the root and shoot, respectively ( Figures 9A and 9B). Expression of osa-miR6250 can be induced by arsenite treatment (Liu, 2012), and arsenate, the oxidation product of arsenite, is a chemical analog of Pi and can be taken up by Pi transporters such as PT8 . It is thus temping to hypothesize that osa-miR6250 could be involved in the regulation of both Pi and arsenate homeostasis in rice. Yet, the only predicted target gene for osa-miR6250 is LOC_Os04g45570, encoding an expressed protein, which is not altered at the transcript level by Pi starvation.
Several miRNAs also showed dramatic decrease upon Pi starvation, such as osa-miR169, osa-miR3979, and osa-miR1433. Of note, the Arabidopsis ortholog of osa-miR169, miR169, as well as osa-miR3979 were recently shown to be involved in the nitrogen starvation response, being strongly downregulated by nitrogen starvation (Jeong et al., 2011;Zhao et al., 2011). Crosstalk between nitrate and Pi homeostasis is well documented in Arabidopsis (Kant et al., 2011), with miR827 and its target gene NLA being required to maintain nitrate-dependent Pi homeostasis. Thus, it is possible that other miRNAs could have a similar function in rice, being involved in both nitrate and Pi homeostasis.
In conclusion, this study presents a high-resolution genomewide transcriptome analysis of the rice responses to Pi starvation and resupply in roots and shoots. Using RNAseq technology, we quantified expression of mRNAs and noncoding RNAs and determined the transcript structure, splicing patterns, and other posttranscriptional modifications of genes. We captured the full transcriptome dynamics associated with Pi stress in roots and shoots, which could not have been performed using traditional microarray technology because of a large number of genes not having probe sets or not being properly annotated. As a consequence, this study revealed the existence of potential key novel regulators of Pi homeostasis in rice. The RABT assembly also identified loci with different transcription start sites upon Pi stress, including PHO2, highlighting an additional level of complexity in the regulation of one of the key regulators of Pi homeostasis, through altered translation activity.
Classification of gene responses allowed us to disentangle the complex signaling pathways involved in Pi stress and identify genes that are directly involved in Pi deficiency and Pi resupply responses. So far, most of the previous studies focused on the genes directly involved Pi starvation, while neglecting the genes involved in response to resupplying Pi to Pi-starved plants. Yet, the latter condition is more comparable to plants grown in their natural environment, where they often confront Pi-limiting conditions and rarely encounter Pi-rich soil patches. Pi resupply induced many genes within a very short period of time, including transcription factors and kinases, suggesting the involvement of a numerous transcriptional and posttranslational modifications, ultimately leading to an increase in protein and amino acid synthesis.
Small RNAs regulate plant nutrition, including Pi homeostasis. Yet, only the study by Jeong et al. (2011) has addressed this question in rice, using plants Pi starved for a short period of time (5 d), thus potentially not detecting essential miRNAs that may be important later. Profiling miRNAs in plants starved for Pi for 21 d in root and shoot enabled us to confirm the well-characterized response of miR399 and miR827 to Pi starvation, but also to identify 20 miRNA families that were not previously associated with Pi starvation, including miR6250, that could be either directly involved in Pi homeostasis or involved in signaling crosstalk with other ions, such as arsenate.
This global analysis of the rice transcriptomic responses to Pi stress identifies candidate genes that can be further characterized to increase our knowledge of Pi homeostasis, bringing us one step closer to generating plants with increased Pi use and/ or acquisition efficiency.

Plant Materials and Growth Conditions
Rice (Oryza sativa cv Nipponbare) was used for all physiological experiments. Hydroponic experiments were performed under controlled conditions (day/night temperature of 30/22°C and a 12-h photoperiod, 200 µmol photons m 22 s 21 ), allowing 0.5 liters of hydroponic solution per plant. The hydroponic solution consisted of a modified solution as described by Yoshida et al. (1976), containing 1.425 mM NH 4 NO 3 , 0.513 mM K 2 SO 4 , 0.998 mM CaCl 2 , 1.643 mM MgSO 4 , 0.075 µM (NH 4 ) 6 Mo 7 O 24 , 0.25 mM NaSiO 3 , 0.009 mM MnCl 2 , 0.019 µM H 3 BO 3 , 0.155 µM CuSO 4 , 0.152 µM ZnSO 4 , and 0.125 mM EDTA-Fe, with or without 0.323 mM NaH 2 PO 4 , resulting in the +Pi and 2Pi conditions. The pH of the solution was adjusted to 5.5, and the solution was renewed every 3 d. Rice seeds were first pregerminated in tap water for 2 d before being transferred into the hydroponic solution, containing 0.323 mM Pi (+Pi) for 2 weeks. Half of the seedlings were then transferred to a solution lacking Pi (0 mM Pi) for 21 d before being resupplemented with 0.323 mM Pi for 24 h, while the other half of the seedlings continuously remained in +Pi conditions (control). During the resupply experiment, some rice seedlings were left in Pi-deficient media to serve as control. Roots and shoots were harvested separately at the following time points: 1 h, 6 h, 24 h, 3 d, 7 d, and 21 d after transfer to -Pi conditions, as well as 1, 6, and 24 h after Pi resupply. Furthermore, all experiment procedures such as media replacement and sample collection were performed at a similar time of the day (2 h after light) to minimize possible circadian effects.

Quantification of Pi Tissue Concentration
Determination of Pi in tissues was measured by releasing the cellular contents into water by repeated freeze-thaw cycles, or by incubation for 1 h at 85°C, and quantifying Pi by the molybdate assay according to the procedure of Ames (1966).

Total RNA Isolation and Library Preparation
The total RNA from the roots and shoots tissues was extracted using TRIzol reagent (Invitrogen), according to the manufacturer's instructions.
The integrity and quality of the total RNA was checked using NanoDrop 1000 spectrophotometer and formaldehyde-agarose gel electrophoresis. RNA was only used when the A 260 /A 280 nm ratio was >1.8.
For RNAseq library synthesis (three biological replicates per condition and one plant per replicate), total RNA was first depleted of rRNA using the Ribo-Zero rRNA removal kit (Plant Leaf and Plant Seed/Root kits; Epicentre). To do so, 1 µg of total RNA from root samples was used as input for rRNA removal, while 2 µg of total RNA was used for shoot samples. Sequencing libraries were generated using the TruSeq RNA sample prep kit (Illumina). An average of eight libraries were multiplexed and loaded on each lane of the Illumina Hiseq flow cell v3. Sequencing was then performed on a Hisequation 1000, as a 2 3 101 paired-end run, according to the manufacturer's instruction (Illumina).
For small RNAseq library generation (thee biological replicates per condition and one plant per replicate), 1 µg of total RNA was used as input for the Small RNA sample preparation kit (Illumina). In brief, specific 59 and 39 Illumina RNA adapters were sequentially ligated to small RNA molecules. Adapter-ligated molecules were then reverse transcribed and PCR amplified before being gel purified (corresponding to 15 to 30 nucleotides long). Libraries were then multiplexed (six libraries per lane) and sequenced for 51 cycles using the Illumina HiSequation 1000, as per the manufacturer's instructions.

Polysome Gradient Preparation
Polysome isolation was performed as previously described (Jabnoune et al., 2013). Briefly, roots from plants starved for Pi for 30 d were harvested, frozen, and ground to powder in liquid nitrogen. To prepare cytoplasmic extracts, 150 mg of powder was combined with 1.2 mL of chilled polysome buffer (100 mM Tris-HCl, pH 8.4, 50 mM KCl, 25 mM MgCl 2 , 5 mM EGTA, 18 mM cycloheximide, 15.5 mM chloramphenicol, and 0.5% [v/v] Nonidet P-40). Debris was removed by centrifugation at 16,000g for 15 min at 4°C. Aliquots of the resulting supernatant were loaded on 20 to 50% (w/w) continuous Suc gradients and centrifuged at 32,000 rpm for 165 min in a Beckman SW41 rotor at 4°C. After centrifugation, fractions were collected from the bottom to the top of the gradient with continuous monitoring of the absorbance at 254 nm, resulting in a polysomal and nonpolysomal fraction. RNA was extracted using Trizol (Invitrogen) following the manufacturer's instructions. A fixed volume of RNA samples corresponding to the polysomal and nonpolysomal fractions was then reversed transcribed and used for subsequent qRT-PCR.

Mapping of RNAseq Reads and Transcript Assembly and Abundance Estimation Using the Tuxedo Suite
All analyses were performed using the Cufflinks package (Trapnell et al., 2012), version 2.0.2 (Bowtie2 v2.0.0-beta7 and TopHat v2.0.3), and using the rice reference genome and gene model annotation file (GFF file) from MSU Rice Genome Annotation Project database v7 (http://rice. plantbiology.msu.edu/). A Bowtie2 index was generated using the rice IRGSP-1.0 genome. To align the RNAseq reads to the genome, TopHat was run with the following options: -b2-sensitive -r 0 -mate-std-dev 100 -g 1 -G, allowing only one multihit and providing a gene model annotation file (MSU v7).
The resulting aligned reads were then used to create a RABT assembly using Cufflinks v 2.0.2. Cufflinks was run in the discovery mode, allowing the identification of previously unannotated genes, as described by Trapnell et al. (2012). For this purpose, only reads corresponding to samples originating from the 21-d time point (+Pi and -Pi) were used. Assemblies were then merged using Cuffmerge, using the MSU v7 gene model annotation file as the reference annotation, resulting in a new transcript assembly, used to quantify transcript abundance.
Finally, transcript abundance (FPKM) and identification of differentially expressed genes was performed using Cufflinks with the default parameters (P < 0.05 and 5% FDR cutoff) with the addition of -N -b, corresponding to upper quartile normalization, and running a bias detection and correction algorithm. Differential transcript abundance at all genes was calculated as the logarithm base-2 of the expression ratio (FPKM minus Pi /FPKM plus Pi and FPKM Pi resupply /FPKM minus Pi ).
Mapped reads were then visualized using AnnoJ genome browser (http://www.annoj.org/), accessible at http://www.plantenergy.uwa.edu. au/annoj/Secco_2013.html. Because of the large amount of tracks loaded by default in AnnoJ, users can remove tracks of low interest to increase the loading speed of the Web page. Three different gene annotations are provided: the MSU v7, the RAP-DB, and the RABT assembly.
qRT-PCR cDNA was synthesized from 1 µg of total RNA using the Bio-Rad iScript reverse transcription kit, according to the manufacturer's instructions. Transcript levels were assayed using the LightCycler 480 (Roche). Absolute expression levels were normalized to the housekeeping gene LOC_Os06g11170 (Narsai et al., 2010), encoding a nucleic binding protein (see Supplemental Table 4 online).

Small RNA Abundance and Estimation
To identify and remove the 39 adapter sequence from raw reads, Cutadapt was used. For this, the first eight bases of the adapter sequence were used. Reads containing the adapter sequence were truncated up to the junction with the adaptor sequence. Reads were then filtered using the University of East Anglia Small RNA workbench to remove reads matching to known rice tRNAs and rRNAs, as well as remove reads that did not match to the rice genome (MSU v7) and discard reads with low complexity. Only reads between 19-and 24-nucleotides long were selected for further analyses. MiRNA abundance and differential expression analysis were performed using miRanalyzer (http://bioinfo2.ugr.es/miRanalyzer/miRanalyzer.php; Hackenberg et al., 2009), allowing zero mismatches to known miRNAs.

Accession Numbers
Illumina reads of all samples have been submitted to the Sequence Read Archive at the National Center for Biotechnology Information (http://www. ncbi.nlm.nih.gov/sra) under accession number SRA097415. Supplemental Data Sets 1 to 4 are deposited in the DRYAD repository: http://doi. org/10.5061/dryad.4480g.

Supplemental Data
The following materials are available in the online version of this article.