Genome-wide identification of soybean microRNAs and their targets reveals their organ-specificity and responses to phosphate starvation

Abstract

Background

Phosphorus (P) plays important roles in plant growth and development. MicroRNAs involved in P signaling have been identified in Arabidopsis and rice, but P-responsive microRNAs and their targets in soybean leaves and roots are poorly understood.

Results

Using high-throughput sequencing-by-synthesis (SBS) technology, we sequenced four small RNA libraries from leaves and roots grown under phosphate (Pi)-sufficient (+Pi) and Pi-depleted (-Pi) conditions, respectively, and one RNA degradome library from Pi-depleted roots at the genome-wide level. Each library generated ∼21.45−28.63 million short sequences, resulting in ∼20.56−27.08 million clean reads. From those sequences, a total of 126 miRNAs, with 154 gene targets were computationally predicted. This included 92 new miRNA candidates with 20-23 nucleotides that were perfectly matched to the Glycine max genome 1.0, 70 of which belong to 21 miRNA families and the remaining 22 miRNA unassigned into any existing miRNA family in miRBase 18.0. Under both +Pi and -Pi conditions, 112 of 126 total miRNAs (89%) were expressed in both leaves and roots. Under +Pi conditions, 12 leaf- and 2 root-specific miRNAs were detected; while under -Pi conditions, 10 leaf- and 4 root-specific miRNAs were identified. Collectively, 25 miRNAs were induced and 11 miRNAs were repressed by Pi starvation in soybean. Then, stem-loop real-time PCR confirmed expression of four selected P-responsive miRNAs, and RLM-5’ RACE confirmed that a PHO2 and GmPT5, a kelch-domain containing protein, and a Myb transcription factor, respectively are targets of miR399, miR2111, and miR159e-3p. Finally, P-responsive cis-elements in the promoter regions of soybean miRNA genes were analyzed at the genome-wide scale.

Conclusions

Leaf- and root-specific miRNAs, and P-responsive miRNAs in soybean were identified genome-wide. A total of 154 target genes of miRNAs were predicted via degradome sequencing and computational analyses. The targets of miR399, miR2111, and miR159e-3p were confirmed. Taken together, our study implies the important roles of miRNAs in P signaling and provides clues for deciphering the functions for microRNA/target modules in soybean.

Background

Non-coding small RNAs can be grouped into small interfering RNAs (siRNAs) and microRNAs (miRNA) in plants[1]. SiRNAs consist of trans-acting siRNA(ta-siRNA), natural antisense transcript-derived siRNA(nat-siRNA), and repeat associated siRNA (ra-siRNA)[1]. Increasing evidences verify that 20–24 nucleotide-long miRNAs play important roles in growth, development, and stress adaptations in planta via modulating gene activity[1, 2]. In the plant cell, primary miRNA (pri-miRNA) is transcribed by RNA polymerase II (PolII)[3], then cut to precursor of miRNA (pre-miRNA) containing the distinctive stem-loop structure by Dicer-like (DCL), and finally pre-miRNA is processed to form miRNA/miRNA* duplexes, with miRNA* ultimately degraded[1, 4]. Hua Enhancer 1 (HEN1) methylates the miRNA/miRNA* duplex on the 3’terminal nucleotide of each strand to protect it from degradation by other small exonucleases. The methylated duplex is transported into the cytoplasm with the help of HASTY (HST)[1] and recruited by ARGONAUTE (AGO) proteins to suppress translation or to cleave the target transcripts mostly in the coding region[1].

Like DCLs, AGOs are highly conserved proteins across plants and animals. In Arabidopsis, ten genes encode AGO, which are grouped into 4 clades[1, 5]. AGO1 is involved in most miRNA biogenesis. AGO4 and AGO6 are responsible for ra-siRNA production and control DNA methylation. AGO7 participates in ta-siRNA formation[5]. Recent studies showed that the 5’ terminal nucleotide appears to determine the fate of miRNAs. For instance, AGO1 preferentially binds miRNAs with a 5’ terminal uridine (U), and AGO2 and AGO4 recruit small RNAs with 5’ terminal adenosine (A), AGO5 binds small RNAs that terminate with cytosine (C) at the 5’ terminus[5].

Recently, a number of miRNAs have been documented to be involved in nutrient signaling. Phosphorus (P) deficiency specifically induces miR399 in Arabidopsis and rice[6–8]. MiR827 induced by P limitation negatively regulates the transcript level of NITROGEN LIMITATION ADAPTATION (NLA) in Arabidopsis. Plus, osa-miR827 expression is strongly enhanced by phosphate (Pi) starvation in both roots and leaves. In situ hybridization indicates that osa-miR827 is expressed in mesophyll, epidermis and ground tissues of roots. Moreover osa-miR827 relays P signaling via negatively regulating the target genes of OsSPX-MFS1 and OsSPX-MFS2[9]. More ambiguous results pertaining to miRNA involvement in P deprivation was obtained with the observation through RNA ligase mediated 5’rapid amplification of cDNA ends (RLM-5’ RACE) that miR2111 cleaves At3g27150, which encodes a kelch domain-containing F-box protein, but Pi starvation induced the levels of both miR2111 and At3g27150[10]. One last set of results in regards to P nutrition shows that miR778, miR398, miR169 and miR408 are also responsive to P limitation[10, 11]. In regards to other nutrients, miR167 and miR393 were found to regulate root development in response to nitrogen (N)[12, 13]. Sulphur (S) starvation stimulated miR395, which targets plastidic ATP sulfurylase (APS), and the high-affinity sulfate transporter SULTR2;1[6, 12, 14]. Copper (Cu) limitation stimulates the level of miR398 and miR402[6, 15, 16]. Once expressed, miRNAs can be transported in phloem but not xylem vessels[12]. The levels of miR395, miR398 and miR399 were strongly augmented in response to S, Cu or Pi starvation in phloem. Under iron (Fe) deficiency, while the levels of miR399 and miR2111 were decreased in phloem, and even undetectable in roots and leaves[12, 16], indicating that Fe limitation alters P homeostasis.

The first study on soybean miRNAs in 2008 identified 35 novel miRNA families[17], and then sixty-nine miRNAs grouped into 33 families and their targets were determined through computational analysis[18]. After that, eighty-seven novel miRNAs were further isolated from soybean roots, seeds, flowers and nodules[19]. Twenty-six new miRNAs were identified via small RNA and degradome-associated deep sequencing in soybean seeds[20]. In addition, mis-expressions of miR482, miR1512, and miR1515 in soybean increased nodulation[21].

Low P availability resulting from its low mobility in soils is a common limiting factor for soybean yield[22]. Although miR399, the classic P-responsive miRNA, was predicted to exist in soybean[18], no experiments to verify the existence of this or other P-responsive mRNAs have been reported in soybean. To identify P-responsive miRNAs, along with, leaf- and root-specific miRNAs in soybean, in this study we sequenced four small RNA libraries from soybean leaves and roots grown under Pi-sufficient (+Pi) and Pi-depleted (-Pi) conditions, respectively. From those libraries, new conserved, less-conserved, and novel miRNAs as well as P-responsive miRNAs in soybean leaves and roots were identified genome-wide. Furthermore the targets of soybean miRNAs were determined via degradome sequencing from -Pi-roots combined with computational analyses. Finally, the existence of four P-responsive miRNAs (miR399, nov_6, nov_9, and nov_10) were verified via stem-loop real time (RT) PCR, and the targets of miR399, miR2111, and miR159e-3p were confirmed via RLM 5’ RACE. With results in hand, the possible functions of miRNA/target modules are discussed.

Results

Deep sequencing of small RNAs in soybean leaves and roots

Relative to control plants, long-term (7 days or 14 days) Pi starvation increased the ratio of roots to shoots (Additional file1A), and decreased the concentration of soluble phosphate (SPi) in leaves and roots (Additional file1B). Pi limitation globally induces the expression of many genes, including AtIPS1 (At3g09922) a well known non-coding Pi-starvation responsive gene, and AtPLDZ2 (At3g05630), which is involved in the conversion of phospholipid to glycolipid and free P[7]. Based on Blast searches of these genes in Phytozome (http://www.phytozome.org) and TIGR (http://www.tigr.org), the sequences of GmIPS1 and GmPLDZ2 were retrieved. The transcript levels of GmPLDZ2 (Additional file1C) and GmIPS1 (Additional file1D) were induced by long-term Pi deprivation. These data verified that the Pi starvation in soybean was affected at both physiological and molecular levels.

Additional file1 indicated that soybean was subjected to Pi starvation at 6 h, 12 h, 7 d and 14 d treatments. Accordingly, the 6 h and 12 h Pi starvation were treated as short-term P limitation, and 7 d and 14 d as long-term P stress. To augment the chance of finding more miRNAs in a single sequencing run, total RNA was extracted from leaves and roots of both +Pi and -Pi plants at 6 h, 12 h, 7 d, and 14 d, and then equal amounts of total RNA from each time point and treatment were pooled and used to construct four small RNA libraries labeled as leaf+Pi (HPL), root+Pi (HPR), leaf-Pi (LPL) and root-Pi (LPR).

Overall, more than 20 million raw reads from the four small RNA libraries were obtained. Clean reads were produced by excluding reads smaller than 18 nt and adaptors. The percentage of clean reads to total reads ranged from 96.92% to 99.59% (Table1). More than 80% and 75% of total RNA reads from the two leaf libraries and two root libraries were mapped to soybean genome, respectively (Table1). Here, unique small RNAs were defined as small RNA species with unique sequence. For unique small RNAs, more than 80% of reads in HPL and LPL libraries, and more than 55% of reads in HPR and LPR libraries could be mapped to the genome, respectively (Table1).

Table 1 Statistics of four small RNA libraries from soybean

Additional file2 summarized the origin profiles of total small RNAs. Fifteen and 14 percent of clean small RNA reads were mapped to known miRNA genes in HPR and LPR, respectively, but over 48% of small RNA reads from HPL and LPL libraries were mapped to miRNA genes. Additional file3 demonstrated that the percentage of clean small RNAs with 20-24 nt was 85.16%, 89.90%, 52.22%, and 46.05% in HPL, LPL, HPR, and LPR, respectively, indicating good quality for the four small RNA libraries. In addition, more than 20% of the clean small RNA reads were assigned to unannotated regions of the Glycine max 1.0 genome, supporting the notion that most small RNAs originate from intergenic regions.

Small RNAs containing 21, 22, or 24 nucleotides had high abundance relative to those of other lengths in all four libraries, while those containing 20 nucleotides were high in leaves (Figure1A; Additional file3). The 39-45% abundance of 21 nt small RNAs in HPL and LPL library was dominant over 22 and 24-nt small RNAs (Figure1A). When only unique reads were considered, 24-nt small RNAs were most abundant, with 21 nt and 22 nt small RNA following in all libraries (Figure1B). Overall, these data are consistent with previous studies[10, 20].

Figure 1

Size distribution of small RNAs through Solexa sequencing. Small RNA size distribution from the libraries of leaves and roots grown under +Pi (HP) and -Pi (LP) conditions. The size distribution was plotted versus frequency (%) of small RNA read length relative to total small RNA reads (A) or unique small RNA reads (B).

The nucleotide preference of small RNAs

As described above (see Background), AGO proteins preferentially recruit miRNAs based upon the 5’ terminal nucleotide[5]. In the sequences obtained from soybean, the representation of the nucleotides in the 5’ position varied among miRNA lengths (Additional file4).

Over 90% of the first nucleotide in 18 nt small RNAs in the two root libraries was C. The percentage of guanosine (G) in the first base of 19 nt small RNAs was the highest in the LPR library, whereas that of A was the highest in the HPR library. Interestingly, greater than 50% of U in the first position of 20 nt, 21 nt, and 22 nt-long small RNAs in all four libraries. Over 60% and 90% of 23 nt small RNAs were 5’-terminated with U in the two root libraries, but more than 90% of 23 nt small RNAs were tagged with G in the two leaf libraries. The percentage of A and U in the first position of 24 nt small RNAs were higher than that of the other nucleotides in leaf libraries, and the percentage of G in the first position of 24 nt small RNAs was higher than other nucleotides in the two root small RNA libraries. The 25 nt small RNAs with 5’ U were dominant in HPL, LPL, and HPR libraries (Additional file4). These results indicate that P availability generally played little roles in the first base bias of small RNAs in soybean.

Identification of conserved and new miRNAs

After excluding the small RNA reads which can be mapped to protein-coding and structural RNA-coding regions, candidate miRNAs were predicted based on published methods[23]. To determine bona fide miRNAs, several strict criteria to scrutinize mature miRNAs were implemented, including: (1) the abundance of putative miRNAs was at least 100 transcripts per million (TPM) in any one of the four libraries; (2) the reads for 5p and 3p strand could be detected and the ratio of 5p/3p or 3p/5p was higher than 0.9; (3) the length of small RNAs was 20 to 24 nt; (4) the minimum folding free energy of the precursor of miRNA was lower than -37 KJ/mol. (5) the minimal free energy index (MFEIs) was higher than 0.85; (6) RNAfold predicted a hairpin secondary structure for pre-miRNA; (7) features of real pre-miRNA, as tested in the online service, plantMIRNAPred[24], verified the candidate as pre-miRNA. Subsequently, candidate mature miRNAs were Blast searched against the soybean miRNAs deposited in miRBase 18.0 (http://www.mirbase.org) and new conserved, less-conserved (fabaceae-specific) and soybean novel microRNAs (soybean-specific) were determined.

As shown in Table2, seventeen soybean miRNA families have been found in other plant families, and six have only been noted in fabaceae. Mature miRNAs ranged from 20 to 23 nt, and 11 out of 23 miRNA families were 21 nt-long. Most of the first nucleotide in all soybean miRNAs were U. More C was present in the 19th nucleotide than any other nucleotide (Table2). This was consistent with earlier studies on soybean and cotton[18]. At the time of conducting this research, more than 240 soybean miRNAs, representing 65 families were curated in miRBase. A total of 126 soybean miRNAs in the current work were identified, 34 of which were annotated in miRBase 18.0. Eighty-nine of the 126 miRNAs, grouped into 21 families were conserved (Table3), and 25 were less conserved miRNAs that were only be found in legumes, namely legume-specific (Table4). The remaining 12 were soybean-specific (novel) miRNAs (Table5).

Table 2 Conserved and less-conserved miRNAs families detected in four small RNA libraries from soybean

Table 3 Conserved soybean miRNAs in four small RNA libraries

Table 4 Less-conserved soybean miRNAs in four small RNA libraries

Table 5 Novel soybean miRNAs in four small RNA libraries

Importantly, ninety-two soybean miRNAs have been identified in this study for the first time (Tables3,4 and5). Sixty-two of these are conserved miRNAs and grouped into 14 families (Table3). Eight are less-conserved and classified into 5 families (Table4). Ten of the miRNAs identified in this study can be found in other plant species, but have not yet been classified into any miRNA families in miRBase. Finally, 12 miRNAs identified herein are as yet soybean-specific novel miRNAs (Table5). More importantly, in regards to plant nutrition, miR399a, miR399b, miR399c, miR399d, and miR399e were detected in the present small libraries. Moreover miR399a, miR399b, and miR399e were localized to chromosome 5, while miR399c and miR399d were localized to chromosome 8 (Table3).

Soybean leaf- and root-specific miRNAs

To determine leaf- and root-specific miRNAs, the abundance of each mature miRNA was compared between leaves and roots. If a miRNA was not detectable in leaves or roots, then it was considered as specific to the other tissue. Under control (+Pi) conditions, ten leaf-specific, and four root-specific miRNAs were found (Figure2A). Under Pi-depleted conditions, there are twelve leaf-specific miRNAs, and two root-specific miRNAs (Figure2B). The two root-specific miRNAs in Pi-depleted plants are novel to soybean (gma-nov_2 and gma-nov_10).

Figure 2

Unique and overlapping soybean miRNAs in leaves and roots under +Pi and -Pi conditions. A: leaf-specific and root-specific miRNAs in +Pi conditions; B:leaf-specific and root-specific miRNAs under -Pi conditions.

A number of miRNAs are notable for their expression patterns. This included most members of miR156, 164 and 167 families, along with 12 individual miRNAs (miR168, miR172b-3p, miR2118a, miR2118b, miR408c, miR1507a, miR1508d, miR1508e, miR1509a, miR1510b-5p, miR1510c, and miR1511) that were found in high abundance (>1000 TPM) in one or both of the HPL or LPL treatments (Tables3,4 and5). Others appeared to be constitutively expressed in leaves and roots, including, for example, members of the conserved miR160, miR164, and miR2118 families (Table3). Within families, expression can vary significantly among family members. For example, expression levels of miR156 family members ranged from 46 to 66744 TPM in leaves, and from 0.01 (actually 0, only for normalization) to 22303 TPM in roots, and expression levels of miR166 family members ranged from from 48 to 16705 TPM in leaves, and from 12 to 20693 TPM in roots (Table3). Overall, a diverse range of responses was observed among treatments, tissues and family members, but variation between leaves and roots appeared to be more prevalent in less-conserved than in conserved miRNAs (Tables3,4 and5).

Identification of P-responsive miRNAs

Soybean P-responsive miRNAs were identified by calculating the log fold change in read counts, log ₂(-Pi/+Pi). If this value was greater than 1, the miRNAs were considered to be induced by P deficiency. Figure3A shows that the expression of 21 out of 126 (16.7%) mature miRNAs was stimulated in leaves by Pi starvation with fold changes ranging from 1 to 15.27. Moreover, in contrast to other -Pi-induced miRNAs, miR169q, miR396j, miR399e, and miR4416a were sharply induced in leaves by Pi limitation, with no expression detected in +Pi and expression over 10 in -Pi (Figure3A; Tables3 and4).

Figure 3

Fold change (log ₂ (-Pi/+Pi)) of miRNA expression analysis using transcript per million (TPM) transformed data from P treatment libraries. A: miRNAs from leaves; B: miRNAs from roots. Black and white bar indicates -Pi-induced and -Pi-repressed miRNAs, respectively.

In roots, twelve out of 126 (9.5%) miRNAs were stimulated by Pi starvation with log fold changes ranging from 1.03 to 14.18 (Figure3B). The levels of miR397a, miR482j, miR482k, and miRnov_2 in roots were induced from no detected expression to considerable expression in +Pi and -Pi (Figure3B; Tables3 and5). On the other hand, miR169r and miR3522a were significantly down-regulated by Pi depletion (Figure3B; Tables3 and4). Figure4A shows that thirty-six miRNAs were Pi-responsive in soybean genome-wide, with 26 and 18 being responsive in leaves and roots, respectively. Eight miRNAs were induced under -Pi conditions both in leaves and roots and none were repressed (Figure4B, C; Tables3,4 and5;). Perhaps more intriguingly, 5 miRNAs were tissue specific and had expression altered by P treatment. Those that were induced in -Pi relative to +Pi are the leaf specific miR396j and miRnov_7, and the root-specific miRnov_2 (Figure4B; Tables3 and5). Those that were repressed are the leaf-specific miR4376a and the root-specific miR169r (Figure4C; Tables3 and4). This tissue specificity and P responsiveness implies the tissue-specific roles of miRNAs in P signaling.

Figure 4

Unique and overlapping soybean P-responsive miRNAs in leaves and roots induced or repressed by -Pi treatment. A: Unique and overlapping soybean P-responsive miRNAs in leaves and roots; B: Unique and overlapping soybean miRNAs in leaves and roots up-regulated by -Pi treatment; C: Unique and overlapping soybean miRNAs in leaves and roots down-regulated by -Pi treatment.

The expression levels of mature miR399, miRnov_6, miRnov_9, and miRnov_10 were determined through stem-loop real-time PCR[25]. Because mature sequences of miR399a, 399b, 399c, and miR399d are 100% identical (Table3), it is impossible to distinguish the levels of the four mature miR399s with stem-loop PCR. The total abundance of miR399a/b/c/d was significantly induced by Pi starvation in both leaves (Figure5A) and roots (Figure5B) after 7 d of Pi starvation, but not at any other tested time point. The expression level of miRnov_6 in leaves (Figure5C) and roots (Figure5D) was decreased on day 14 of Pi starvation (P <0.05). The expression of miRnov_9 was repressed by P deficiency on day 14 in roots, but not in leaves (Figure5E and F). Interestingly, the level of miRnov_10 was repressed in leaves on day 14 (Figure5G) and induced in roots at 6 h (Figure5H), which might indicate that miRnov_10 plays an important role in local and systematic P signaling pathways.

Figure 5

Expression analysis of miRNA399a/b/c/d, miRnov_6, miRnov_9, miRnov_10 under Pi-starvation and others different nutrient deficiency. A-H: expression patterns of miRNA399a/b/c/d/ (A and B), miRnov_6 (C and D), miRnov_9 (E and F), and miRnov_10 (G and H) in leaves and roots in +Pi and -Pi; I-L: expression patterns of miRNA399a/b/c/d (I), miRnov_6 (J), miRnov_9 (K), and miRnov_10 (L) in leaves and roots under different nutrient deficiency conditions. Stars above bar indicate the significant difference.

To determine whether the responses of the above-mentioned miRNAs are specific to Pi stress, the responses to N, potassium (K), and S starvation were explored. Expression alterations of marker genes GmNRT2(Glyma13g39850), GmHAK1(Glyma19g45260), GmSult1 (Glyma06g111500) indicated that the tested soybean seedlings were really subjected to N, K, and S starvation, respectively (data not shown). The miRNA miR399a/b/c/d responded to Pi starvation in leaves and roots, but not to other nutrient deficiencies (Figure5I). Neither miRnov_6 nor miRnov_9 responded to any nutrient deficiencies at 6 h and 7 day either in leaves or roots (Figure5J, K), which is consistent with the requirement of 14-day Pi starvation for a response (Figure5C, D, E, F). Therefore, we could not determine whether miRnov_6 or miRnov_9 was specifically responsive to Pi. The expression of miRnov_10 in roots was induced by Pi starvation and K deficiency at 6 h (Figure5L), which was consistent with the previous test and indicative of a nonspecific response.

AtIPS1 attenuates miR399s activity via the formation of the three-nucleotide bulge in the highly complementary region where cleavage occurs[26]. Upon searching in soybean transcriptome (http://www.tigr.org), four GmIPS (GmIPS1-4) members were found. Among them, GmIPS1 nearly perfectly matched miR399a, b, c, d, and e over the center region, thus forming a three-nucleotide bulge (Additional file5), and thereby implying that soybean miR399 activity might be negatively modulated by GmIPS1 as in Arabidopsis.

Soybean root RNA degradome library sequencing

To determine the targets of soybean miRNAs, an LPR small RNA degradome library was constructed based on described methods[20]. A total of 28, 557,354 high quality reads containing more than 25 million clean reads were obtained from the root degradome library (Additional file6). In addition, slightly more 21- than 20-nt sequences were obtained (data not shown). After excluding reads smaller than 18 nt and adaptors, approximately 88% of clean reads were analyzed (Additional file6). Furthermore 83.93% of clean reads (21,185,857 out of 25,241,382) and 67.72% of unique reads (2,406,610 out of 3,553,894) were mapped to the soybean genome (Additional file6).

Target prediction of Glycine max conserved, less-conserved, and novel miRNAs

Pairfinder was employed to analyze the dedradome sequencing data as described in Methods. Degradome data showed that 51 genes were the targets of 19 conserved miRNAs (Additional file7), 10 genes were the targets of 4 less-conserved miRNAs (Additional file8), and 11 genes were the targets of 8 novel miRNAs (Additional file9). One gene, Glyma08g21620 was detected to be the target of both miR166g and miRnov_2. Hence, a total of 71 genes were determined to be the targets of 126 miRNAs by degradome sequencing.

Since some miRNA targets were not detected in degradome sequencing, psRNAtarget (http://plantgrn.noble.org/psRNATarget/) was employed as a complementary approach to predict miRNA target genes. In combination with degradome sequencing, a total of 154 genes were predicted or detected to be the target of 126 miRNAs. Putatively, 98 genes were attacked by 89 conserved miRNAs, 37 genes were the targets of 25 less-conserved miRNAs, and 20 genes were targeted by 12 novel miRNAs. Previous studies indicated that conserved miRNAs attach to targets in the CDS. Consistent with this, for the conserved and less-conserved miRNAs, 94.8% (182 out of 192) and 85.2% (40 out of 47) of cleavage sites were in CDS regions, respectively. For novel miRNAs, the percentage was 100%. Interestingly, the 28 transcripts targeted by the 12 novel soybean miRNAs included 19 transcripts for 7 miRNAs (miRnov_1a, miRnov_1b, miRnov_2, miRnov_3, miRnov_5a, miRnov_5b, miRnov_6, miRnov_7) in the root degradome library, and 9 transcript targets for four miRNAs (miRnov_4, miRnov_8, miRnov_9, miRnov_10) predicted through computational analysis (Additional file9).

As for the conserved miRNAs, the number of targets for conserved miRNA (2.31 targets per miRNA) was higher than that of less-conserved miRNA (1.53 targets per miRNA) (Additional files7 and8). Among the conserved miRNAs, five (miR162c, miR166h-3p, miR166j-3p, miR390b,and miR408b-5p) were identified to have only one target, while the rest have two or more targets, most notably miR156d, miR167e, miR167g, miR167j, miR167k, miR172b-3p, and miR172h-5p (Additional file7). In contrast, 15 out of 25 less-conserved miRNAs had only one target. However, the previously unreported less-conserved miR1508e and miR3508 each have seven targets (Additional file8). Glyma08g21610 and Glyma16g34300, which encode a AGO proteinendoing and a HD-ZIP protein respectively, were predicted target of miR166g and miR168 (Additional file7). A NF-YA gene Glyma10g10240 was the putative target of miR169c (Additional file7). Importantly, these predicted cleavage sites for the three targets are consistent with previous RLM-5’ RACE results[20, 27, 28], supporting the reliability of our computational analysis.

In general, miRNAs that are conserved across plants, such as miR156, miR164, miR167 and miR169, target transcription factors (TFs), whereas less-conserved miRNAs target fewer TFs (Additional file8). Overall, 54% (53 out of 98) of conserved miRNA target genes were TFs (Additional file7), while only 2.7% (1 out of 37) and 10.0% (2 out of 20) of less-conserved and novel miRNA target genes, respectively, were TFs (Additional files8 and9). This is in accordance with previous study[20].

In regards to nutrient stress, the targets of root P-responsive miRNAs and their cleavage sites are highlighted in Table6. A PHO2, two Pi transporter transcripts (Glyma10g04230 and Glyma14g36650), and an AP2 protein gene (Glyma01g13410) were predicted to be targets of miR399a, b, c, d, and e (Table6). In fact, miR2111 was detected in small RNA libraries in this study, but the level of it was too low to be filtered out based on the strict criteria. The target of miR2111 was predicted. As in Arabidopsis, a kelch repeat-containing F-box protein, Glyma16g01060, was the putative target of miR2111. The targets of the Pi starvation-induced miR169o, miR397a, miR408c, miR4416a, miRnov_2, and miRnov_7 were detected in the degradome library or computationally predicted, as were the targets of the P limitation-repressed miRNAs miR159e-3p, miR169r, miR2109, miR4376, miRnov_5a, miRnov_5b, and miRnov_6a (Table6). It must be pointed out that most of the targets listed in Additional files7,8 and9 need to be experimentally confirmed with RLM-5’ RACE or transit expression analysis.

Table 6 Target genes of P-responsive miRNAs in soybean

Determination of targets of miR399, miR2111, and miR159e-3p through RLM-5’ RACE

To test the accuracy of target gene cleavage site location results, RLM-5’ RACE was compared with predictions. The cleavage site of miR399 in PHO2 is predicted to occur at position 1050, and in Glyma10g04230 at position 407 (Additional file7). Figure6A shows the putative cleavage site of miR159e-3p in the target gene Glyma13g25716 at position 1340 as determined through degradome sequencing. In addition, Figure6B shows the cleavage sites of miR2111 in two target genes, Glyma19g25770 and Glyma16g06160, as predicted with WMD3 (http://www.weigelword.com). RLM-5’ RACE has been successfully employed to determine cleavage sites of miRNAs in soybean[20, 27]. In this study, RLM-5’RACE using RNA from Pi-depleted roots confirmed the cleaved fragments of GmPHO2 (Glyma13g31290) and GmPT5 (Glyma10g04230) mRNA predicted previously (Figure6C, D). The experimentally determined cleavage site of miR159e-3p in Glyma13g25716 (Figure6E) matched the predicted site, and the cleavage site in Glyma16g06160 predicted with WMD3 (Figure6F). These results indicate that degradome sequencing and computational predictions can reliably predict miRNA interactions with target genes.

Figure 6

Confirmation of the targets of miR399, miR2111, and miR159e-3p by RLM-5’ RACE. Degradome sequencing analysis indicating a cleavage site at position 1050 from the 5’ end of the transcript of Glyma13g25716.1 is targeted by miR159e-3p (A). Target prediction of miR2111 with WMD3 (http://www.weigelworld.com) showing Glyma19g25770 and Glyma16g06160 as putative targets of miR2111 (B). Cleavage sites from the 5’ end for Glyma13g31290 (PHO2-1) by miR399 (C), GmPT5 by miR399 (D), Glyma13g25716 by miR159e-3p (E), and Glyma16g06160 by miR2111 (F). The arrow indicates the cleavage position in the transcript of target genes, and the number above the arrow shows the sequenced clone numbers of PCR products.

Possible functions of soybean miRNAs’ targets

A total of 154 target genes were identified through degradome sequencing and computational predictions (Additional files7,8 and9). To better understand the biological functions of these genes in soybean, GO analysis[29] was employed to classify target genes based on their involved biological processes.

As shown in Additional file10, a total of 154 target genes are positively or negatively involved in many biological processes in soybean. For instance, target genes positively regulate the following processes: (1) nucleoside, nucleotide, and nucleic acid metabolic process (45/154, GO:0019219); (2) RNA metabolism (35/154, GO:0051252); (3) gene expression (49/154, GO:0010468); (4) macromolecule biosynthetic process (46/154, GO:0010556); (5) meristem development (10/154, GO:0048509). On the other hand, target genes negative regulate these processes: (i) seed development (31/154, GO:0048316); (ii)shoot development (25/154, GO:0048367);(iii) root development (19/154, GO:0048364); (iv) meristem initiation (10/154, GO:0010014). These results indicate the important roles of target genes in soybean in response to Pi starvation.

A more narrow focus on the targets of -Pi induced and repressed miRNAs returns functions that can be grouped into several categories based on the involved processes (Table6). These include: (i) protein synthesis or degradation; (ii) P uptake and transport; (iii) stress-related processes; (iv) cell division; and (v) ROS homeostasis. These results indicate that complex networks of P signaling are present in soybean. Figure7 outlines the possible functions of P-responsive miRNAs and their target genes.

Figure 7

Possible functional networks for P-depletion responsive miRNAs in soybean. Relationships between 25 P-depletion induced and 11 P-depletion repressed miRNAs and their target genes shown based upon putative physiological functions. Bold font and normal font indicate -Pi induced or repressed miRNAs, respectively.

cis-element analysis in the promoter of P-responsive miRNAs

AtPHR1 binds the promoter region of miR399 in Arabidopsis[30], indicating control of miRNA expression by TFs. cis-element analysis shows that TATA-box and TSS elements in the pre-miRNA upstream region can be found in over 90% of miRNA genes (Additional file11). A total of 270 TSSs and 230 TATA-boxes were found in the promoters of 126 soybean miRNA genes (Additional file11), with more TSSs and TATA-boxes found within 1.0 kb upstream of pre-miRNA than in the 1.0 to 2.0 kb upstream region (Additional file12).

Several well-defined Pi-responsive cis-elements were selected as references to identify P-responsive cis-elements potentially regulating the currently studied miRNAs[31]. A total of 377 Pi-responsive elements were detected for 126 miRNA genes, the average number of cis-elements was 3.02 (Additional file13). Among them, miR156w, miR156x, miR166g, miR168c and miR168d all contained 8 cis-elements, and miR166a-5p harbored nine. Genes for miR399a, miR399b, miR399d and miR399e had PHR1 binding sites, but miR399c only had a W-BOX binding site (Additional file13), indicating the expression of miR399a, b, d, and e might be regulated by PHR1 in soybean. The average frequency of PHR1, PHO-like, Pi-responsive, W-box, TATA-box, and TC elements in -Pi-induced miRNA promoter regions was higher than that in the -Pi-depressed miRNA promoter region, while the frequency of PHO and NIT2 elements was lower in the Pi-induced miRNA promoter regions than that in the -Pi- depressed miRNA promoter regions (Figure8A).

Figure 8

Frequency of P-responsive cis -elements in the promoter regions of P-responsive miRNAs. Frequency of cis-elements types (A), and average number of cis-elements per upstream region (B) from 25 up-regulated and 11 down-regulated miRNA genes differentially expressed in -Pi relative to +Pi.

Discussion

Deep sequencing of small RNAs and RNA degradome in soybean

Soybean is an important crop that provides oils and proteins to human and animals. Little is known about the involvement of miRNAs in soybean leaf and root development, and P signaling. Although high-throughput sequencing has been employed to reveal soybean miRNAs[20], there are likely many miRNAs to still be found. SBS has been successfully employed to find miRNAs genome-wide in rice[32] and Glycine max[20]. Figure1 shows a major peak at 21 nt in the total small RNAs from four small RNA libraries. A peak at 19 nt found in Arabidopsis root libraries[10] was not found in the current study. On the other hand, our data identify 24 nt unique small RNA species as dominant over other kinds of small RNAs in four small RNA libraries, which is consistent with previous studies[10, 20, 23]. One possible reason for the prevalence of 24 nt small RNAs is that most precursors of siRNAs are processed into 24 nt-long siRNAs by DCL3[1]. In regards to nutrient effects, it has been reported that P availability does not change total small RNA profiles in rice leaf libraries[23]. This is consistent with the results reported here (Figure1), but it does not address the fact that a number of individual miRNAs are differentially expressed among P treatments (Figures2,3 and4). Hsieh et al.[10] (2009) reported that more than 24% of small RNA reads mapped to tRNA in the Arabidopsis genome in two root libraries, while the percentage in leaf libraries is lower[10]. This high percentage of tRNA in root libraries was not observed in soybean (Additional file2). In conclusion, the profiles of total small RNA in plants are dependent on plant species, tissue type and environmental conditions.

Since 2008, degradome sequencing has been employed to screen Arabidopsis miRNA targets[33]. A total of 174 genes targeted by 87 unique miRNAs were identified in rice cultivar 93-11 from a young panicles degradome library[32]. Here, degradome sequencing detected 71 genes to be cleaved by conserved, less-conserved, and novel soybean miRNAs. RLM-5’ RACE has confirmed the targets of miR166g, miR168 and miR169c in previous work[20, 27, 28], as well as the targets of miR399, miR2111 and miR159e-3p in this study. Moreover, the cleavage sites were in accordance with predictions or degradome sequencing (Additional file7). This indicates that degradome analysis is a powerful tool to find soybean miRNA targets. The limitation of this study is that just one degradome small RNA library was produced from Pi-starved roots. MicroRNA-degraded fragments existing specifically in leaves or other organs above the roots were not sequenced. Nevertheless, using a combination of experimental and computational approaches, 95 new targets for conserved, less-conserved, and novel miRNAs were putatively identified here. Accordingly, identification of more targets of soybean miRNAs will shed fresh light on the functions of soybean miRNAs in the near future.

Discovery of conserved, less-conserved, and novel miRNAs in soybean

Presently, a total of 395 soybean mature miRNAs (362 precursors) were curated in miRBase 18.0. In this study, we adopted very strict criteria to determined candidate miRNAs as outlined above. Interestingly, in contrast to previous studies[20, 34, 35], 62 of the conserved miRNAs (Table3), 18 of the less-conserved miRNAs (Table4), and 12 of the novel soybean-specific miRNAs (Table5) were new discoveries in soybean. A possible explanation for these discoveries is the sampling of tissues from multiple time points and treatments. RNA was sampled from leaves and roots at different time points and P treatments, which increases the odds of finding previously unreported miRNAs in soybean.

Our results significantly improve our understanding of miRNA families. The first to mention is the miR156 family, which is big and conserved across plants, including 12 miR156 members in both Arabidopsis and rice, 11 members in maize, and 8 in Medicago (http://www.phytozome.org). Here, the soybean miR156 family members have been expanded from 15 to 29 (Table3). In addition, 11 miR166s were identified for the first time (Table3), indicating that miR166 is also a very big miRNA family in soybean. We note that as an ancient polyploidy descendent, the genome of Glycine max was duplicated two times, 59 and 13 million years ago[36]. These events likely led to significant increases in the two miRNA families listed above, as well as, potentially increasing membership in other families. Furthermore, miR1507, miR1509, and miR1510 were identified, and to date, have only been found in Medicago and soybean (Table2). In short, the known sizes of soybean miRNA families were significantly expanded through this study. However, we did not detect miR157, miR393,miR395, or miR398 as previously reported. Possible explanations are that (i) strict criteria were implemented in this study for identification of mature miRNAs; and (ii) the RNA samples were limited to low P-stressed leaves and roots, along with their control, so -S-induced miR395 and -Cu-induced miR398 were not likely to be detected.

miRNAs in soybean leaves and roots

MiRNAs play crucial roles in leaf and root development[1]. Figure2 shows 112 miRNAs detected both in leaves and roots, along with leaf- and root-specific miRNAs extracted from +Pi or Pi treatments. Comparable to the reported very high levels of miR156f and miR156h in soybean leaves[35] were the currently reported very high expression levels of miR156d, miR156h, and miR156o in soybean leaves (Table3). This implies crucial roles for all of these miRNAs in leaf development. The high expression of 11 miRNAs (gma-miR164, miR167, miR168b, miR319a, miR396a, miR482b, miR482b*, miR2118a, miR2118b, miR1508a, and miR1509a) in soybean leaves has been verified by microarray analysis, as were low expression levels of miR169a, miR390c, miR1507c, and miR1510a[35]. In this study, the abundance of several miR169 family members (miR169c, miR169p, miR169q, and miR169r) was low yet still detectable in soybean leaves (Table3), which indicates some roles for these miRNAs in leaves, and stands in contrast to the report that gma-miR169 is abundant in soybean shoot apical meristem (SAM) and undetectable in mature soybean leaves[35]. The higher expression of miR164 in mature soybean leaves than in SAM has also been reported[35], and consistent with this, the levels of miR164e, f, g, h, and i in leaves were much higher than in roots (Table3). Over-expression of rice miR172 results in loss of spikelet determinacy and flower abnormalities[37]. Interestingly, miR172b-3p, miR172h-5p, and miR172k were present at higher levels in leaves than in roots (Table3), conversely, miR172b-3p was still found in significant numbers in roots. Transcripts of miR396k in leaves were higher than that in roots. This miRNA targets the cytokinin (CTK) oxidase/dehydrogenase encoded by Glyma12g01390 (Additional file7). As CTK regulates meristem activity in leaves[38], miR396k appears to regulate soybean leaf meristem activity. One interesting study that this work brings to mind is the exploration of the roles of the 10 and 12 leaf-specific miRNAs under +Pi and -Pi conditions, respectively (Figure2) in soybean compound leaf development.

Many miRNAs have been reported to regulate root development[13, 39, 40]. The expression of miR160-resistant AUXIN RESPONSE FACTOR 17 (ARF17) leads to a shorter primary root and few lateral roots in Arabidopsis[39]. Ata-miR164 mediates lateral root development through attacking NAC1, and miR167 modulates adventitious rooting via targeting ARF6 and ARF8[40]. Auxin receptors TIR1, AFB1, AFB2, and AFB3 are targets of miR393[13]. Relative to other miR166s, gma-miR166g, miR166k, miR166h-3p and miR166j-3p were abundant in soybean roots (Table3). AUX/IAAs are cleaved by miR167s (Additional file7), and the levels of miR167e, miR167j, and miR167l are high in roots (over 400 TPM) (Table3). Under +Pi conditions, miR169q, miR169r, miR4416a and miR1512b were specifically expressed in roots (Tables3 and4), and under -Pi conditions, miRnov_2 was specifically induced in roots (Table5), indicating these miRNA are crucial regulators in root development and adaptations to P starvation. It was reported that five conserved miRNAs (miR159, miR162, miR166, miR390, and miR399) presented similar expression levels in root apexes and nodules, but miR169, miR171, miR393, and miR396 enriched in root tips[41]. However, few miRNAs were reported to control root development and nodulation in soybean, future studies should focus on the field. All together, our results demonstrate that miRNAs are intimately involved in many aspects of plant development, and much of these roles remain to be elucidated.

P-responsive miRNAs in soybean leaves and roots

In this study, we found 26 P-responsive miRNAs in leaves and 18 P-responsive miRNAs in roots (Figure4A). Recently, more P-responsive miRNAs have been found in Arabidopsis[10], common bean[42], and white lupin[43]. MicroRNA chip experiments showed that eight miRNAs (miR156/157, miR167, miR168, miR319, miR159, miR894, miR1507, and miR1509) were induced by Pi starvation in soybean leaves, and seven miRNAs (miR159, miR894, miR1507, miR1509, miR396, miR474, and miR482) were induced in soybean roots by low P[31]. In Arabidopsis roots, miR156a, b, c, d, e, and f are moderately induced by P deficiency, and the ata-miR156 family is also induced by -N, -K[10]. In this study, no miR156 was found to be induced by low P in leaves or roots (Table3). A possible reason is that RNA samples were pooled over multiple time points spanning short to long durations of stress. Hsieh et al.[10] (2009) extracted RNA from Arabidopsis roots and leaves treated for 7 days with low P[10].

In Arabidopsis, miR399s induced by Pi starvation are transported over long distances through phloem[11, 30]. The expression of miR399 has been predominantly found in vascular tissues, especially in root phloem companion cells in Arabidopsis[12]. MiR399 accumulated in Medicago and tobacco roots during arbuscular mycorhizal symbiosis[44], and targeted a putative phosphate transporter (Mendtr5g076920.1)[45]. Concordantly, soybean miR399s (miR399a, miR399b, miR399c, and miR399d) were induced by Pi starvation both in leaves and roots, but the levels in leaves are around 2.5 times higher than in roots (Table3). Whether they are only accumulated in leaves and later transported to roots via the phloem, or if they play active roles in shoot P transport remains an open question. Grafting experiments involving root stocks and shoots with divergent miR399 sequences or expression patterns might answer this question. The qRT-PCR data in the present study confirms that miR399a/b/c/d specifically responds to Pi deprivation (Figure5), while RLM-5’ RACE confirms that soybean miR399 cleaves PHO2 and GmPT5, as in Arabidopsis, rice, and common bean[42]. It is worth exploring the conservation and variation of miR399s and target genes between soybean and other plant species.

Soybean miR396k and miR397a were also added as -Pi induced miRNAs in this study (Table3). Arabidopsis miR398a, miR398b, miR398c, and miR408 were repressed in leaves under Pi-depleted conditions[10], while Cu depletion induced the expression of miR408, miR399, and miR2111 in Arabidopsis, rice and Brassica[12]. In this study, miR408c was stimulated by low P in soybean leaves (Table3, Figure3), indicating divergent roles for this family among plant species, or complex interactions between P and Cu signaling in plants that remains to be adequately outlined. In contrast to miR482 expression in Arabidopsis[10], five miR482s were found to be induced by -Pi treatment (Table3; Figure3). Additionally, miR778 and miR827 are reported to be specifically induced by -Pi in Arabidopsis[10], but neither was detected under -Pi conditions in soybean. To date, no miR827 has been reported in legumes (http://www.mirbase.org), which, together with the present results, indicates that the miR827 gene may have ceased to function or has dramatically diverged over the course of legume evolution.

Past studies suggested that only miR399 and miR395 are transported from shoots to roots via phloem but not xylem vessels[12, 16]. With miRNAs only being found in phloem, and not xylem vessels[16], then it might be possible to differentiate local and systematic miRNA mediated responses of soybean to low P availability by profiling miRNAs in phloem sap in future.

A total of 125 putative cis-elements in 24 soybean-Pi responsive miRNA genes were found and P-responsive motifs exist in the promoter regions of 54 non-phosphorus responsive miRNAs[31]. Interestingly, the total P-responsive element frequency in P-responsive miRNAs was higher than that in non-responsive miRNAs (Figure8B). Whether this difference in cis-element frequency is the direct reason for induced expression of P-responsive miRNA under P stress in soybean is, potentially a worthwhile and enlightening project.

Possible functions of miRNAs/target modules in soybean

Many miRNAs function in growth, development and stress adaptations through the regulation of TFs[2]. Consistent with this, a total of 51 TF genes were targeted by conserved (Additional file7) and less-conserved miRNAs (Additional file8), indicating crucial roles for these miRNAs in soybean.

Fifteen miR156s target SPL s in Arabidopsis[46]. Among these, SPL3, SPL4, and SPL 5 function in controlling flower time and phase change, while SPL9 and SPL15 play roles in leaf initiation. Over-expression of SPL3 and SPL9 stimulates flowering and over-expression of miR156 delays flowering via down-regulation of SPL activity[47]. In Glycine max, there are 45 genes encoding SPL proteins (http://www.phytozome.org). In this study, 16 SPL transcripts are targets for 14 miR156 family members (Additional file7). This data also reveals that 6 soybean CUC-like genes are targets of miR164 (Additional file7), suggesting the functions of miR164/CUC modules in soybean leaf development and growth as they do in Arabidopsis[48].

Several HD-ZIP class III genes such as PHV, REV, and AtHB15 are negatively regulated by miR165/166[49]. Consistently, four HD-ZIP transcription factors were predicted to be cleaved by miR166 family miRNAs (Additional file7). NF-YAs, the targets of miR169, participate in N and drought stress responses[50]. Soybean 5 NF-YA transcripts are potentially targeted by miR169 members (Additional file7). In addition, miR169c and miR169r were down-regulated by -Pi in leaves and roots respectively, and miR169q was up-regulated in leaves by low P (Table6). Hence, miR169s might act as integrators of P and N signaling. Ata-miR169s are down-regulated by low P, consistently,some of which have PHR1 binding sites in the upstream region[10]. PHO-like binding sites and/or W-box elements were detected in promoter regions of soybean 169c, 169o, and 169q (Additional file12). BAM3 (At4g20270) encodes a receptor kinase-like protein that functions in shoot and flower meristem development[51]. Soybean miR390b, d, e and f all target Glyma02g45010, which shows very high similarity with BAM3.

In contrast to miR399 activity in Arabidopsis and rice[10], soybean miR399 targets PHO2 in the 5’ UTR but not in the coding region (Figure6C). GmPT5 (Glyma10g04230) was reported to be responsible for P homeostasis in nodule development[52]. RLM-5’ RACE data verified predictions that miR399a/b/c/d/e cleaves GmPT5 (Additional file7, Figure6D), implying potential regulation of P nutrition in nodules by miR399. Moreover, GmIPS1 might act as a RNA mimic to attenuate miR399 activity as AtIPS1/At4 does in Arabidopsis[26] (Additional file5). Then RNA mimcry might be a useful tool to decipher other miRNA functions in soybean. Among those with unknown functions is soybean miR2111, which attacks a kelch repeat-domain containing F-box protein gene (Figure6F), as it does in Arabidopsis[10]. However, the function of miR2111 and its target in Arabidopsis is still unclear. Exploring the role of miR399 and miR2111 in soybean in the near future at genetic and biochemical levels promises to yield useful insights into how these miRNAs function in plants, and their effects on associated networks.

AtAGO1(At1g48410) is cleaved by miR168a and miR168b in Arabidopsis[53]. Glyma16g34300 and Glyma09g29720, homologues of AtAGO1 are also targets of miR168 in soybean (Additional file7). Pentatricopeptide repeat (PPR) proteins are predicted to be involved in RNA editing and metabolism in mitochondria and are essential for 5’end processing of transcripts[54]. MiR1508d targets one PPR transcript, while miR1508e targets 7 PPR transcripts (Additional file8). In Arabidopsis, three PPR genes, At1g06580, At1g62720, and At1g62670, were predicted to be the targets of miR161[55]. These results indicate that regulation of PPR by miRNAs is conserved across plants, though it may be accomplished by different miRNAs among plant species. In Arabidopsis, SGS3 is involved in trans-acting siRNAs generation, and thus participates in the post-transcriptional gene silencing and natural virus resistance[56]. Interestingly, miR2118a and miR2118b both target an X1-like transcription factor encoded by Glyma05g33260, which is a homologue of AtSGS3. OsBIRH1, a DEAD box RNA helicase, regulates defense responses[57]. Four DEAD helicase were predicted to be the targets of miRnov_6 (Table6). Taken together, these results suggest that RNA metabolism is also tightly regulated by miRNAs in soybean.

Accumulation of reactive oxygen species (ROS) is one strategy for plants to cope with early stages of abiotic stress. However, higher ROS levels will damage proteins, nucleic acids, and lipids. P, K and Cu deficiency, as well as, drought stress boost ROS levels[58]. Under drought conditions, miR408 is strongly induced in photosynthetic tissues in Medicago[58]. Table3 demonstrates that miR408c was highly expressed in leaves in +Pi and induced even further in -Pi. These results imply that miR408c might be an integrator of stresses. Superoxide dismutase (SOD), peroxidase (POD), glutathione S-transferase (GST), and cytochrome P450 genes are strongly induced by long-term Pi starvation[59], further supporting the notion that ROS participate in root responses to low P stress. Interestingly, POD, NADP oxidase, and cytochrome P450 genes were targets of miR1507a, miRnov_1, and miRnov_9 (Table6; Additional files8 and9).

Glyma17g05970.1 was targeted by miR4376-5p, and miR4376-5p was undetectable in roots (Additional file8; Table4). The homologue of Glyma17g05970.1 in Arabidopsis regulates root hair growth, trichome development, and organelle trafficking[60]. Loss of function of AtENO1 (At1g74030), an Arabidopsis phosphoenolpyruvate enolase, results in distorted trichomes and fewer root hairs[60]. An enolase,Glyma19g37520 was the target of gma-miRnov_7, which is also specifically expressed in leaves. In roots, the abundance of miR4416a is higher than in leaves, and its target is a nodulin-like protein gene, Glyma10g06650 (Additional file8). This implies that miR4416a is potentially involved in nodule development in interactions with rhizobial symbionts.

The cell division activity gradually decreases around the meristem, leading to the determinate growth of Arabidopsis[61]. P deprivation increases the expression of Cyclin D3[62]. The data herein revealed that two Cyclin D3 genes were putative targets of gma-miRnov_10 (Table6), and miRnov_10 was down-regulated in leaves stressed by -Pi (Figure3; Table5). Stunted growth of P stressed soybean might be mediated through the actions of miRnov_10.

Recently, 28 TFs in maize roots have been documented to be induced by P limitation, while 14 TFs are repressed[63]. Although putative P-responsive cis-elements were found in soybean P-responsive miRNA genes (Figure8), the TFs responding to Pi starvation, and controlling root or leaf development in soybean remain unknown. In this study, many TFs were predicted to be targets of soybean miRNAs (Additional files7,8 and9). Identifying and understanding the activities of those TFs will likely provide insight into how plant growth and development respond to P deficiency.

Conclusions

A total of 126 miRNAs were identified in soybean through deep sequencing, including 92 previously unidentified (Tables3,4 and5). Among these, leaf- and root-specific miRNAs were determined (Figure2), and P-responsive miRNAs in leaves and roots were identified (Figure4A). These 126 soybean miRNAs target 154 genes as revealed via degradome sequencing and computational predictions (Addititonal file7,8 and9). Use of qRT-PCR verified the expression of four P-responsive miRNAs and 5’ RACE confirmed targets of miR399, miR2111, and miR159e-3p. Finally, cis-element analysis indicates the existence of P-responsive motifs in the promoter region of soybean miRNA genes. Taken together, these findings provide useful information for plant scientists to decipher soybean miRNA functions and establish a framework for exploring P signaling networks regulated by miRNAs. More extensive analysis of these miRNAs across time and spatial scales, and transgenic studies will facilitate exploring their roles in P signaling and leaf or root development.

Methods

Growth of soybean

The seeds of the sequenced soybean (Glycine max L. Merrill cv.Williams 82) were germinated in paper pouch for 2 days at 28°C in darkness and then grown in light for further 2 days, the uniform seedlings were transplanted into full nutrient solution (pH=5.9),which contained 250 μ M KH₂PO₄ and grown 5 days; soybean seedlings with the first fully developed trifoliate leaves were transferred into phosphate (Pi)-sufficient (250 μ M KH₂PO₄, +Pi) or Pi-deplete (0 μ M KH₂PO₄, -Pi) nutrient solution,respectively.Leaves and roots were separately sampled at 0 h, 6 h, 12 h, 7 day and 14 day after treatment. For nitrogen deprivation, NH₄NO₃ was omitted and KCl was substituted for KNO₃. For potassium depletion, KNO₃ was omitted. For sulfate starvation, all${\text{SO}}_4^{2-}$ was substituted with Cl^-1. All soybean seedlings were cultured in 24-l plastic boxes containing nutrient solution as indicated. The nutrient solution was automatically aerated 15 min every 3 h, and was replaced with fresh solutions every two days. Soybean plants were grown in a green house with 16 h light cycles at the Root Biology Center of South China Agriculture University.

Determination of ratio of root to shoot and concentration of soluble phosphate (SPi)

The dry weight (DW) of roots and aerial parts (including shoots and leaves) were determined with standard methods. The fresh root and leaf samples at different time points after +Pi and -Pi treatments were weighted separately, rinsed in distilled water and dried, frozen, and ground in liquid nitrogen. Approximately 200 mg of grounded samples were suspended in 2 ml distilled water overnight, and centrifuged to pellet the cellular debris. Subsequently the supernatant was assayed for SPi using phosphomolybdate colorimetric assay[64].

Quantitative real-time PCR

Total RNA were extracted from roots and leaves with a miRcute^TM miRNA Isolation Kit (Tiangen,http://www.tiangen.com) according to the manufacturer’s protocol. RNA samples were treated with RNase-free DNase I (Invitrogen,http://www.invitrogen.com) to avoid amplification from genomic DNA. The first cDNA strand was synthesized from total RNA using the PrimerScript RT Enzyme (TaKaRa,http://www.takara-bio.com/). The housekeeping gene Glyma17g239000 (GmEF1a) was used as an endogenous control to normalize the samples[65]. Quantitative real-time PCR (qRT-PCR) was performed using SYBR Premix EX Taq^TM (TaKaRa). cDNA sequences and genomic sequences of GmPLDZ2, GmIPS1, GmNRT2, GmHAK1, and GmSult1 were downloaded from Phytozome, and specific primer pairs were designed with PerlPrimer[66] as listed in Additional file14. All reactions were run on a Rotor-Gene 3000 (Corbett Research, Australia). Reaction conditions for thermal cycling were based on standard methods.

Construction of soybean small RNA libraries and degradome library and sequencing

The samples after harvested were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction. After RNA isolation, equal amounts of total RNA from roots and leaves at 6 h, 12 h, 7 d, and 14 d of Pi-replete or Pi-depleted treatments were mixed. Fragments of 18-30 bases were purified from 10 μ g total RNA mixture using a Novex 15% TBE-Urea gel. The 5’ and 3’ adaptors (Illumina) were added to the ends of fragments. Reverse transcription PCR (RT-PCR) was performed using a RT-PCR kit (Invitrogen). PCR products were purified and quantified for Illumina sequencing with a Solexa sequencer in Shenzhen Huada Gene Sci-Tech Company (Shenzhen, China).

The soybean Pi-depleted roots RNA degradome library was constructed as previously described[20, 67, 68]. The total RNA was extracted as described above. In brief, poly(A) RNA was extracted from 200 μ g of total RNA using miRcute^TM miRNA isolation Kit (Tiangen). A 5’ RNA adapter containing a MmeI recognition site was ligated to the poly(A) RNA possessing a 5’-phosphate with T4 RNA ligase (Ambion). The ligation products were purified and amplified by 20 PCR cycles and gel-purified for SBS sequencing in the Shenzhen Huada Gene Sci-Tech Company.

Identification of miRNAs

Small RNA reads and degradome reads were both generated through Illumina Genome Analyzer II high- throughput sequencing. After sequencing, low quality reads and clip adapter sequences from raw data were removed via software developed by BGI (http://www.genomics.cn). Small RNA libraries were then constructed from, small RNAs ranging from 18-31 nt, and then mapped to the soybean genome using SOAP2[69]. Unique RNA sequences that perfectly matched the genome were subjected to subsequent analysis. The RNA reads showing sequences identical to known miRNAs from miRBase (http://www.mirbase.org) were selected as the miRNA dataset for soybean. Further sequences matching non-coding rRNA, tRNA, snRNA and snoRNA in the Rfam database were removed. Then the reads overlapping with exons of protein coding genes were excluded to avoid mRNA contamination. The remaining sequences were considered to be candidate miRNA for further bioinformatics analyses.

Because miRNA precursors have a hairpin structure, 150-200 nt of the sequence flanking the genomic sequences of small RNAs was extracted from Phytozome (http://www.phytozome.org). The MIREAP pipelinewas then used to analyze structural features to identify new miRNA candidates (https://sourceforge.net/projects/mireap/). The resulting structures, with minimal matched nucleotide pairs of miRNA and miRNA* exceeding 16 nt and with maximal size differences of miRNA and miRNA* up to 4 nt, were retained as new miRNA candidates. The filtered pre-miRNA sequences were folded again using MFOLD and checked manually[70].

Stem-loop quantitative real time PCR

Stem-loop specific reverse transcription was carried out according to the previously described Methods[25, 71]. Reverse transcription reactions were performed using total RNA from soybean roots, and leaves grown under different nutrient conditions as indicated above. The gma-miR156b was used as a reference miRNA gene to normalize samples as previuously described[71]. Stem-loop specific reverse transcription was basically carried out according to the methods described[71]. The reactions contained 1 μ g of total RNA, and each reaction was primed with a pool of 0.25 μ M 5 gene-specific stem-loop primers. The RNA and primers were mixed with RNase-free water up to 10 μ l and incubated at 70°C for 5 min followed removed to ice-cooling immediately. Then, 6 μ l 5RT-Buffer, 1 μ l 5 mM dNTP, 0.5 μ l RNA Inhibitor, and 1 μ l 200 U MML-V RT Enzyme (Promega) were added and supplemented up to a final volume 30 μ l with RNase-free water. Synthesis was performed at 42°C for 30 min on a Veriti Thermal Cycler (Applied Biosystem), and inactivation of the enzyme was performed at 85°C for 5 min. Samples were then held at 4°C. All cDNA samples were 50-fold diluted with RNase-free water before being used as templates in RT-qPCR analysis. All primers used in stem-loop RT-PCR are listed in Additional file14.

RLM-5’ RACE

To determine the cleavage sites of miRNA on target genes, RLM-5’RACE was employed. Total RNA was extracted from Pi-depleted soybean leaves and roots with the miRcute miRNA Isolation Kit (TIANGEN,http://www.tiangen.com) as described by the manufacturer. Then, extracted total RNAs were ligated with 5’RACE oligo adaptors, and the reverse transcription was carried out based on the GeneRacer kit (Invitrogen,http://www.invitrogen.com). The first found PCR was carried out with 5’RACE general PCR and outer gene-specific PCR, and the next PCR was run using the diluted initial PCR reaction and inner 5’ RACE and gene-specific primers (GSP). PCR products were cloned and sequenced according to standard methods. Additional file14 lists 5’ RACE and gene-specific PCR primers.

Identification of target genes for miRNAs

After excluding low quality reads, reads with 5’ primer contaminants, reads without 3’ primer, reads without the insert tag, and reads shorter than 18 nt were removed. The remaining 20-21 nt-long reads with high quality were collected for subsequent analyses. Degradome tags were analyzed for expression and distribution on the soybean genome using SOAP2[69]. Raw sequences were first normalized to reads per 10 million (RP10M), and identical degradome sequences with single base over 0.7 percentage in the clean reads was classified as polyN. Then, distinct reads that perfectly matched soybean cDNA sequences were further analyzed. The Pairfinder.pl script developed by the BGI Degradome group was used to align clean sequences to soybean known miRNAs from miRBase and miRNAs identified in this study. All alignments with scores up to 7 and no mismatches at the cleavage site (between the 10th and 11th nucleotide) were considered candidate targets.

Analysis of miRNA promoter and cis-acting elements

Based on the pre-miRNA sequences of soybean miRNAs identified and perfectly mapped to the soybean genome (http://www/phytozome.net/soybean) in this study, and not including pre-miRNA sequences curated in miRBase but not found in the present study, 2 kb sequences upstream of the pre-miRNAs were downloaded from Phytozome using the soybean genome browser. These sequences were used to predict transcription start sites (TSSs) and TATA-boxes, and analyze P-responsive cis-elements[31].

Statistical analysis of data

All data for ratio of root to shoot, soluble phosphorus, and relative expression level of genes and mature miRNA are from experiments with three biological replicates. All data were analyzed using Origin 7.5 (OriginLab Corporation, USA) for calculating means and SEs, and SAS 6.2 (SAS Institute, USA) for ANOVA analyses.