Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing.

Recent studies have demonstrated the important role of plant microRNAs (miRNAs) under nutrient deficiencies. In this study, deep sequencing of Arabidopsis (Arabidopsis thaliana) small RNAs was conducted to reveal miRNAs and other small RNAs that were differentially expressed in response to phosphate (Pi) deficiency. About 3.5 million sequence reads corresponding to 0.6 to 1.2 million unique sequence tags from each Pi-sufficient or Pi-deficient root or shoot sample were mapped to the Arabidopsis genome. We showed that upon Pi deprivation, the expression of miR156, miR399, miR778, miR827, and miR2111 was induced, whereas the expression of miR169, miR395, and miR398 was repressed. We found cross talk coordinated by these miRNAs under different nutrient deficiencies. In addition to miRNAs, we identified one Pi starvation-induced DICER-LIKE1-dependent small RNA derived from the long terminal repeat of a retrotransposon and a group of 19-nucleotide small RNAs corresponding to the 5' end of tRNA and expressed at a high level in Pi-starved roots. Importantly, we observed an increased abundance of TAS4-derived trans-acting small interfering RNAs (ta-siRNAs) in Pi-deficient shoots and uncovered an autoregulatory mechanism of PAP1/MYB75 via miR828 and TAS4-siR81(-) that regulates the biosynthesis of anthocyanin. This finding sheds light on the regulatory network between miRNA/ta-siRNA and its target gene. Of note, a substantial amount of miR399* accumulated under Pi deficiency. Like miR399, miR399* can move across the graft junction, implying a potential biological role for miR399*. This study represents a comprehensive expression profiling of Pi-responsive small RNAs and advances our understanding of the regulation of Pi homeostasis mediated by small RNAs.

Plants acquire mineral nutrients that are channeled into diverse biochemical and metabolic pathways essential for the completion of their life cycle. Phosphorus (P), one of the essential macronutrients of plants, is often not readily accessible to plants because of the low availability of phosphate (Pi) in the soil, the major acquired form of P (Marschner, 1995). To overcome this limitation, plants have evolved a number of responses to conserve and mobilize internal Pi and to enhance the acquisition of external Pi (Raghothama, 1999;Poirier and Bucher, 2002;Lin et al., 2009). Recently, the molecular basis underlying these adaptive responses has been an intense subject of research.
It is now recognized that small regulatory RNAs contribute to the complexity of gene regulation. MicroRNAs (miRNAs) and small interfering RNAs (siRNAs) of 21 to 24 nucleotides long are two major classes of small regulatory RNAs in plants, functioning as negative regulators of gene expression (Brodersen and Voinnet, 2006;Jones-Rhoades et al., 2006;Mallory and Vaucheret, 2006). miRNAs are processed from single-stranded RNA precursors capable of forming imperfectly complementary hairpin structures by the RNase III enzyme DICER-LIKE1 (DCL1) or DCL4. In most cases, plant miRNAs down-regulate gene expression by direct cleavage of their target transcripts, although translational repression has been proposed to be common as well (Brodersen et al., 2008). In contrast, siRNAs are generated from perfectly complementary long double-stranded RNA precursors, requiring the activity of RNA-dependent RNA polymerase.
Endogenous siRNAs can be further classified into trans-acting siRNAs (ta-siRNAs), natural antisense transcript-derived siRNAs, and heterochromatin siRNAs according to their genomic origin and the interaction of distinct components (Vaucheret, 2006;Vazquez, 2006). ta-siRNAs are phased 21-nucleotide RNA molecules whose synthesis is triggered by the cleavage of miRNA on the TAS transcript (Peragine et al., 2004;Vazquez et al., 2004;Allen et al., 2005). They act in-trans to down-regulate the expression of unrelated loci from which they are produced.
The biological functions of miRNA cover a broad range of physiological processes. In addition to their involvement in normal plant growth and development, the expression of many miRNAs is regulated by various biotic and abiotic stresses, which may contribute to the development of adaptive responses to overcome unfavorable conditions (Sunkar and Zhu, 2004;Navarro et al., 2006;Chiou, 2007;Li et al., 2008a;Lu et al., 2008). Recent findings have also demonstrated the important role of several miRNAs under nutrient deficiencies. miR395, miR398, and miR399 are up-regulated under sulfur (S), copper (Cu), and Pi deficiency, respectively. miR395 is involved in S assimilation and allocation by adjusting the activity of ATP sulfurylase (APS) and restricting the spatial expression of a sulfate transporter (AtSULTR2;1) responsible for root-to-shoot transport of sulfate (Jones-Rhoades and Kawashima et al., 2009). Other than a role in alleviating oxidative stresses (Sunkar et al., 2006), miR398 modulates the expression of Cu/zinc (Zn) superoxide dismutase according to the availability of Cu (Yamasaki et al., 2007). Moreover, cell-specific expression of miR169 is crucial in the development of nitrogen (N)-fixing nodules in Medicago truncatula (Combier et al., 2006). miR399 regulates Pi homeostasis by controlling the expression of PHO2 encoding a ubiquitin-conjugating E2 enzyme, UBC24 (Fujii et al., 2005;Aung et al., 2006;Bari et al., 2006;Chiou et al., 2006). It acts as a positive regulator to activate the Pi uptake and translocation when Pi supply is insufficient, while PHO2 functions as a negative regulator to suppress these activities in order to prevent Pi excess when external Pi is ample. Notably, recent results from reciprocal grafting suggest that miR399 may serve as a long-distance signal from shoots to suppress PHO2 expression in roots as a way of communication between the Pi status of shoots and the Pi transport activity of roots Pant et al., 2008). The long-distance movement of miR399 could be an early response to Pi deficiency. The shoot-derived miR399 is responsible for the degradation of PHO2 in the roots, where miR399 is not readily expressed .
New technologies of high-throughput sequencing have become powerful tools to disclose the large inventory of small RNA species in plants. These deep sequencing strategies may discover new miRNAs or novel small RNA classes and provide quantitative profiling of small RNA expression (Lu et al., 2005;Rajagopalan et al., 2006;Nobuta et al., 2007;Heisel et al., 2008;Zhu et al., 2008). In this study, Solexa sequencing was employed to profile the changes of small RNA expression in response to Pi deficiency. Very recently, an independent study was reported in which P-responsive miRNAs were identified using quantitative reverse transcription (RT)-PCR of primary transcript of miRNAs followed by small RNA sequencing (Pant et al., 2009). In our study, we not only profiled the expression of small RNAs but also looked into their regulatory mechanisms and potential functions. Importantly, we found two unique classes of small RNAs up-regulated by Pi deficiency and uncovered an autoregulatory mechanism of PAP1/MYB75 via miR828 and TAS4-siR81(2) that regulates the biosynthesis of anthocyanin. Changes in the abundance of specific groups of Pi-responsive small RNAs may mediate the activity or expression of target genes that are important in the development or regulation of adaptive responses, enabling plants to survive in a low-Pi environment.

RESULTS AND DISCUSSION
Overview of Small RNA Profiles in Response to Pi Deficiency Four small RNA libraries from Pi-sufficient roots (R+Pi), Pi-deficient roots (R2Pi), Pi-sufficient shoots (S+Pi), and Pi-deficient shoots (S2Pi) were constructed by Solexa high-throughput sequencing technology. After adaptor trimming, the numbers of sequence reads and unique sequence signatures from the raw data were calculated and then mapped to the Arabidopsis (Arabidopsis thaliana) genome (Supplemental Table S1). In summary, the sequencing yielded 4.8 to 5.1 million raw reads from each library, and approximately 90% of the reads remained after trimming the adaptor sequences. Approximately 3.5 million sequence reads (approximately 77%-83% of trimmed sequence reads) corresponding to 0.6 to 1.2 million unique sequence signatures could be mapped perfectly onto the genome (The Arabidopsis Information Resource 8). Only those perfectly mapped sequences were analyzed further. The sequence reads that mapped to single and multiple genomic locations were defined, respectively, as single-hit reads and multiple-hit reads. One-third of mapped sequences were single-hit reads in the root and one-half were single-hit reads in the shoot.
The length of the mapped sequences ranged from 16 to 27 nucleotides. We first examined the correlation between the length of small RNAs and the proportion of total sequence reads or unique sequence signatures (Fig. 1). Similar distributions were observed in +Pi and 2Pi libraries. Consistent with the earlier reports, there were two major peaks at 21 and 24 nucleotides in the total sequence reads of shoot libraries (Fig. 1A). Unexpectedly, an additional peak at 19 nucleotides was observed in the root libraries. This peak has not been noticed in the previous deep sequencing databases, probably because the pure root sample has not been examined. Further analysis revealed that the majority of these 19-nucleotide sequences originated from tRNAs (see below). When the unique sequence signatures were examined, the patterns of all four libraries were nearly identical (Fig. 1B). The 24-nucleotide small RNAs were dominant in either their reads or unique sequences, indicating that they are rich in sequence diversity.

Overrepresentation of tRNA-Derived Small RNAs in the Pi-Starved Roots
Mapped sequences were classified according to the annotation of the genome (Supplemental Table S2). Most sequences were mapped onto the nonannotated intergenic regions (30.12%-42.30%) and to a lesser extent onto MIRNA genes and protein-coding genes. Strikingly, a significantly higher proportion of small RNAs was mapped to tRNA genes in the root libraries (24.29%-34.32%) than in the shoot libraries (3.33%-5.58%). The percentage in the 2Pi root library (34.32%) was even higher than that in the +Pi root library (24.29%). When the origin of these small RNAs was further analyzed (Supplemental Table S3), we did not find any correlation between the abundance of small RNAs from specific tRNA species and their codon usage. We were surprised to find that a 19-nucleotide sequence cleaved from the 5# end of Gly-tRNA TCC represented over 80% of tRNA-derived small RNAs in the roots and accounted for up to 18.44% and 27.70% of total sequence reads in the +Pi and 2Pi root libraries, respectively, compared with only 1.00% to 1.79% in the shoot libraries (highlighted in Supplemental Table S3). For validation, we employed RNA gel analyses using the probes from three different segments corresponding to the 5# and 3# ends and the central anticodon loop of tRNA (Fig. 2B). The small RNAs derived from Asp-tRNA GTC were also examined using the 5# end probe. Consistent with the sequencing Figure 1. Global profiling of small RNAs. Correlations between the small RNA length and the proportion of total sequence reads (A) or unique sequence signatures (B) are shown. nt, Nucleotides. Figure 2. Overrepresentation of tRNAderived small RNAs in the Pi-starved roots. A, Small RNA gel analyses of tRNA-derived small RNAs in +Pi or 2Pi roots or shoots. B, The tRNA regions corresponding to the 5# and 3# end probes are designated by the black line, and the anticodon probe is designated by the gray line. The black arrowheads in A and B indicate the 5# cleaved products and corresponding cleavage site in tRNA, respectively. The vertical gray lines in A indicate the products of tRNA halves from the cleavage at the anticodon loop (gray arrowhead in B). nt, Nucleotides. data, the 19-nucleotide RNAs from Gly-tRNA TCC and Asp-tRNA GTC corresponding to the 5# end of tRNA were highly accumulated in the roots but much less in the shoots ( Fig. 2A). The 2Pi roots showed the maximum accumulation. When probing with the 5# or 3# end probe, additional bands at 30 to 40 nucleotides likely representing tRNA halves cleaved at the anticodon loop were observed preferentially in the roots ( Fig. 2A). These results revealed a spatial and temporal expression pattern of small RNAs derived from the specific cleavage on tRNA molecules rather than random degradation and may represent a novel processing of small RNAs.
Studies in microorganisms or mammalian cells have suggested that the tRNA-derived small RNAs are accumulated in a developmentally regulated manner and become dominant in specific tissues or under specific stress conditions (Lee and Collins, 2005;Haiser et al., 2008;Jochl et al., 2008;Li et al., 2008b;Thompson et al., 2008;Pant et al., 2009;Thompson and Parker, 2009). These small RNAs could be associated with the quality control of protein synthesis or regulation of gene expression. Moreover, a recent report showed that tRNA half molecules present in the phloem sap of pumpkin (Cucurbita maxima) were able to inhibit translational activity in vitro (Zhang et al., 2009). The authors proposed that these phloemdelivered tRNA halves may be a long-distance signal to coordinate the metabolic status between source and sink tissues. The differential accumulation of tRNAderived small RNAs between roots and shoots observed here may represent the consequence of such long-distance movement.

Expression Profiling of miRNAs in Response to Pi Deprivation
Forty-five and 55 differentially expressed MIRNA genes with greater than 1.5-fold relative change in sequence counts were identified in roots and shoots, respectively ( Fig. 3; Supplemental Table S4). MIR399s, MIR778, MIR827, and MIR2111 were clearly upregulated, whereas MIR169s, MIR395s, and MIR398s were down-regulated, in 2Pi roots and/or shoots. The up-regulation of various species from the MIR399 family by 2Pi was consistent with our previous observations , suggesting the reliability of our sequencing data. Among all MIRNA genes, the MIR399 family showed the highest degree of induction in both tissues.
The expression of miRNAs with relative changes greater than a ratio of 3 (highlighted in boldface in Fig.  3) was validated by RNA gel analysis. To examine whether the response of these miRNAs was specific to 2Pi or not, the RNA isolated from different nutrient deficiencies, including N, Pi, potassium (K), S, Cu, and iron (Fe), was inspected at once (Fig. 4). miR156 was highly expressed but unaffected by nutrient deficiencies in the shoots; however, it was up-regulated by 2Pi, 2N, or 2K, with the highest induction under 2N. miR156 regulates a group of SQUAMOSA PROMOTER-BINDING PROTEIN-LIKE (SPL) transcription factors involved in juvenile-to-adult vegetative phase transition and other development processes (Wu and Poethig, 2006;Gandikota et al., 2007;Schwarz et al., 2008;Wang et al., 2008). The physiological role for the up-regulation of miR156 in the 2Pi, 2N, or 2K root remains to be explored.
On the other hand, miR399, miR778, and miR827 were specifically up-regulated by 2Pi (Fig. 4A). It is interesting that several distinct small RNA sequences were generated by cleavage at different positions of the miR778 precursor (Supplemental Fig. S1). Because the corresponding miRNA* sequence could be detected for each cleaved RNA duplex, we suggest to rename miR778 as miR778-5p and the other two forms derived from the miR778 precursor as miR778-3p.1 and miR778-3p.2, respectively. Multiple processing of the miR778 precursor has not been reported previously, probably due to the low recovered sequence reads. All miR778 species were up-regulated by 2Pi, with the highest sequence count of miR778-3p.1 (Supplemental Table S4). The target genes of miR778-5p encode SET domain-containing proteins, which are involved in the regulation of histone methylation (Ebbs and Bender, 2006). Although the target genes of miR778-3p.1 and miR778-3p.2 await identification, generation of multiple miRNAs from the same precursor has increased the complexity for the miR778 regulation networks.
The target genes of miR827 encode SPX domaincontaining proteins. In yeast cells, SPX domaincontaining proteins participate in Pi transport or sensing (Lenburg and O'Shea, 1996). One of the target genes, At1g02860, encodes a nucleus-localized ubiquitin E3 ligase containing the RING and SPX domains, which was recently characterized as NITROGEN LIMITA-TION ADAPTATION (NLA) involved in the N recycle during N limitation (Peng et al., 2007) or as BENZOIC ACID HYPERSENSTIVE1 (BAH1) involved in immune responses (Yaeno and Iba, 2008). Here, we validated the targeting of miR827 to At1g02860 mRNA by RNA ligase-mediated (RLM) 5#-RACE analysis (Supplemental Fig. S2A) and observed the down-regulation of At1g02860 under 2Pi (Supplemental Fig. S2B). At1g02860 may be regulated by many factors and have multiple functions in response to diverse external stimuli. The miR827-mediated regulation of At1g02860 under 2Pi conditions may imply an interaction between N and P, although miR827 only responds to 2Pi but not 2N (Fig. 4A).
The expression of several down-regulated miRNAs was also confirmed (Fig. 4B). miR169 was downregulated in 2Pi, 2N, or 2S roots but to a lesser extent in the shoots. Down-regulation of miR169 was much more severe under 2N. miR395 was downregulated in 2Pi or 2N roots and shoots. miR398 was down-regulated by 2Pi, 2N, 2K, or 2Fe in both roots and shoots. In contrast, miR395 and miR398 were upregulated by 2S and 2Cu, respectively, which is in agreement with previous findings (Jones-Rhoades and Yamasaki et al., 2007). These observations suggest that miRNAs can coordinate cross talk among different nutrient-deficient responses. Noticing the upregulation of APS4 and SULTR2;1 under 2Pi (Supplemental Fig. S2C), we speculate that the suppression of miR395 under 2Pi may somewhat contribute to this effect. Up-regulation of APS4 and SULTR2;1 may result in the increase of sulfate assimilation and translocation for sulfolipid biosynthesis, a compensation for the reduced phospholipids during 2Pi. miR169 and its target gene, HAP2/NFYA transcription factor, are critical in the development of nodules (Combier et al., 2006). Additionally, miR169a and NFYA5 have been shown to be involved in drought tolerance (Li et al., 2008a). NFYA5 functions upstream to activate several drought-responsive genes, including those involved in oxidative stress.
The target genes of miR398 encoding Cu/Zn superoxide dismutase are also involved in the response to oxidative stress (Sunkar et al., 2006). Several nutrient deficiencies, such as N, Pi, K, and S, can trigger oxida- Change of target gene expression to Pi deficiency. The results were based on the microarry analysis from Genevestigator (https://www.genevestigator.ethz.ch/gv/index.jsp) except for the target genes of miR399 , miR2111-5p and TAS4-siR81(2) (Figs. 5 and 7), miR395, and miR827 (Supplemental Fig. S2), whose expression was analyzed by quantitative RT-PCR. Changes larger than 2-fold are indicated.
f Potential role of target gene under Pi deficiency. g Because of its low abundance, the expression of miR402 was validated by quantitative RT-PCR of primary transcript. tive stresses (Shin et al., 2005;Schachtman and Shin, 2007). Changes of miR169 and miR398 expression could be involved in adjusting a plant's ability to survive the oxidative stress caused by nutrient deficiencies. Interestingly, the opposite expression of miR398 in response to Cu and Fe deficiencies may reflect an elegant regulation of the preferential expression of different superoxide dismutase, Cu/Zn-superoxide dismutase or Fe-superoxide dismutase, coupling with different metal cofactors. These Pi-responsive miRNAs and the potential functions of corresponding target genes relative to adaptive responses to 2Pi are summarized in Table I. Up-regulation of miR399, miR778, and miR827 and down-regulation of miR398 by 2Pi were also observed in a recent study in which a quantitative RT-PCR platform was established to profile the primary transcripts of miRNAs (Pant et al., 2009). The changes of these primary transcripts of miRNAs agreed with those of their mature miRNAs in response to 2Pi, suggesting that they are mainly regulated at the transcriptional level. However, changes in miR156 and miR395 were not discovered in their system, probably because of differences in growth conditions or sample sources (Table I). To examine whether the regulation of these miRNAs depends on the PHOSPHATE STAR-VATION RESPONSE1 (PHR1) MYB transcription factor, which activates a group of Pi starvation-induced genes, including MIR399s (Rubio et al., 2001;Bari et al., 2006), we analyzed the existence of P1BS cis-element, the biding site of PHR1, within the 2-kb sequences upstream of miRNA hairpins. In addition to MIR399s, the P1BS element could be found in the upstream sequences of several species of MIR156, MIR169, MIR398c, and MIR827 genes (Table I), suggesting the possible regulation of these genes by PHR1.

Characterization of the Pi-Responsive miRNA2111 Family
We identified a miRNA gene family consisting of two genes (miR2111a and miR2111b) from the nonannotated sequence clusters, which were not reported as miRNAs during the time of our analysis. They were Figure 5. Characterization of miR2111 in response to Pi deficiency. A and B, Small RNA gel analyses of miR2111-5p and miR2111-3p under different nutrient-deficient conditions (A) and in small RNA biogenesis mutants (B). Col, Ecotype Columbia; Ler, ecotype Landsberg erecta; wt, wild type. C, Hairpin secondary structures of miR2111a and miR2111b precursors. D, RLM 5#-RACE analysis of miR2111-5p cleavage on At3g27150 mRNA. E, Quantitative RT-PCR analysis of At3g27150 mRNA in response to 2Pi. One of two biological replicates is presented, and the error bars indicate the SD of two technical replicates. specifically up-regulated by 2Pi but not by any of the other nutrient deficiencies tested (Fig. 5A). They possessed a potential hairpin secondary structure and could produce a small RNA duplex with twonucleotide overhang at the 3# end (Fig. 5C), the primary criterion of miRNA definition (Meyers et al., 2008). Because a substantial amount of both complementary strands of the RNA duplex was detected, it was difficult to differentiate miRNA from miRNA*; thus, we named them miR2111a-5p, miR2111a-3p, miR2111b-5p, and miR2111b-3p. The sequences of miR2111a-5p and miR2111b-5p were identical, while those of miR2111a-3p and miR2111b-3p differed by four nucleotides.
Although this family was also identified by other groups very recently (Fahlgren et al., 2009;Pant et al., 2009), its biogenesis and target genes have not been studied. The biogenesis of these miRNAs was analyzed using various mutants defective in the production of different classes of small RNAs (Fig. 5B). Here, we showed that the accumulation of these miRNAs was significantly reduced in the Pi-starved hen1 and dcl1 mutants. Nevertheless, some accumulation of miR2111 was still observed in the dcl1 mutant, suggesting that miR2111 can also be processed by another DCL, probably DCL4, in the absence of DCL1.
The gene At3g27150, encoding a Kelch domaincontaining F-box protein, was predicted to be a target of miR2111-5p. The cleavage of the mRNA of At3g27150 was confirmed by RLM 5#-RACE analysis (Fig. 5D). Interestingly, we did not observe the negative correlation between the expression of miR2111-5p and At3g27150. Instead, the steady-state level of At3g27150 mRNA was moderately induced by 2Pi (Fig. 5E) regardless of the enhanced expression of miR2111-5p. Several reports have shown a positive temporal correlation between the expression of a miRNA and its target gene in spite of the negative regulatory role of miRNAs (Llave et al., 2002;Achard et al., 2004;Kawashima et al., 2009). The consequences of different miRNA regulatory circuits were classified into four categories, spatial restriction, mutual exclusion, dampening, and temporal regulation, in which the final expression level of the miRNA and its target gene may not be correlated (Voinnet, 2009). We hypothesize that miR2111-5p may regulate the spatial expression of At3g27150 by restricting its expression in specific cells and/or fine-tuning the expression of At3g27150 to maintain an optimal expression level. Alternatively, the possibility of translational repression cannot be ruled out. Nevertheless, regulation by miR2111 may be only one of many regulatory routes of At3g27150, because multiple regulatory pathways integrated from different levels are usually involved in gene regulation. It is important to note that the target genes of miR399, miR827, and miR2111-5p encode proteins involved in ubiquitin-mediated protein degradation, suggesting a crucial role of posttranslational regulation of protein abundance in the control of 2Pi responses.
The homologs of the miR2111 gene could be identified in many other dicotyledon plant species, such as Brassica rapa, M. truncatula, Lotus japonicus, and Vitis vinifera (Supplemental Fig. S3), but they could not be found in monocotyledon plants. This observation suggests that either the miR2111 gene had become lost in monocotyledons or, more likely, that it had emerged after the separation of the dicotyledons from the monocotyledons. The sequence of miR2111-5p is more conserved than that of miR2111-3p. The four miR2111 genes in M. truncatula are clustered in the same direction within a range of 40 kb. Interestingly, there are only 32 nucleotides between the middle two hairpins, which are likely derived from a polycistronic transcript. Conservation of the miR2111 species im-plies the importance of miR2111-mediated regulation for survival under the 2Pi environment.

Pi-Responsive Small RNAs Derived from LTR of Retrotransposon
A cluster of small RNAs that originated from the long terminal repeat (LTR) of Copia95 retrotransposon was identified as showing a differential expression pattern in response to 2Pi (Fig. 6A). This small RNA cluster is located in the intergenic region of At5g27990 and At5g28000 containing a multiple array of truncated LTRs either as inverted or tandem repeats, suggesting the occurrence of several duplication events (Fig. 6C). The arrangement of inverted repeats allows the formation of a hairpin structure. The length of this hairpin is relatively long, at approximately 250 nucleotides. Unlike those miRNA-like small RNAs recently found in rice (Oryza sativa) grains (Heisel et al., 2008;Zhu et al., 2008), which are cleaved from a long hairpin, the small RNAs generated from this LTR hairpin are not phased. smRPi1 LTR , the most abundant small RNA from the hairpin stem, accumulated spe- Figure 7. Autoregulatory mechanism of PAP1/MYB75 via miR828 and TAS4-siR81(2). A, Small RNA gel analysis of TAS4-siRNAs under different nutrient-deficient conditions. B, Small RNA gel analysis of TAS4-siR81(2) and miR828 in the shoots of wild-type (wt) plants and tas4, mir828, and PAP1/MYB75 activation (pap1-D) mutants. C, Expression of PAP1/MYB75, PAP2/ MYB90, and MYB113 by quantitative RT-PCR analysis in tas4 and mir828 mutants. One of two biological replicates is presented, and the error bars indicate the SD of two technical replicates. D, Alteration of anthocyanin accumulation in tas4, mir828, and pap1-D mutants. Error bars represent the SD (n = 5). FW, Fresh weight; OD, optical density. E, A proposed model of the autoregulatory pathway of PAP1/ MYB75 via miR828 and TAS4-siR81(2). cifically in 2Pi roots but not in other nutrient-deficient tissues (Fig. 6A). DCL1 is essential for the biogenesis of smRPi1 LTR , and its accumulation was not reduced in the dcl2,3,4 triple mutants or in the rdr1, rdr2, and rdr6 mutants (Fig. 6B). smRPi1 LTR is likely generated from the cleavage of a single-stranded rather than a double-stranded RNA precursor, because the passenger strand of its duplex was identified from the hairpin stem with two mismatched base pairs. Because we were unable to detect smRPi1 LTR in the Landsberg accession, we hypothesize that smRPi1 LTR is a newly evolved small RNA due to rapid rearrangement of LTR and may represent an intermediate small RNA species from the transition of siRNAs to miRNAs. AtCopeg1 (for Copia evolved gene1; At2g04460), a protein-coding gene originating from a truncated AtCopia95 retrotransposon, was recently shown to be highly up-regulated by 2Pi . It will be interesting to elucidate the relationship between AtCopeg1 and smRPi1 LTR in terms of coevolution or cross-regulation.

Up-Regulation of miR828-and TAS4-Derived siRNAs by Pi Deficiency and Autoregulation of PAP1/MYB75
The sequencing results revealed the up-regulation of ta-siRNAs derived from TAS4 and miR828, which initiates the cleavage of TAS4 RNA in 2Pi shoots (Supplemental Table S4). These ta-siRNAs were upregulated in the shoots by 2Pi as well as 2N (Fig. 7A). TAS4-siR81(2) targets a group of MYB transcription factors, PAP1/MYB75, PAP2/MYB90, and MYB113, that are involved in the biosynthesis of anthocyanin (Rajagopalan et al., 2006). It is well documented that plants accumulate anthocyanin under 2Pi, in which these MYB transcription factors are up-regulated (Misson et al., 2005). The negative regulatory role of siRNA seemed to be contradictory to the positive correlation between the expression of TAS4-siR81(2) and its target genes under the same conditions. This prompted us to carry out further investigation. TAS4 and MIR828 T-DNA knockout lines and a PAP1/ MYB75 activation mutant, production of anthocyanin pigment1-Dominant (pap1-D; Borevitz et al., 2000), were examined. As predicted, there was no detectable TAS4-siR81(2) in the tas4 mutant, and neither miR828 nor TAS4-siR81(2) was detected in the mir828 mutant (Fig. 7B). The accumulation of miR828 and TAS4-siR81(2) was already observed in the +Pi shoots of pap1-D mutant and increased under 2Pi (Fig.  7B). This suggests that PAP1/MYB75 can positively regulate MIR828 and/or TAS4 genes either by activating their expression directly or indirectly or by suppressing a negative regulator of MIR828 and/or TAS4 via a different feedback loop. Analyses of the 1-kb upstream sequences of MIR828 and TAS4 genes revealed multiple MYB binding sites (predicted by the PLACE database; http://www.dna.affrc.go.jp/PLACE/ index.html) and one PAP1 cis-regulatory element, (C/T) (A/C)NCCAC(A/G/T)N(G/T) (Dare et al., 2008), at 305 and 747 bp upstream of the MIR828 hairpin precursor and TAS4, respectively. Thus, PAP1 likely acts directly on the MIR828 and TAS4 promoters to activate their transcription. On the other hand, the transcript level of PAP1/MYB75, PAP2/MYB90, and MYB113 was elevated in the tas4 and mir828 mutants under both +Pi and 2Pi conditions (Fig. 7C). This could be explained by the relief of targeting of TAS4-siR81(2). Consistent with the increased expression of these MYB transcription factors, the anthocyanin accumulation was enhanced in the shoots of tas4 and mir828 mutants under both +Pi and 2Pi conditions (Fig. 7D).
These results uncovered an autoregulatory mechanism of PAP1/MYB75 via miR828 and TAS4-siR81(2) in which an adequate expression level of PAP1/ MYB75, PAP2/MYB90, and MYB113 and a proper accumulation of anthocyanin are maintained during 2Pi (Fig. 7E). In this model, Pi deficiency results in up-regulation of those MYB transcription factors, which subsequently activate the biosynthesis of anthocyanin. On the other hand, the increased level of PAP1/MYB75 can also lead to production of TAS4-siR81(2) via the activation of miR828 and/or TAS4, which antagonizes its own expression and other MYB transcription factors. Whether the up-regulation of miR828 can be bypassed by the activation of PAP1/ MYB75 remains to be studied. Because anthocyanin is a common stress pigment, we hypothesize that this autoregulation machinery may not be restricted to 2Pi and could also be applied to other stress conditions in which anthocyanin is accumulated.

miRNA399* and Its Movement
In addition to miR399s, we observed a considerable amount of miR399* accumulated in 2Pi tissues. The read number of different miR399* species was about half that of miR399 (e.g. miR399f*) or even higher than that of miR399 (e.g. miR399a* and miR399d*). Such high accumulation of miR399* could not be neglected because it was more abundant than many other miRNAs. RNA gel analyses confirmed the existence of miR399* with the high abundance in miR399d* and miR399f*, predominantly in the roots (Fig. 8A), which was consistent with the high level of their primary transcripts and of mature miR399d and miR399f .
We and others have previously shown that miR399 can move from shoots to roots and function as a longdistance signal in response to 2Pi Pant et al., 2008). We found that like miR399f, miR399f* could be detected in the wild-type rootstock grafted with the miR399f-overexpressing scions (Fig. 8B). This result indicates that both miR399f and miR399f* can move across the graft junction from scions to rootstocks. We hypothesize that miR399* may have certain physiological relevance, possibly functioning in concert with miR399.
There may be several possible functions of miR399*. First, miR399* may be loaded into the RNA-induced silencing complex to target different genes. Many miRNA* species from Drosophila were reported to associate with the Argonaute protein and have inhibitory activity (Okamura et al., 2008). Several genes were predicted to be the putative targets of miR399* species; however, we were unable to validate them by RLM 5#-RACE analysis. Thus, whether these miR399* species possess gene-silencing activity remains uncertain. Second, given the ability of miR399 and miR399* to move across the graft junction, we suspect that miR399* may assist the long-distance movement of miR399 by forming a miR399/miR399* duplex in the phloem stream. Several miRNA* species were also detected in the phloem saps of Brassica napus Pant et al., 2009). Although the authors concluded that they remain in the single-stranded form, based on the in vitro RNase digestion assay, they could not exclude the possibility of separation from a double-stranded RNA duplex during the RNA extraction. An in vivo nondestructive analysis needs to be developed to differentiate the double-stranded small RNA molecules from the single-stranded molecules. Alternatively, miR399* may regulate or buffer the targeting efficiency of miR399 on PHO2 mRNA via the formation of an RNA duplex.

CONCLUSION
In summary, we have presented an extensive survey of the small RNAs showing differential expression in response to 2Pi. These Pi-responsive small RNAs and their target genes are likely involved in the development or regulation of adaptive responses to 2Pi. Interestingly, while miR399, miR778, miR827, and miR2111 responded specifically to 2Pi, miR156, miR169, miR395, and miR398 also responded to other nutrient stresses. These responses involve upregulation or down-regulation of distinct miRNAs, which coordinate and balance the demands of different nutrients. We hypothesize that metabolic adjustment, carbon assimilation, or oxidative stresses may have an effect on the cross talk among different nutrients, as we observed here. Although there is still much more to learn about the biological significances of these Pi-responsive small RNAs, this work has opened a new avenue for the functional study of small RNAmediated gene regulation in 2Pi responses.  Borevitz et al., 2000;Alonso et al., 2003) were grown hydroponically as described Chiou et al., 2006). Seventeen-day-old seedlings were subjected to 2Pi treatment by transferring them from +Pi (250 mM KH 2 PO 4 ) into Pi-free half-strength Hoagland solution. Roots and shoots were collected separately after 5 d of treatment. Control samples were maintained in +Pi and collected at the same time. For treatments with different nutrient deficiency, the specific nutrient element was withdrawn from the nutrient solution for 5 d. To compensate for the ion originally associated with the omitted nutrient element, it was supplied in the form of chlorite salt at equivalent concentration except in the 2Cu and 2Fe treatments, in which they were omitted from the solution.

RNA Preparation and Detection and High-Throughput Small RNA Sequencing
Total RNA was isolated using Trizol reagent (Invitrogen). Four RNA samples from +Pi or 2Pi roots or shoots of wild-type plants were used to generate small RNA libraries by the Solexa sequencing technology (Illumina). Small RNA gel-blot analysis was carried out as described  with modification according to Pall et al. (2007). Quantitative RT-PCR analysis was performed as described . The sequences of probes and PCR primers are listed in Supplemental Table S5.

Analysis of Sequencing Data
The identical adapter-trimmed reads were grouped as unique sequences with associated counts of the individual sequence reads. Each unique sequence was aligned to the Arabidopsis genome (The Arabidopsis Information Resource 8; http://www.arabidopsis.org/), and only perfectly mapped sequences were retained and analyzed further (for details, see Supplemental Materials and Methods S1).
We defined a sequence cluster as a group of same-strand sequences that were mapped to similar genomic locations, and the union of the mapped regions of the individual sequences of the cluster can form a consecutive nongapped region. For every cluster, the number of sequence reads from each of the four samples was counted separately and used for differential expression analysis.
We used two approaches to perform the differential expression analysis, a cluster-based approach for finding novel differentially expressed gene models and a gene-model-based approach for analyzing known gene models, especially for miRNA genes and ta-siRNA genes. All sequence clusters or genes of interest were compared between two libraries for obtaining differentially expressed ones using a statistics test.
For each of the differentially expressed sequence clusters that were not annotated, the sequence including the cluster union and 200-bp flanking genomic sequence extending from each side was extracted and folded by RNAfold (Hofacker, 2003). A sequence cluster was considered to be a novel miRNA if it met the following four criteria: (1) the structure forms a hairpin; (2) the sequence cluster is located in the duplex region of the hairpin structure; (3) there is another sequence cluster located on the other strand of the duplex; and (4) the most abundant sequences from each of the two clusters can form a duplex with two-nucleotide overhang at the 3# end. The targets for the novel miRNAs were predicted using methods described elsewhere (Allen et al., 2005). Validation of the target gene was carried out by RLM 5#-RACE analysis using the GeneRacer kit (Invitrogen).

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. Multiple processing of the miR778 precursor.
Supplemental Figure S2. Expression of the target genes of miR395 and miR827.
Supplemental Table S1. Statistics of small RNA sequences.
Supplemental Table S2. Categories of small RNA origins.
Supplemental Table S5. Sequences of probes and primers.
Supplemental Materials and Methods S1.