Novel MiRNA and PhasiRNA Biogenesis Networks in Soybean Roots from Two Sister Lines That Are Resistant and Susceptible to SCN Race 4

The soybean cyst nematode (SCN), Heterodera glycines, is the most devastating pathogen of soybean worldwide. SiRNAs (small interfere RNAs) have been proven to induce the silencing of cyst nematode genes. However, whether small RNAs from soybean root have evolved a similar mechanism against SCN is unknown. Two genetically related soybean sister lines (ZP03-5373 and ZP03-5413), which are resistant and susceptible, respectively, to SCN race 4 infection were selected for small RNA deep sequencing to identify small RNAs targeted to SCN. We identified 71 less-conserved miRNAs-miRNAs* counterparts belonging to 32 families derived from 91 loci, and 88 novel soybean-specific miRNAs with distinct expression patterns. The identified miRNAs targeted 42 genes representing a wide range of enzymatic and regulatory activities. Roots of soybean conserved one TAS (Trans-acting siRNA) gene family with a similar but unique trans-acting small interfering RNA (tasiRNA) biogenesis profile. In addition, we found that six miRNAs (gma-miR393, 1507, 1510, 1515, 171, 2118) guide targets to produce secondary phasiRNAs (phased, secondary, small interfering RNAs) in soybean root. Multiple targets of these phasiRNAs were predicted and detected. Importantly, we also found that the expression of 34 miRNAs differed significantly between the two lines. Seven ZP03-5373-specific miRNAs were differentially expressed after SCN infection. Forty-four transcripts from SCN were predicted to be potential targets of ZP03-5373-specific differential miRNAs. These findings suggest that miRNAs play an important role in the soybean response to SCN.


Introduction
Soybean (Glycine max) is an agronomically important crop that is rich in human dietary protein.The soybean cyst nematode (SCN), Heterodera glycines Ichinohe, is an obligate sedentary endoparasite that causes extensive damage to soybean worldwide and accounts for over one billion dollars of crop loss annually in the US [1].These obligate parasites start their life cycles as nonfeeding, mobile infective second-stage juveniles (J2) in soil that are able to locate and then penetrate into host roots [2].The J2 then initiate the formation of specialized feeding sites called syncytia, which function as metabolic sinks to nourish the nematodes.In susceptible cultivars, nematodes depend entirely on functional syncytia to acquire nutrients to develop into reproductive adult males or females.The J2 also penetrate roots of resistant cultivars and initiate syncytia.However, resistance soon manifests as degeneration of the young syncytia and failure of the nematode to develop further [3].Syncytium formation and maintenance are mediated through nematode signaling and are accompanied by changes in plant gene expression [4].Identification of host plant genes and nematode genes that change expression, and may therefore be involved in plant-nematode interactions, would increase our understanding of the molecular mechanisms involved in this complex interaction, which will lead to the development of durable crop protection strategies.With the recent discovery of gene expression control of parasitism proteins via siRNA molecules [5], and recent advances in genomics, small RNAs (sRNAs), which are involved in the molecular mechanism of the soybean-SCN system, are now the focus of much research.
Endogenous sRNAs are known to be important regulators of gene expression at the transcriptional and post-transcriptional levels.In plants they are divided into several classes: trans-acting siRNAs (tasiRNAs), heterochromatin-associated siRNAs, natural antisense siRNAs (nat-siRNAs) and miRNAs [6].These classes of non-coding RNAs are distinguished by their biogenesis pathways and the types of genomic loci from which they arise [7].TasiRNA biogenesis from TAS loci depends on the miRNA-directed cleavage of their transcripts [8,9]; indeed, three tasiRNA pathways have been characterized in Arabidopsis [8,10].Although miRNAs constitute only a small fraction of the sRNA population [11,12], miRNA-guided post-transcriptional gene regulation is one of the most conserved and well-characterized gene regulatory mechanisms [11,13,14].There is increasing evidence that miRNAs negatively regulate their target genes, which function in a wide range of biological processes, including organogenesis, signal transduction and stress responses [15,16].
Plant miRNAs are generated from hairpin-structured noncoding transcripts and processed by Dicer such as DCL1 (DICER-LIKE 1), which cleaves a short (21-bp) duplex from the stem region [17].The duplex is incorporated into an AGO complex and the miRNA* strand is subsequently degraded.The mature miRNA strand guides the AGO complex (RNA-induced silencing complex, RISC) to either protein-coding RNAs, which are cleaved by AGO at a specific position [18], or translational arrest [19].Due to their evolutionary conservation, miRNAs have been found to exist in both plants [20,21] and animals [22][23][24].Conserved miRNA molecules can also be found in ferns, mosses and fungi [11].To date, many miRNAs have been identified and deposited in miRBase V20.0 (http://www.mirbase.org/).Of these, 25,141 are mature miRNA products, from a total of 193 species.Comparative analysis indicates that some of the miRNA families are highly conserved among all plant species while others have diverged and evolved, generating abundant family-and speciesspecific miRNAs [11,25,26].These dynamic and evolving miRNAs could serve as a driving force for the selection of improved and novel traits in plants.
As non-conserved or species-specific miRNAs are often expressed at a lower level than conserved miRNAs, many species-specific miRNAs have not been identified in small-scale sequencing projects.However, high-throughput sequencing technologies allow identification of many species-specific miRNAs in several species [10,[27][28][29][30][31].Elucidating the function of these molecules requires effective approaches to identifying their targets.Recently, a new method called degradome sequencing, which combines high-throughput RNA sequencing with bioinformatic tools, has been used to screen for miRNA targets in Arabidopsis [32][33][34].Using degradome sequencing, many of the previously validated and predicted targets of miRNAs and tasiRNAs have been verified [32,33,35,36], indicating that this is an efficient strategy for identifying sRNA targets in plants on a large scale.
To determine if soybean has evolved sRNAs that repress the development and growth of SCN, and their potential targets, we selected the sister lines ZP03-5373 and ZP03-5413, which are resistant and susceptible, respectively, to SCN race 4 infection, and performed a comprehensive analysis of root miRNAs by deep sequencing, computational prediction and molecular approaches.Novel and conserved soybean miRNAs, tasiRNAs and phasiR-NAs, and their partial targets were identified.Small RNAs upregulated by SCN infection were identified and the molecular regulation mechanism was discussed.

sRNA population in soybean root
To investigate the role of soybean miRNAs in response to SCN infection, two genetically related soybean sister lines (ZP03-5373 and ZP03-5413) were subjected to deep sequencing.The sister lines shared the same parents and displayed a different resistance to SCN race 4. Line ZP03-5373 exhibited high resistance to SCN race 4, whereas ZP03-5413 was susceptible to race 4. Two sRNA libraries from the roots of the sister lines were constructed and sequenced using Illumina GAIIx.A total of 15,101,204 sRNA raw reads were generated.After removing adaptor sequences, filtering out low quality reads and cleaning up sequences derived from adaptor-adaptor ligation, 7,903,242 and 5,931,837 reads, respectively, were obtained.These sRNAs consisted of 4,979,640 unique sequences (Table S1), which were matched to the public soybean genomic database (Soybean Genome V9.0, http://www.phytozome.net/index.php)using the SOAP program, leading to 3,409,866 genome-matched unique reads.These reads were subjected to further analysis (Table S1).The 20-24-nt sRNAs constituted over 80% of the identified soybean sRNAs, and the 21nt class of sRNAs was the most abundant in both lines (Figure 1A).Notably, the expression of the unique 24-nt sRNAs was markedly higher than the 21-nt class in both lines (Figure 1B).

Conserved and less-conserved miRNA families and their expression in soybean root
The reads (3 million) that mapped perfectly to the soybean genome were subjected to miRNA identification.miRBase 20.0, which contains 555 soybean mature miRNAs, was searched for known soybean miRNAs.As a result, a total of 420 known soybean mature miRNAs were identified from the two libraries, of which 364 were sequenced in both libraries; 26 miRNAs were detected in only ZP03-5373 and 30 in only ZP03-5413 (Table S2).Expression levels of the known miRNAs, as reflected by normalized reads (reads per million genome-matched reads, RPM), varied substantially among families in both lines.The highest read abundance (31,416 RPM and 20,776 RPM) was detected for gma-miR159, and was 2-25-fold greater than other relatively abundant miRNA families, including gma-miR396, gma-156, miR168, and miR166, whose total abundance ranged from 1,000 to 15,000 RPM (Table S2).Gma-miR862a was expressed in ZP03-5413 but not in ZP03-5373 (Table S2).Substantial variation was observed for gma-miR393,      gma-miR398 and gma-miR399, whose abundance in ZP03-5373 was 10-fold greater than in ZP03-5413 (Table S2).
After excluding sRNA reads that perfectly matched known soybean miRNAs, the remaining 21 to 22 nt were subjected to rigorous secondary structural analysis of their precursors using RNAfold software (http://mfold.rna.albany.edu/).The minimum free energy (MFE) of the hairpin structure of the miRNA precursor was set to less than 225 kcal/mol.Those precursors with a canonical stem-loop structure were further analyzed by means of a series of stringent filter strategies to ensure that they met common criteria [37].Precursors carrying both the miRNA-5p and miRNA miRNA-3p were then selected for further analysis.A total of 258 new soybean miRNA candidates that were not previously reported, including miRNA-5p and miRNA-3p, were identified from the two libraries, of which 132 were sequenced in both libraries.Eighty-seven miRNAs were detected in only ZP03-5373 and 39 in only ZP03-5413.Novel miRNA candidates were further assigned to miRNA families using sequences similar to other known miRNAs in the miRBase database (two or fewer mismatches).These miRNAs or families had been reported previously in some plant species or families, but are not conserved among angiosperm and coniferophyta lineages [26].They were referred to as less-conserved miRNAs in this study.A total of 71 miRNAs-miRNAs* counterparts belonging to 32 families derived from 91 loci (Table 1 and Figure S1), that had previously been identified and reported in at least one plant species or family [11] were identified from the 258 miRNA candidates.A canonical predicted stem-loop structure could be identified in all 32 less-conserved miRNA families (Figure S1).Overall, all lessconserved miRNAs displayed lower expression levels than the conserved miRNAs, with the exception of gma-miR482C2, which was expressed at abundances of 4,000 RPM and 8,000 RPM in ZP03-5373 and ZP03-5413, respectively (Table 1).However, as with the conserved miRNAs, some of the less-conserved miRNAs were expressed differentially between ZP03-5373 and ZP03-5413.For example, ZP03-5413-biased expression was observed for gma-miR395C1, while ZP03-5373-biased expression was apparent for gma-miR393C1 and gma-miR2109C1 (Table 1).To validate the miRNA RPM data, we performed stem-stoop-based qRT-PCR analysis for selected miRNAs representing conserved, lessconserved and soybean-specific (discussed below) examples in the two lines.We found that while the qRT-PCR results for most of the miRNAs (miR1509a, miR1509b, miR2111 (Figure 2A) and miR395C1 (Figure 2B), miRC2, miRC6, miRC20 (Figure 2C), etc.) were reflective of the relative abundances of the sequenced RNAs in the two lines, others displayed varying degrees of divergence between the two analyses.For example, the miRC18 RPM value for ZP03-5373 was fourfold higher than for miRC10, while the abundances of miRC18 and miRC10 were in agreement, based on qRT-PCR results (Figure 2C).For miR482C2 the opposite pattern between qRT-PCR and miRNA sequencing was observed (Figure S2), which may have resulted from deep-sequencing deviation.

Prediction and validation of novel soybean-specific miRNAs
Because numerous species-specific miRNAs considered to be of a more recent evolutionary origin [11] have been identified in other species, soybean is likely to have evolved unique miRNAs.After excluding sRNA reads homologous to known miRNAs or families (two or fewer mismatches, miRBase 20), the remaining 21-22-nt sRNAs were subjected to rigorous secondary structural analysis of their precursors using the RNAfold software (http:// mfold.rna.albany.edu/).Those precursors with a canonical Cleavage site; b Category 0: .1 raw read at the position, abundance at position is equal to the maximum on the transcript and there is only one maximum on the transcript.Category 1: .1 raw read at the position, abundance at position is equal to the maximum on the transcript, and there is more than one maximum position on the transcript.Category 2: .1 raw read at the position, abundance at position is less than the maximum but higher than the median for the transcript.Category 3: .1 raw read at the position, abundance at position is equal to or less than the median for the transcript.Category 4: only one raw read at the position.P-value should not exceed 0.05.doi:10.1371/journal.pone.0110051.t003 stem-loop structure were further analyzed by means of a series of stringent filter strategies to ensure that they met the criteria established by the research community [37].A total of 74 miRNA family candidates derived from 88 loci (Table 2 and Figure S3) met the screening criteria, of which all had miRNA star (miRNA*) sequences identified from the same libraries.We termed these soybean-specific miRNAs.Of the 88 soybean-specific miRNAs, 75 belonged to the 21-nt class and 13 to the 22-nt class (Table 2 and Figure S2).In general, the soybean-specific miRNAs were less abundant than the conserved and less-conserved miRNAs in the two lines examined.For example, only gma-miRC20 displayed a read abundance above 600 RPM in ZP03-5373, while 50% of the 88 miRNA family candidates yielded levels below 10 RPM (Table 2).This low level of expression was confirmed by stemloop qRT-PCR analysis.(Figure 2C).As reported above for conserved miRNAs, the RPM values of some soybean-specific miRNAs corresponded to their relative abundance determined by miRNA qRT-PCR (gma-miR2 and gma-miR20, etc.), but several exhibited divergence (e.g., gma-miRC4, gma-miRC14 and gma-miRC32 (Figure 2C and Figure S2)).

Identification of the targets of miRNAs by degradome analysis
To identify the targets of the conserved and soybean-specific miRNAs reported here, we performed degradome sequencing to generate a total of 12.8 million short reads representing the 59 ends of uncapped, poly-adenylated RNAs.About 77.66% of the unique reads were perfectly aligned to the soybean genome (Soybean Genome V9.0, http://www.phytozome.net/search.php).These reads were subsequently screened and analyzed using the Cleaveland 3.0 software [38].A total of 42 targets in five categories (0 to 4) were identified (Table 3 and Figure S4), with 42 targets for 76 conserved and soybean-specific miRNAs belonging to 21 families (Table 3 and Figure S4).
Among these targets for the conserved miRNA families, eight fell into category 0, which represented the most abundant degradome tags corresponding to the cleavage site and matching cognate transcripts, and one of them into category 2, whose cleavage abundance was higher than the median but below the maximum.The number of identified gene targets varied among the miRNAs, from one to four (Table 3).However, miRNAs that targeted members of a gene family usually had more targets.For example, miR156 could target four members of the squamosa promoter-binding-like protein family (Table 3).Although most of the genes (36 of 42) identified were the conserved targets of these miRNAs across a wide range of plant species, some (6 of 42) had not previously been reported in other species.For example, miR169, which is known to target NF-YA (nuclear factor-Y subunit alpha) in other species, was found to target the genes encoding flavonol synthase and sulfate transporter.Similarly, miR393, which exclusively targeted mRNAs for the F-box auxin receptors TIR1 (Transport Inhibitor Response Protein 1), and several members of auxin signaling F-box protein, the growth regulating factor (AFB) gene family in plants also targeted the ribosomal protein L20 gene (Table 3).It was noted that a few identified soybean-specific gene targets fell into category 4, a lowconfidence group, and so should be further validated experimentally.Therefore, the targets falling into category 4 were not listed in the results.
Gene targets were also identified for six soybean-specific miRNAs.Of the 13 gene targets identified, 4 belonged to category 0 and 2 to category 1, while the remainder was classified into category 3 (Table 3).The soybean-specific miRNAs, like the conserved miRNAs, targeted genes of diverse functions.For example, gma-miRC23 targeted the gene encoding the auxin response factor, while gma-miRC33 targeted the gene encoding auxin-signaling F-box.Gma-miRC26 and gma-miRC61 each targeted members of gene families that encode the transcription factor TCP and the GRAS family transcription factor, respectively.Furthermore, gma-miRC15 targeted up to four members of the SBP-domain-containing-protein gene families.Hence, these soybean-specific miRNAs may be involved in the regulation of an array of metabolic and biological processes and signaling pathways.
miRNAs triggered secondary siRNAs biogenesis pathway in soybean root TAS transcripts are directed by miRNAs to produce tasiRNAs, which then guide the cleavage of other transcripts.To date, four TAS gene families have been characterized in Arabidopsis, of which the miR390-TAS3 and miR828-TAS4 pathways are conserved in plants [39,40].Here we identified TAS3 soybean orthologous genes (Glyma09g03731.1 and Glyma15g14675.1), together with their corresponding trigger miRNAs-miR390.These two genes also contained two complementary sites for gma-miR390, and the signatures were detected only at the 39 target site (Figure S5).We also found similar siRNA biogenesis patterns in the cleaved TAS3 (Table S3).Together, these data indicate that miR390-TAS3 biogenesis pathways and functions are at least partially conserved in soybean root.Because auxin signaling and modulation are essential for diverse biological processes in soybean, especially root development and seed ripening [41,42], miR390-TAS3 biogenesis-derived tasiARFs in roots could orchestrate auxin signaling that might be directly relevant to seed growth and development.In addition, four gma-miR393 target transcripts and three gma-miR1510 target transcripts in both ZP03-5373 and ZP03-5413 were identified as producing secondary siRNAs (Figure S5).Gma-miR393 and gma-miR1510-triggered secondary siRNA biogenesis pathways have been reported in soybean [15].
The targets of these phasiRNAs were identified by analysis of the soybean degradome (Table S3).Besides the ARF4, a further five novel targets of miR390-TAS3 were found.Moreover, we identified six novel targets for the five phasiRNAs derived from gma-miR393 targets, eight novel targets for the five phasiRNAs derived from gma-miR1507 targets, 29 novel targets for the 22 phasiRNAs derived from gma-miR1510 targets, eight novel targets for the seven phasiRNAs derived from gma-miR1515 targets, five novel targets for the four phasiRNAs derived from gma-miR171 targets and 15 novel targets for the nine phasiRNAs derived from gma-miR2118 targets (Table S3).

Verification of miRNA-guided cleavage of target mRNAs in soybean
To verify the miRNA-guided target cleavage, RLM-59 RACE experiments were performed to detect cleavage products of the four predicted gma-miRNAs.As shown in Figure 4, all four gma-miRNA guided target cleavages occurred at nucleotide 10 or 11 (Figure 4).Thus, all four predicted targets had specific cleavage sites corresponding to the miRNA complementary sequences.

SCN-infection-associated miRNAs
The sequencing frequencies for miRNAs in the two libraries were used as an index for estimation of the relative abundance.The expression levels in SCN-resistant soybean root and SCNsensitive soybean root were compared based on the ''reads per million'' genome-matched reads (RPM) of miRNAs.Using ZP03-5373 (RPM)/ZP03-5413 (RPM) values .5 or ,5, a total of 34 miRNAs belonging to 27 families were identified to be significantly differentially expressed.The results are shown in Table S4.Most of the differentially expressed miRNAs were up regulated in the roots of the SCN-resistant line ZP03-5373 (Table S4).Although the absolute expression level of miRNA is useful, the identification of differential expression profiles at the wholegenome level in response to endogenous cues or stresses is often desirable to detect miRNA function in particular cell processes.In order to examine if the miRNAs might play a role in SCN resistance, the expression pattern of 14 miRNAs that were expressed specifically in ZP03-5373 were analysed using qRT-PCR in SCN-infected and uninfected ZP03-5373 plants.Seven miRNAs were up regulated significantly after the SCN-infection (Figure 5), and therefore appeared to be important in SCN infection and re-generation.A search of the SCN genome sequences identified 44 potential target genes in SCN by these seven SCN-inducible soybean miRNAs, suggesting a possible function of these miRNAs in regulating the expression of these SCN genes (Table 4).
Forty-four transcripts were predicted to be potential targets of differentially expressed miRNAs.A large number of the identified targets were function proteins (Table 4), including NADH dehydrogenase, SSR2, FLP and CBN-ATP-4; i.e., the relative expression level of gma-miRC6 in ZP03-5373 was markedly higher than that of ZP03-5413 (Table 2 and figure 5).Nineteen targets of gma-miRC6 were identified.One of the gma-miRC6 targets, the SSR2 gene, encodes a translocon-associated protein subunit beta, which is associated with protein translocation across  the endoplasmic reticulum (ER) membrane.Another miRC6targeted FLP-16 protein was potentially involved in the neuropeptide signaling pathway and negatively regulated striated muscle contraction.The acyl carrier protein (ACP) is an important component in the fatty acid synthase complex.MAGOH mutations in the mago nashi (grandchildless) gene produce progeny with defects in germplasm assembly and germline development [43,44].PTCHD3 plays a role in reproduction development [45], while others are related to protein phosphorylation (AGC family protein kinase) and ATPase activity (Protein VHA-14).All of these genes are important in the development and regeneration of SCN.Ten genes were potential targets of gma-miRC18, which is involved in both the folding and transportation of proteins, and degradation pathways (Table 4).The gma-miRC32 targets encode a hypothetical NADH dehydrogenase, which is the first enzyme required in the respiratory chain pathway.

Discussion
Soybean is an important economic crop.Recently, highthroughput sequencing of sRNAs and RNA degradome has been successfully used to reveal large numbers of soybean miRNAs and their targets.A number of miRNAs have been reported to be involved in organ development [46], nutrient signaling [47], biotic and abiotic stress [48,49].These studies imply the important roles of miRNAs in soybean development and interaction with environment, which provide clues for deciphering the functions for microRNA/target modules in soybean.SCN is a significant plant pathogen responsible for an estimated $2 billion annually in yield losses worldwide.The planting of resistant soybean cultivars is the key to managing SCN population levels in the field.Despite some resistant cultivars having been developed and used, there remains a lack of understanding of the molecular basis of the resistance to this pathogen because only two major loci, rhg1 and rhg4 have been cloned [50]; the remaining quantitative trait loci (QTLs) are distributed on the other 16 linkage groups (LG) (except LG D1b and F) (soybase.net)and remain unknown.Progress in understanding the effectiveness and durability of natural plant resistance and enabling the design of novel strategies for resistance through biotechnological approaches has, therefore, been limited.Comparison of the gene expression profiles of soybean-SCN interactions has revealed distinct differences in gene expression between the resistant and susceptible reactions.Therefore, it is important to select suitable soybean lines to detect differently expressed genes.
In this study, to develop a better understanding of the molecular events associated with resistance to SCN race 4, we employed the sister lines, ZP03-5373 and ZP03-5413, in a comparative analysis of sRNA expression using deep sequencing.ZP03-5373 and ZP03-5413 have similar agronomic traits except for resistance to SCN race 4. Our previous study showed that ZP03-5373 was resistant but ZP03-5413 was susceptible to SCN 4, suggesting that some differentially expressed genes may have negative impacts on syncytium development and maintenance.
SCNs are highly evolved sedentary plant endoparasites that can enter soybean roots to successfully parasitize plants.RNA interference (RNAi) involving host-induced gene silencing in parasites has been reported [5].A potential mechanism underlying the involvement of miRNAs in controlling cyst nematodes is proposed.Here, the candidate targets of differentially expressed miRNAs in SCN were predicted (Table 4).Our results predicted the existence of a novel miRNA-mediated regulatory cascade involved in the SCN life cycle in soybean root.These observations demonstrate the relevance of the targeted genes of SCN during the nematode life cycle and, potentially more importantly, suggest that an effective resistance to cyst nematodes in soybean may be achieved using this technology.But which should be confirmed by experiment in the future.

Conclusions
This study describes large scale cloning and characterization of two genetically related soybean sister lines miRNAs, phasiRNAs and their potential targets, we also found that the expression of 34 miRNAs differed significantly between the two lines.Seven ZP03-5373-specific miRNAs were differentially expressed after SCN infection.Forty-four transcripts from SCN were predicted to be potential targets of ZP03-5373-specific differential miRNAs.These findings suggest that miRNAs play an important role in the soybean response to SCN and providing the foundation for further characterization of their roles in the regulation of diverse physiological processes.

Plant materials
Two genetically related soybean lines, Zhongpin03-5373 (ZP03-5373) and Zhongpin03-5413 (ZP03-5413), which are resistant and RNA extraction and preparation of sRNA and degradome cDNA libraries for Solexa sequencing Soybean root total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA).Total RNA was size-fractionated by 15% denaturing polyacrylamide gel electrophoresis, after which small RNA fragments of 18-28 nt were isolated from the gel and purified.The small RNA molecules were then sequentially ligated to a 59 adaptor and a 39 adaptor and converted to cDNA by RT-PCR following the Illumina protocol.The concentration of the sample was adjusted to ,10 nM and a total of 10 mL were used in a sequencing reaction.The purified cDNA library was sequenced on an Illumina GAIIx.
The degradome library was constructed as described previously [33].For the short RNA libraries, the degradome cDNA library was sequenced on an Illumina GAIIx.

Bioinformatic analyses
After masking adaptor sequences and removal of contaminated reads, the clean reads were filtered for miRNA prediction.First, reads that matched rRNA, tRNA, snRNA, snoRNA, repeat sequences, and other ncRNAs deposited in Rfam (http://www.sanger.ac.uk/software/Rfam) [51] and the GenBank noncoding RNA database (http://www.ncbi.nlm.nih.gov/) were discarded.The retained 18-28-nt reads were mapped onto the genome of soybean, using V 9.0 (http://www.phytozome.net) by the bowtie2 software.All perfectly matched sRNAs were retained for miRNA prediction.After rigorous screening, all retained sequences of 18-28 nt with a frequency of three or more copies were considered potential miRNAs.We then attempted to align the predicted miRNAs to all soybean known mature miRNA sequences in miRBase, version 19.0 [51] to identify novelty.Finally, secondary structure prediction of individual miRNAs was performed with the MFOLD software (Version 2.38, http://mfold.rna.albany.edu/?q=mfold/RNA-Folding-Form) using the default folding conditions [52], and novel miRNAs were predicted using the psRobot software [53].The identification of phased transcripts in soybean was performed by a method described previously [54].
The degradome analysis and the classification of target categories were performed using CleaveLand 3.0 [38].Small RNA target prediction was run against the transcriptome of interest.The alignment scores (using the [8] rubric) for each hit up to a user-defined cutoff were calculated, full RNA-RNA alignments were printed, and the 'cleavage site' associated with each prediction was also calculated.The cleavage site is simply the 10 th nt of complementarity to the aligned sRNA.For randomized queries, no alignments were retained; however, concise records of each predicted target for the random queries were retained, including the predicted cleavage sites.We also used the psRobot software to identity the targets of phasiRNAs.miRNA targets were predicted in SCN using the microTar software (http://tiger.dbs.nus.edu.sg/microtar/)[55].

Figure 1 .
Figure 1.Length distribution of redundant and unique sRNA sequences.The length distribution of redundant and unique sRNAs in ZP03-5373 (a) and ZP03-5413 (b).The 21-nt of redundant is the predominant sRNA species and the 24-nt of unique is the most abundant.doi:10.1371/journal.pone.0110051.g001

Figure 3 .
Figure 3. Five novel phasi-acting siRNA biogenesis pathways in soybean root.The abundance of each secondary siRNAs is plotted (left).The phasing secondary siRNAs corresponding to the miRNA cleavage sites are highlighted in red.The miRNA complementary sites are shown with red arrows.The length distribution is plotted on the right (middle).The phasing radial graph is represented next to this (right).Each spoke of the radial graph represents 1 of the 21 phasing registers, with the total number of sRNAs mapping to that register plotted as distance from the center.A, sense transcript; AS, antisense transcript.doi:10.1371/journal.pone.0110051.g003

Figure
Figure S1 Secondary structures of 71 putative lessconserved soybean miRNAs and miRNAs.Pink section represents miRNA-5p; yellow section represents miRNA-3p.(PDF) Figure S2 qRT-PCR results.qRT-PCR confirming express pattern of miRNAs in ZP03-5373 and ZP03-5413.The expression levels were normalized against the U6 RNA.(JPG) Figure S3 Secondary structures of 75 putative soybeanspecific miRNAs and miRNAs counterparts.Pink section represents miRNA-5p; yellow section represents miRNA-3p.(DOCX) Figure S4 degradome T-plot.We used reads in plotting the cleavages on target mRNAs, which were referred to as 'target plots' (t-plots) by German et al [17].Signature abundance throughout the length of the indicated transcripts is shown.miRNA:mRNA alignments along with the detected cleavage frequencies are shown.The frequencies of degradome tags with 59ends at the indicated positions are shown in black, with the frequency at position 10 of the inset miRNA target alignment highlighted in red.(PDF) Figure S5 The small RNAs corresponding to the miRNA targets.The abundance of each secondary siRNAs is plotted (A).The phasing secondary siRNAs corresponding to the miRNA cleavage sites are highlighted in red.The miRNA complementary sites are shown with red arrows.The length distribution is plotted on the right (B).The phasing radial graph is represented next to

Table 1 .
New members of conserved and less-conserved soybean miRNAs.

Table 3 .
Identification of soybean miRNAs targets using the degradome.

Table 4 .
Potential targets in SCN for miRNAs expressed at a high level in the ZP03-5373 line.

Table 4 .
Cont.The target gene is the transcript identified from the SCN ESTs.(http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=51029); b E-value was calculated according to Blast and should be less than 1.00E-5.doi:10.1371/journal.pone.0110051.t004susceptible, respectively, to SCN race 4 were used in this study.The two sister lines, ZP03-5373 and ZP03-5413 were developed from the cross of two SCN resistant parents ''Jin 1265'' ?''Hartwig''.The former was resistant and the latter was susceptible to SCN race 4. Elite line Jin 1265 was derived from cultivar Hupizhi Heidou for its resistance.Thus, ZP03-5373 and ZP03-5413 have the same genetically pedigrees but different resistance to SCN race 4, which provided an opportunity to gain further insight into the underlying genetic control of resistance.Soybean were grown in a glasshouse at 22-25uC with a 16 h light/8 h dark photoperiod and light intensity of .8000lx.Roots from 3-weeksold seedlings were collected and used for RNA extraction.And was used for small RNA expression and degradome analysis. a Table S2 Known miRNAs identified in G. max.Table S4 Differential expressed soybean miRNAs.(XLS) Table S5 miRNA and primer sequences.(XLSX)