Phosphorus remobilization from rice flag leaves during grain filling: an RNA‐seq study

Summary The physiology and molecular regulation of phosphorus (P) remobilization from vegetative tissues to grains during grain filling is poorly understood, despite the pivotal role it plays in the global P cycle. To test the hypothesis that a subset of genes involved in the P starvation response are involved in remobilization of P from flag leaves to developing grains, we conducted an RNA‐seq analysis of rice flag leaves during the preremobilization phase (6 DAA) and when the leaves were acting as a P source (15 DAA). Several genes that respond to phosphate starvation, including three purple acid phosphatases (OsPAP3, OsPAP9b and OsPAP10a), were significantly up‐regulated at 15 DAA, consistent with a role in remobilization of P from flag leaves during grain filling. A number of genes that have not been implicated in the phosphate starvation response, OsPAP26, SPX‐MFS1 (a putative P transporter) and SPX‐MFS2, also showed expression profiles consistent with involvement in P remobilization from senescing flag leaves. Metabolic pathway analysis using the KEGG system suggested plastid membrane lipid synthesis is a critical process during the P remobilization phase. In particular, the up‐regulation of OsPLDz2 and OsSQD2 at 15 DAA suggested phospholipids were being degraded and replaced by other lipids to enable continued cellular function while liberating P for export to developing grains. Three genes associated with RNA degradation that have not previously been implicated in the P starvation response also showed expression profiles consistent with a role in P mobilization from senescing flag leaves.


Introduction
Senescence is the final stage of plant development (Kong et al., 2013). It is typically characterized by visible leaf yellowing and is controlled by both environmental cues and internal elements such as phytohormones and gene regulators (Breeze et al., 2011;Dong et al., 2012). A range of biotic and abiotic environmental cues have been implicated in inducing senescence, including pathogen infection, drought, nutrition limitation and oxidative stress (Buchanan-Wollaston et al., 2003). The yellowing of senescing leaves is related to the degradation of macromolecules like chlorophyll, proteins and RNA, and this degradation ultimately results in a decline in photosynthetic activity (Buchanan-Wollaston, 1997;Guo et al., 2004).
Another key component of senescence in annual crop plants is the remobilization of mineral nutrients from vegetative organs to seeds (Kong et al., 2015). Phosphorus (P) is remobilized to developing seeds with particularly high efficiency, resulting in P harvest indices well above the carbon harvest indices in cultivars of most major crop species (Ara ujo and Teixeira, 2003;Batten, 1992;Rose et al., 2007Rose et al., , 2010. Although this may be ecologically advantageous because high seed P content improves seedling competitiveness in soils that are naturally low in bioavailable P (White and Veneklaas, 2012), in the context of agriculture it results in the removal of vast amounts of P from fields in harvested products.
The loading of P into grains is of critical importance for the global P cycle because the accumulation of P in the grains of the major cereals, oilseed and pulse crops globally removes the equivalent of over 50% of the P applied as fertilizer each year (Lott et al., 2000). Provided seed germination and seedling vigour can be maintained (see Pariasca-Tanaka et al., 2015 and discussion therein), reducing the amount of P stored in the grains of major crops may be a viable option to reduce the amount of P lost from the P cycle within agricultural production systems (Raboy, 2009;Rose and Wissuwa, 2012;Rose et al., 2013). Although breeding crop cultivars with reduced P levels in grains is an attractive solution, practical advances are hampered by our limited understanding of the physiology, and genetic and molecular regulation of P loading into developing grains, including P remobilization from senescing vegetative tissues (Wang et al., 2016).
The degradation of lipids, nucleic acids and proteins is a key component of the senescence process, and it is possible some of the genes involved in recycling P from phospholipids, RNA and proteins under P starvation in vegetative growth may also play a role in the mobilization of P from senescing tissues, as discussed in a number of recent reviews (Smith et al., 2015;Stigter and Plaxton, 2015). However, to the best of our knowledge, the only genes that have been definitively linked to P loading in developing grains are a putative sulphate transporter which has an endosperm-specific impact on grain P in barley (Raboy et al., 2014;Ye et al., 2011), the ATP-binding cassette (ABC) transporters ZmMRP4 and AtMPR5/AtABCC5 which, when inactive, result in low seed P in maize and Arabidopsis, respectively (Nagy et al., 2009;Shi et al., 2007), and the purple acid phosphatase (PAP) AtPAP26 which is involved in remobilization of P from senescing leaves to developing seeds in Arabidopsis (Robinson et al., 2012). NAC transcription factors are also involved in nutrient remobilization during senescence (Stigter and Plaxton, 2015); for example, the NAC transcription factor NAM-B1 plays a role in P loading into wheat grains, but it appears to play a broad role in nutrient remobilization to grains not limited to P (Uauy et al., 2006).
Given the key role of gene families such as P transporters (PTs), PAPs and SPX domain proteins (named after proteins SYG1/ PHO81/XPR1) in P recycling under P starvation stress Wu et al., 2013), and the role of RNases in remobilization of P from senescing hakea leaves (Shane et al., 2014), we hypothesized some members of these gene families may be involved in the remobilization of P from senescing leaves during grain filling. To test this hypothesis, we first identified a total of 115 P starvation-related (PSR) genes based on recent literature (Lin et al., 2009;Liu et al., 2011;MacIntosh et al., 2010;Secco et al., 2013) (see Table S1). We then investigated the expression of these genes in flag leaves during grain filling in the model cereal rice (Oryza sativa L.) using an RNA sequencing (RNA-seq) approach.
Several studies have investigated gene transcript abundance during leaf senescence in crops such as cotton and maize (Lin et al., 2015;Zhang et al., 2014) and in a range tissues during grain filling in rice (e.g. Duan and Sun, 2005;Zhu et al., 2003), including flag leaf tissue, and transcript data are publically available. However, without concurrent data on tissue P dynamics, it is not possible to put these data in context. Here, we quantified P dynamics in flag leaves during grain filling and subsequently examined gene expression at key time points when P concentrations in the flag leaf were stable compared with when the flag leaf was acting as a P source.

Remobilization of phosphorus from flag leaves during grain filling
Flag leaf P concentrations were relatively stable during early grain filling, but declined significantly from 9 days after anthesis (DAA) ( Figure 1). The decline in leaf P concentration from around 3.5 mg/g at anthesis to around 2.5 mg/g at 12 DAA suggests the leaf was acting as a 'P source' by 12 DAA. Two time points were selected for the RNA-seq study: 6 DAA, when the P concentrations were stable (flag leaf was not acting as a P source), and 15 DAA, when the flag leaf was clearly acting as a P source as shown by a 33% reduction in P concentration relative to 0 DAA. The 15 DAA time point was selected because it was the first time point that had a significantly lower P concentration than the 6 DAA time point.
Global gene expression of RNA-seq data Three biological replicates for each time point (6 DAA and 15 DAA) were analysed using RNA-seq. These six libraries generated 156 million 101-bp paired-end reads after quality control using FASTQC software; 83 million and 73 million from 6 DAA and 15 DAA, respectively (Table 1). Reads were filtered for adapter sequences, poly-N stretches and low-quality reads which resulted in 78 million (94%) and 70 million (95.6%) high-quality reads from 6 DAA and 15 DAA, respectively (Table 1). Of these highquality reads, around 60 million (77.5%) from 6 DAA and 56 million (79.7%) from 15 DAA were mapped to the reference genome (Table 1). There were 26 788 expressed genes at 6 DAA and 26 747 at 15 DAA (Figure 2a) and 5160 differentially expressed genes (DEGs) between the two time points using a cut-off based on P < 0.05 and false discovery rate (FDR) < 0.05 values.

GO enrichment analysis of DEGs
To classify the function of DEGs between two time points, gene ontology (GO) enrichment analysis was performed looking at DEGs at 15 DAA relative to 6 DAA. After filtering for DEGs based on P and FDR values (i.e. <0.05) as well as log 2 fold change (i.e. >1.5), 1643 DEGs were retained. Using parametric analysis of gene set enrichment (PAGE), 1180 DEGs were mapped to the GO term database resulting in a total of 263 enriched GO terms.
Up-regulated DEGs within the term biological process were largely associated with P metabolism, protein modification and transport, which was congruent with up-regulated DEGs within the molecular function term that were found to be associated with phosphotransferase, kinase or transport activity. Interestingly, a high number of up-regulated DEGs could be categorized to nine GO terms related to nucleotide, nucleoside or ribonucleotide binding proteins, although the genes found in each of the respective GO terms are largely identical. No up-regulated DEGs were enriched within the cellular component GO term ( Figure 3). Down-regulated DEGS within biological processes were associated with transcription and translation, while, congruently, down-regulated DEGs under molecular functions were associated with RNA binding and ribosome constitution. Several genes related to photosynthesis and photosystems were downregulated, which was further reflected in a marked downregulation of cellular component genes associated with plastids and organelles.
Expression of 115 PSR genes in the flag leaf at 6 DAA and 15 DAA Of the 115 PSR genes, 100 were expressed at 6 DAA and 97 at 15 DAA in flag leaves ( Figure 2b). The most highly expressed gene of the 115 PSR genes was the P transporter (OsPT) OsPT21 (Os01g0279700) with expression levels of 1613 FPKM (fragments per kilobase of transcript per million mapped reads) at 6 DAA and 1168 FPKM at 15 DAA. The second most highly expressed of the 115 PSR genes was OsPT24 (Os09g0570400) with 694 FPKM at 6 DAA and 334 FPKM at 15 DAA (Table S1). Only 14 of the 115 PSR genes were not expressed (FPKM = 0) at either 6 DAA or 15 DAA ( Figure 2b). We identified 38 PSR genes that were differentially expressed between the two time points (Table 2). Twenty-six of the 38 genes were up-regulated at 15 DAA when the flag leaf acts as a source of P for developing grains (Table 2). These 26 genes included OsPHO1 and OsPHO2, which are primary factors in P starvation signalling, as well as OsSPX1 and OsSPX2 domain  proteins that play a key a role in P homeostasis during vegetative growth Wu et al., 2013). Further, three PAPs (OsPAP1d, OsPAP3c and OsPAP10a) which are involved in the release of inorganic P (Pi) from organic P  and OsSQD2 which is involved in sugar transport to diacylglycerol (DAG) were among the 26 genes. Three of the 26 known P transporters (Liu et al., 2011), OsPT5, OsPT19 and OsPT20, were also among the 26 PSR genes up-regulated at 15 DAA in the flag leaf ( Table 2).
Expression of OsPAP, OsPT and OsSPX domain gene families in the flag leaf at 6 DAA and 15 DAA Given that several genes in the OsPT, OsPAP and OsSPX gene families were differentially expressed between 6 DAA and 15  , in addition to the three OsPAP genes included in the 115 PSR genes that were up-regulated at 15 DAA, a further three OsPAP genes (OsPAP9b, OsPAP15 and OsPAP26) had significantly up-regulated expression (log 2 fold change >1.3) at 15 DAA, although the expression of OsPAP15 was low compared to the other OsPAP genes (Figure 4a). Twenty-one of the 26 P transporters in rice (OsPT1 to OsPT26) were expressed in flag leaves. Three of these-OsPT5, OsPT19 and OsPT20-were differentially expressed at 15 DAA (Figure 4b). Expression levels of OsPT19 and OsPT20 were relatively low, and the expression level of OsPT5 was extremely low (Figure 4b). Interestingly, while OsPT21 and OsPT24 were not significantly differentially expressed between the time points, their overall expression was the highest (>400 FKPM) at both 6 DAA and 15 DAA (Figure 4b).
Six SPX genes (OsSPX1-OsSPX6) (Wang et al., 2009) and four SPX-MFS (major facility superfamily) genes (OsSPX-MFS1-4) are known to be regulated by OsPHR2 under P deficiency in rice (Lin et al., 2010;Wang et al., 2012;Wu et al., 2013). Examination of all OsSPX and OsSPX-MFS expression levels and their regulator OsPHR2 revealed that in addition to OsSPX1 and OsSPX2 (see above), OsSPX-MFS1 and OsSPX-MFS2 were also significantly upregulated at 15 DAA with log 2 fold changes of 2.55 and 2.80, respectively. The absolute expression level of OsSPX-MFS2 was greater than 200 FPKM at 15 DAA ( Figure 4c). Interestingly, expression of the transcription factor OsPHR2 was not significantly up-regulated at 15 DAA.

Discussion
While previous studies have profiled gene expression in a range of tissues during rice grain filling (e.g. Duan and Sun, 2005;Zhu et al., 2003) and leaf senescence in other crops such as cotton 10 000  (Lin et al., 2015) and maize (Zhang et al., 2014), lack of concurrent data on P remobilization precludes conclusions being made about the potential regulation of P remobilization from vegetative tissues during grain filling. In the present study, we identified critical developmental stages during grain filling when flag leaf P concentration was stable (6 DAA) and when the flag leaf was a P source (15 DAA) as the basis for subsequent RNA-seq studies.
Given the dearth of knowledge on specific genes involved in P remobilization from senescing leaves during grain filling (Wang et al., 2016), we investigated whether genes associated with P starvation during vegetative growth were also associated with P remobilization from vegetative tissues during grain filling. Of the 115 PSR genes identified in the published literature, 14 were not expressed in flag leaves at either developmental stage, which is not surprising given that some of these genes have tissue-specific expression. OsRNS1 and OsRNS7, for example, are specifically expressed in root tissue (MacIntosh et al., 2010). As expected, the absolute expression of the 115 PSR genes was low relative to the most highly expressed genes as many of the highly expressed genes (e.g. OsLIR1, photosystem subunits I, II and chloroplast precursor) are key genes for photosynthesis and are highly expressed in flag leaves during grain filling (Reimmann and Dudler, 1993;Sugano et al., 2010).
GO enrichment analysis is a common and powerful approach in gene expression analysis (Subramanian et al., 2005). By categorizing treatment or development-dependent global transcript changes into standardized GO terminology, metabolic shifts are revealed, allowing for the interpretation of possible changes that may occur at the molecular and cellular level (Du et al., 2010). Down-regulated DEGs were highly enriched in GO terms related to photosynthesis (GO:0015979) and plastids (45% of the cellular component category), indicating a reduction of photo-assimilatory anabolism at 15 DAA. Furthermore, down-regulated DEGs were enriched in GO terms related to transcription and translation (GO:0010467, GO:0006396, GO:0003723, GO:0003735, GO: 0006412, GO:0005840 and GO:00529), while up-regulated genes were enriched in GO terms related to protein modification  . This suggests that at 15 DAA, regulation of metabolism was achieved through the modulation of existing proteins rather than synthesis of new proteins, which is consistent with a change from anabolism towards catabolism and the remobilization of nutrients at 15 DAA. Resources are recycled and removed rather than re-invested within the flag leaf, which is supported by the finding that up-regulated DEGs were enriched in GO terms related to transport (GO:0006810, GO:0022857, GO:0022804, GO:0015291, GO:0005215, GO:0015297). Collectively, this clearly indicates a shift towards senescence at 15 DAA. Markedly enriched up-regulated DEGs were found in GO terms related to P metabolism (GO:0006796, GO:0006793, GO:0016310, GO:0016772 and GO:0016773), which corresponds to the observation of enhanced P remobilization from flag leaves at 15 DAA. However, none of the identified 115 PSR genes were found among the enriched genes within P metabolism-related GO terms. We therefore employed a targeted approach by individually comparing the expression levels of the PSR genes between 6 DAA and 15 DAA.
A number of PSR genes and other genes from PT, SPX and PAP families had expression profiles consistent with involvement in the remobilization of P from senescing leaves during grain filling. One of these genes, OsPAP26, is the rice orthologue of AtPAP26 which plays a role in P remobilization from senescing leaves to developing Arabidopsis seeds (Robinson et al., 2012). Interestingly, OsPAP26 has not been implicated in plant responses to P starvation at the vegetative stage, and therefore, its role in P remobilization may be specific to leaf senescence during grain filling (Stigter and Plaxton, 2015). In contrast, the other three OsPAP genes with highly up-regulated expression at 15 DAA in flag leaves, OsPAP3c, OsPAP9b and OsPAP10a, show increased expression in rice under P starvation at the vegetative stage (Secco et al., 2013;Zhang et al., 2011) and may have a broader role in P remobilization from leaves.
Two genes in the SPX-MFS subfamily OsSPX-MFS1 and OsSPX-MFS2 were highly up-regulated at 15 DAA. Wang et al. (2012) observed OsSPX-MSF1 mutant plants accumulated Pi in leaves which led them to suggest OsSPX-MFS1 is a P transporter in leaves during vegetative growth. Up-regulation of these same genes during senescence is consistent with them functioning as P transporters during remobilization of Pi from leaves during senescence. OsSPX1 and OsSPX2 were two of the 26 PSR genes that were up-regulated at 15 DAA, and although up-regulated expression may have arisen in response to the reduced levels of available Pi in the leaves, the similarity in structure and expression levels to OsSPX-MFS1 and OsSPX-MFS2 suggests these genes may have been in part responsible for the reduced levels of available Pi in the leaves. While both OsSPX-MFS1 and OsSPX-MFS2 are regulated by the transcription factor OsPHR2 under P starvation , up-regulation of OsSPX-MFS1 and OsSPX-MFS2 at 15 DAA in the present study did not correspond to increased expression of OsPHR2. This suggests that, as with the regulation of the Arabidopsis AtPAP26 (Stigter and Plaxton, 2015), key regulatory genes in the P starvation response are not necessarily involved in regulation of P remobilization from vegetative tissues during grain filling.
Three P transporters (OsPT5, OsPT19 and OsPT20) were upregulated while the expression of another four (OsPT8, OsPT14, OsPT16 and OsPT22) was down-regulated at 15 DAA. RNAi evidence has suggested OsPT8 is associated with P translocation from vegetative organs to reproductive tissues (Jia et al., 2011) which is contrary to the data presented here. A similar role has been assigned to AtPHT2; 1, an Arabidopsis orthologue of OsPT14 that was also down-regulated at 15 DAA. The expression level of two P transporters, OsPT21 and OsPT24, were similar at the two stages although their expression was much higher than all other P transporters at both 6 DAA and 15 DAA. It has been suggested that OsPT21-OsPT26 transport P between the cytosol and plastids in leaves and roots (Guo et al., 2008) so high expression of these genes may be indicative of the need to maintain Pi required for photosynthesis (Sivak and Walker, 1986). Despite a recent review of P remobilization during leaf senescence concluding that Pi transporters facilitate the translocation of P liberated from organic P sources to developing grain (Stigter and Plaxton, 2015), we found no evidence for up-regulation of any of the 26 known Pi transports in rice in transport of Pi from senescing flag leaves.
Several genes that participate in dephosphorylation such as OsPAP1d, OsPAP3c, OsPAP10c, OsPLDz2 and glycerophosphoryl diester phosphodiesterase and genes related to lipid metabolism, such as digalactosyldiacylglycerol (DGDG) and monogalactosyldiacylglycerol (MGDG) synthases, were up-regulated at 15 DAA. Enhanced activity of lipid phosphatases may be related to the activity of plastid-located MGDG and DGDG synthases for maintenance of membrane lipids in these organelles, which has been observed under low temperature stress (Campos et al., 2003). In particular, OsPLDz2 has been implicated in dephosphorylation of plasma membranes and endoplasmic reticulum phospholipids, generating DAG and releasing phosphate upon P starvation (Cruz-Ram ırez et al., 2006) with the phosphate becoming available for translocation to the developing grains. Moreover, three genes whose products are directly related to de novo synthesis of DAG, glycerol 3-phosphate dehydrogenase, dihydroxyacetone kinase and glycerol 3-phosphate O-acyltransferase were down-regulated, indicating that this process is reduced during senescence and highlighting the importance of phospholipid degradation for the release of DAG and Pi. Sulphoquinovosyl diacylglycerols substitute for phosphatidylglycerol to maintain the proportion of anionic lipids under P starvation (Yu and Benning, 2003). The replacement of phospholipids with sulpholipids and galactolipids is a key strategy used by plants from the Proteaceae family to maintain photosynthesis in P-limited environments . These observations therefore support the notion that the genes involved in recycling P from phospholipids during P starvation in leaves during vegetative growth are also involved in recycling P from leaf phospholipids during leaf senescence which could be translocated to the filling grain. Whether the replacement of phospholipids with other lipids successfully enables flag leaves to maintain high levels of photosynthesis, as occurs when phospholipids are replaced in P-efficient Proteaceae , or whether there is a resultant decline in photosynthesis, is not known.
RNA-P is a relatively large component of the metabolically active P pool in plants , and it was hypothesized RNS2 is part of the P scavenging system under P starvation stress (Abel et al., 2002). Shane et al. (2014) demonstrated RNase activity is highly up-regulated during leaf senescence in Hakea prostrata, and transcriptome analyses indicated that the S-like RNase, AtRNS2, plays a role in the P starvation response during vegetative growth, but is also induced by senescence in Arabidopsis (Hillwig et al., 2011;Taylor et al., 1993 that two genes related to RNA degradation were significantly upregulated at 15 DAA (log 2 = 1.8 and 0.9; Table 4). The genes coded by Os10g0556700 and Os10g0556600 correspond to CCR4-NOT transcription complex subunit 1, a deadenylase participating in poly (A) shortening (Temme et al., 2014), while Os07g0249600 codes for an mRNA-decapping enzyme 1B (Lykke-Andersen, 2002). These two enzymes act on two of the RNA protecting mechanisms, capping and polyadenylation, indicating that during rice leaf senescence, RNA degradation occurs leading to the release of P for its redistribution to the seeds. So far, neither of these two enzymes have been identified as induced by P starvation or leaf senescence; thus, our results widen our knowledge of alternative RNases participating in the liberation of P during leaf senescence. Interestingly, at 15 DAA, the expression of Os07g0661000 and Os04g0680400, that code for an adenosine deaminase and an allantoinase, respectively, were also upregulated, opening the possibility to associate the activity of the corresponding enzymes with that of CCR4-NOT. Deadenylation by the latter would lead to an increase in cytoplasmic adenosine which would be eventually degraded by the consorted activity of adenosine deaminase and allantoinase, to finally generate CO 2 and NH 3, which could be reassimilated via the glutamine oxoglutarate aminotransferase (GOGAT) pathway (Moffatt and Ashihara, 2002), indicating that nucleotide degradation also serves to recirculate N during leaf senescence.

Conclusions
Several genes involved in the phosphate starvation response were found to have gene expression profiles consistent with a role in the remobilization of P from senescing flag leaves. While some of these genes have previously been implicated in the P starvation response, it was evident that some of these genes (e.g. OsSPX-MFS1 and OsSPX-MFS2) are under different regulatory control during senescence. The fact that a number of genes that have not been implicated in the P starvation response were identified provides potential novel targets for manipulation to reduce the quantity of P loaded into cereal grains, thus reducing the amount of P removed from farmers' fields at harvest. Further studies are evidently needed to confirm the involvement of genes identified in the present study and to explore whether further novel genes or signalling molecules that have not been implicated in the P starvation response play a role in P remobilization from senescing leaves to developing grains.

Soil and plant material
Soil from the 0-to 10-cm horizon was collected from Southern Cross University's Brookside field site on the Lismore campus (28°49 0 3.74″ S, 153°18 0 21.49″ E). A detailed description of the soil is given in Rose et al. (2015). Briefly, the soil was pH 5.81 (1:5 H 2 O); total carbon 2.21%; electrical conductivity 0.48 dS/m; Bray P 0.6 mg/kg; and effective cation exchange capacity 13.6 cmol + / kg. The soil was air-dried and passed through a 2 mm sieve before being added at 8 kg/pot to 10 L, nondraining, black plastic pots. Immediately prior to transplanting rice seedlings (see below), nitrogen (N, as urea) and P (as superphosphate) were mixed into the soil for each pot at rates equivalent to 100 kg N/ ha and 50 kg P/ha, respectively, on a pot soil surface area basis, before soils were flooded with tap water.
Seeds of the widely grown rice (mega) variety IR64 were sterilized using HClO 3 for 2 min before germination in Petri dishes in the dark at 30°C for 2 days. Germinated seeds were transferred to a mesh floating on a solution containing 0.1 mM calcium (Ca) and 36 lM iron (Fe). After 10 days, the solution was changed to half strength Yoshida solution (Yoshida et al., 1976) and plants were grown for another 2 weeks, with nutrient solution replaced every week. Seedlings were then transplanted into soil with three evenly sized seedlings planted in each pot.
Plants were grown under controlled glasshouse conditions at Southern Cross University (Lismore, NSW, Australia) with a mean day/night air temperature of 29°C/21°C and relative humidity (RH) of 75%. Additional N (50 kg/ha as urea) was top-dressed at tillering stage (60 days after transplanting) to ensure low soil N did not induce premature leaf senescence.
Individual panicles in each pot were tagged at anthesis, with anthesis defined as when 50% of florets on a panicle had flowered. This occurred at around 85 days after sowing (DAS). At anthesis, and every 3 DAA until maturity, panicles of the same age (excluding panicles of the main tiller) from three separate pots were harvested and separated into grain, husk, rachis, stem, flag leaf and other leaves. This occurred until 33 DAA, when plants reached physiological maturity, with no plant sampled more than once. After weighing, fresh tissue samples were divided into half: one subsample was immediately frozen in liquid nitrogen and kept at À80°C for later RNA extraction and RNAseq analysis, and the second subsample was dried for 7 days at 40°C in a drying room.

Tissue biomass and P measurements
Phosphorus concentration in flag leaf tissue was measured by digesting 0.2 g tissue with 5 mL of nitric acid (HNO 3 ) in a MARS Xpress microwave oven (CEM Corporation, NC, Charlotte). After digestion, each sample was diluted by addition of purified water (Milli-Q, Millipore, Billerica, MA, US) to a final volume of 25 mL and P concentration measured by inductively coupled plasma mass spectrometry (ICP-MS) (Perkin Elmer, Melbourne, Victoria, Australia).

RNA extraction and library preparation
RNA extraction was undertaken on flag leaf tissue samples from two time points selected from the P audit (above). Total RNA was extracted from rice flag leaf tissue using the RNeasy Mini Kit (Qiagen, Victoria, Australia, Chadstone) according to the manufacturer's instructions. After extraction, RNA was quantified using the Nanodrop (ND1000; Labtech, Paris, France) and RNA quality was then examined using a 2100 Bioanalyzer (Agilent Technologies, CA, Santa Clara). High-quality RNA samples for library construction were selected based on 260/280 nm ratio and RNA integrity number (RIN) above 2.0 and 8.0, respectively. Samples that did not meet these quality criteria were re-extracted. The TruSeq mRNA stranded kit (Illumina, Scoresby, Vicotria, Australia) was used to prepare Illumina RNA-seq libraries for each sample from 1.3 lg of total RNA. The concentration of each library was measured by qPCR using the KAPA library quantification kit (KAPA Biosystems, MA, Wilmington) to determine the required number of PCR cycles for library amplification. The final concentration of the amplified library was measured using Qubit â Fluorometric Quantitation (Life Technologies, CA, Carlsbad). Library sequencing was undertaken with the Illumina Hi- National University (ANU), ACT, Australia. The raw sequencing data have been uploaded to ENA (European Nucleotide Archive) database (accession number: PRJEB11899).

RNA-seq data analysis
Raw sequencing reads in FASTQ format were first filtered for quality using FASTQC (Andrews, 2010) followed by removal of adapter sequences, poly-N stretches and low-quality reads using the BBDuk module of the BBMap software package (http://sourceforge.net/projects/bbmap/). All subsequent analyses were based on high-quality sequencing reads. Bowtie v2.2.4 (Langmead and Salzberg, 2012) was used to index the genome. Retained high-quality paired-end reads were mapped against the rice genome IRGSP 1.0 using TopHat (Trapnell et al., 2009), which is a splice aware aligner of RNA-Seq reads. The Ensembl Plants (http://plants.ensembl.org/Oryza_sativa/Info/Annotation) O. sativa cv. Nipponbare (ssp. japonica) reference genome annotation was utilized. TopHat identified the exon-exon junctions and produced the read vs genome alignment in BAM (Binary Alignment Map) format. Cufflinks (Trapnell et al., 2012) then used the TopHat-generated alignment to assemble a set of reference-based transcripts. Finally, the CuffDiff module of Cufflinks was used to identify differentially expressed genes between samples and CummeRbund (http://compbio.mit.edu/ cummeRbund/) R package was used for subsequent analyses.

Involvement of PSR genes in phosphorus remobilization from flag leaves during grain filling
To investigate whether any PSR genes were involved in remobilization of P from leaves during grain filling, we identified a total of 115 PSR genes, including 26 OsPT genes and eight RNS genes, based on recent literature (Lin et al., 2009;Liu et al., 2011;MacIntosh et al., 2010;Secco et al., 2013) (Table S1). First, we assessed the absolute FPKM value of the known 115 PSR genes at the two selected time points. We also compared their expression to the expression of the 30 most highly expressed genes to assess their relative level of expression. Secondly, we assessed the relative change in their expression between the two time points when the flag leaf had a stable P concentration compared with when the leaf acted as P source tissue.

GO analysis of DEGs using RNA-seq data
To analyse the GO enrichment analysis, firstly, significantly differential expressed genes were retrieved by cutting off based on P and FDR values <0.05, respectively, and log 2 fold change >1.5. Total of 1643 DEGs were retrieved of 5160 DEGs. GO terms for each rice gene were obtained using PAGE analysis in AgriGO public web tool (http://bioinfo.cau.edu.cn/agriGO/analysis.php? method=PAGE) which can be accepted an arbitrarily large input list with fold change (Du et al., 2010).
Expression of OsPTs, OsPAPs and OsSPX domains in the flag leaf at 6 DAA and 15 DAA Several OsPTs, OsPAPs and OsSPXs were among the 115 PSR genes with differential expression between the two time points. We therefore investigated the expression of all known members of these gene families at 6 DAA and 15 DAA, including those that were not among the 115 PSR genes. This included the OsSPX-MFS genes that are classified as SPX domain-possessing proteins that are regulated by Pi starvation in rice and preferentially expressed in the leaves (Lin et al., 2010;Wang et al., 2012).

Biological pathway analysis
We examined whether metabolic pathways that include RNase or phospholipase genes were more active when the flag leaf was acting as a P source (15 DAA) using the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway tool. We assumed that genes differentially expressed between 6 DAA and 15 DAA, particularly with elevated expression 15 DAA, might be involved in P remobilization. We therefore retrieved from the KEGG pathway database all the activating enzymes/genes related with lipid and nucleotide metabolism, such as glycerophospholipid metabolism, glycerolipid metabolism, sphingolipid metabolism, glycosphingolipid metabolism and RNA transport/degradation. Those retrieved genes were then mapped to the RNA-seq data to verify their involvement based on log 2 fold change and significance.