Identification of Soybean MicroRNAs Involved in Soybean Cyst Nematode Infection by Deep Sequencing

Soybean cyst nematode (SCN), Heterodera glycines, is the most devastating pathogen of soybean worldwide. MicroRNAs (miRNAs) are a class of small, non-coding RNAs that are known to play important role in plant stress response. However, there are few reports profiling the miRNA expression patterns during pathogen stress. We sequenced four small RNA libraries from two soybean cultivar (Hairbin xiaoheidou, SCN race 3 resistant, Liaodou 10, SCN race 3 susceptible) that grown under un-inoculated and SCN-inoculated soil. Small RNAs were mapped to soybean genome sequence, 364 known soybean miRNA genes were identified in total. In addition, 21 potential miRNA candidates were identified. Comparative analysis of miRNA profiling indicated 101 miRNAs belong to 40 families were SCN-responsive. We also found 20 miRNAs with different express pattern even between two cultivars of the same species. These findings suggest that miRNA paly important role in soybean response to SCN and have important implications for further identification of miRNAs under pathogen stress.


Introduction
The soybean cyst nematode (SCN), Heterodera glycines Ichinohe, is an obligate sedentary endoparasite that causes extensive damage to soybean, Glycine max (L.) Merr., worldwide and accounts for over one billion dollars loss annually in the US [1]. Current nematode control strategies include nematicides, crop rotation and resistant cultivars, but each has serious limitations. The lack of understanding of nematode virulence has hampered the ability to devise novel management tactics [2]. Advances in genomics during the last years provide opportunities for us to investigate molecular mechanism of the G. max-H. glycines system. Nematode genome sequencing projects and comparative analysis contribute to insight the nature of the evolution of plant parasitism [3,4,5]. Laser capture microdissection (LCM), a refined approach to obtain homogeneous cell samples, are enabling the discovery of the genes and pathways contributing to the formation of unique, highly specialized feeding cells induced within plant roots by parasitic nematodes [6,7]. Emerging tools, such as microarray, RNA-seq, RNA interference hold tremendous potential to uncover the mechanism of soybean-SCN interactions [8,9,10,11,12]. Previous studies indicate that plants reprogram gene expression at the transcriptional, post-transcriptional and post-translational levels to reduce the damage to stress. Increasing evidences unraveling that host endogenous smallRNAs are essential in this gene expression reprogramming process [13].
MicroRNAs (miRNAs) are a class of 21-24 nucleotide long noncoding RNAs that negatively regulate gene expression in animals and plants [14]. MiRNAsrole in plant was first demonstrated in plant growth, development, hormone signaling, and RNA metabolism [15,16]. Recently, it was reported that miR393 was induced by a bacterial flagellin-derived PAMP, Flg22 in Arabidopsis. MiR393negatively regulates the expression level of the F-box auxin receptor genes and restricts Pseudomonas syringae growth. Combined with race-specific resistance test, the results suggest that miR393 has a role in imparting basal resistance but not race-specific resistance [17]. Researchers find 10 miRNAs are down-regualted in loblolly pine response to rust fungus [18]. Two miRNAs, bra-miR158 and bra-miR1885, were greatly upregulated when Brassica rapa was infected by Turnip mosaic virus [19]. Genome-wide small RNA profiling analysis from nodules or from roots inoculated with symbiotic rhizobia identified a number of nodulation-associated miRNAs from different legumes [20,21,22,23,24]. Mtr-miR166, Mtr-miR169 of M. truncatula and three soybean miRNAs, Gma-miR482, Gma-miR1511, Gma-miR1512 have been functionally analysed and associated to rhizobial symbiosis [25,26,27]. These increasing evidences show that miRNAs are involved in plants response to biotic stresses.
To date, many conserved or novel miRNAs have been identified in soybean by computational analysis or high-throughput sequencing [28,29,30,31,32,33]. MiRNAs related to Bradyrhizobium japonicum, Phytophthora sojae response are obtained by comparative analyzed [20,34]. However, there is no report of miRNAs associated with SCN infection. Here, using solexa analyzer, we identified a diverse set of small RNAs which response to SCN infection in both resistance and susceptibility soybean roots. Totally, 21 new miRNAs belong to 19 families were obtained. Moreover, we also found that 40 soybean miRNA families showed differential expression level in response to soybean cyst nematode infection. Some potential targets of these miRNAs were stress responsive genes.

Plant material and RNA isolation
Glycine max, cultivar ''Harbin xiaoheidou'' (HB) were chosen as a soybean cyst nematode race 3 resistant standard and ''Liaodou 10'' (L10) were used as SCN race 3 sensitive standard. The plants were grown in greenhouse, cyst inoculation samples were propagated on SCN race 3 infected-soil which was collected from the pilot field of Shenyang Agriculture university, control samples were grown on autoclaved blank soil. Root samples were collected 30 days after seedlings emerge. All the tissues were harvested and immediately frozen in liquid nitrogen and stored at 280uC. 80 mg of plant tissues were ground to powder with liquid Nitrogen. Total RNA was isolated from soybean root samples with TRIzol reagent (Invitrogen, USA) according to the manufacturer's protocol. Total RNA was treated with DNaseI (NEB, USA) for 30 min at 37uC, and then purified by ethanol precipitation. The RNA samples were qualified by Agilent 2100.

Small RNA libraries construction and SBS sequencing
Small RNAs isolation and library construction were performed as described by Hafner Markus, et al [35]. Firstly, total RNA was size fractioned on 15% TBE urea polyacrylamide gel and 18-30 nt fraction was collected. The 59 RNA adapter was ligated to the RNA pool with T4 RNA ligase. Ligation products were gelpurified. Then ligated with 39RNA adapter and purified. Small RNAs with both side adaptors were subjected to reverse transcription, 15 cycles of PCR reaction was performed to produce sequencing libraries, and the products were gel-purified.
Small RNAs samples were sequenced using Illumina Genome Analyze II at Beijing genomics Institute (Shenzhen). 10 pM of each sample was used for cluster generation. After hybridization of sequencing primer, 35 cycles of base incorporation were carried out on the 1G analyzer. More than 16 million 35 nt reads were generated for each sample.

Data analysis and new miRNAs identification
Image analysis and base calling were performed using Illumina Pipeline. Low quality reads (tags with unknown nucleotide 'N'), empty tags (sequence with only adaptor sequence), tags less than 17 nt, and low complexity tags were trimmed with our own perl script. Removed 39 and 59 adapter, purified tags were mapped to the soybean genome sequence (http://www.phytozome.net/cgibin/gbrowse/soybean) using SOAP v1.11 (BGI). Tags that perfectly matched soybean genome were obtained for subsequent analysis. These tags were used to search for Rfam database (http://rfam.sanger.ac.uk/), unique reads mapped to tRNA (transfer RNA), rRNA (ribosomal RNA), snRNA(small nuclear RNA), snoRNAs (small nucleolar RNAs), known miRNAs (microRNAs) were removed. Repeat overlapping sequences were annotated as repeat-associated small RNAs. Tags matched to gene exon region were excluded from further analysis.
Uni-tags resulted from the pre-processing steps, with copy number more than one and less than 24 nt were tested by Mireap developed by BGI (Shen Zhen) to predict novel miRNAs. The algorithm includes two dependent parts. First of all, candidate miRNA sites are screened out from breakpoints defined by small RNAs mapping. At each genomic locus, two sequences covering the read were extracted for secondary structure analysis, one sequence extending 160 nt upstream and 20 nt down-stream from the read, and the other covering 20 nt upstream and 160 nt downstream of the read. Then, a minimal stringent criterion was used to select miRNA candidates. The secondary structure must have a hairpin with at least 18 paired nucleotides in its stem region. The free energy of the hairpin should be less than -18 kCal/mol. MiRNA and miRNA* had to reside in different arms of a hairpin structure, each with no more than 6 unpaired bases. The maximum bulge over the miRNA/miRNA* duplex was not more than 4 bases, and asymmetry of the miRNA/miRNA* duplex equal to or less than 4. In addition, sequencing of both miRNA and miRNA* required that the miRNA/miRNA* duplex had 39 overhangs at both ends, a typical feature of Drosha and Dicer processing [36]. The resulting set of sequences and their respective RNA structures were further analyzed to distinguish genuine miRNA precursors from other RNAs containing similar RNA structures.
Mature miRNA sequences were used as queries to search for potential target mRNAs in the Glycine max database (http:// compbio.dfci.harvard.edu/index.html, DFCI gene index release 16) using the web-based computer psRNATarget Server (http:// biocomp5.noble.org/psRNATarget/). Results from these analyses were inspected on the Phytozome to obtain the loci and protein annotation. The rules used for miRNAs target gene prediction were based on suggestions by Allen et al. and Schwab et al. respectively [37,38]. No more than 4 mismatches between the sRNA and target were allowed. The minimum free energy (MFE) of the miRNA/target duplex should be . = 74% of the MFE which the miRNA bound to its perfect complement [39,40].

Differential expression analysis of miRNAs
The frequency of miRNA was normalized to TPM (number of transcripts per million clean tags) in order to compensate for variable numbers of tags generated for each sample. The foldchange between SCN-infected and control sample were calculated. To avoid divided by 0, we use 1 instead of tag that was not detected in sample. Then the statistically analysis was performed according to Poisson distribution.

Quantitative PCR analysis
The profiles of 16 different expressed miRNAs were assayed by Real-time PCR. About 200ng qualified total RNA for each sample were used for this analysis.
The miRNA gene specific reverse transcript stem-loop primers were designed using Primer 5.0 following the method described by Varkonyi-Gasic et al [41]. The reverse transcript reaction was performed with Superscript II (Invitrogen, USA) following the the manufacturer's protocol. Reverse transcript products were used as template for Real-time PCR and all reaction were assayed in triplicates. The reactions were performed in Stepone RT-PCR machine (Applied Biosystems, USA) using SYBR Premix ExTaq kit (Takara, China). PCR cycling began with template denaturation and hot start Taq activation at 95uC for 1 min, then 40 cycles of 95uC for 5 sec, and 60uC for 30 sec performed and data collected during each cycle at the 60uC extension step. The U6 snRNA was selected as a reference gene for normalization. Relative quantitation of the miRNAs amplified was performed by the comparative DDC T method. The threshold cycle (C T ) indicated the fractional cycle number at which the amount of amplified target reaches a fixed threshold.

High-throughput sequencing of soybean small RNAs
To investigate the role of soybean miRNAs in response to soybean cyst nematode (SCN) infection, two soybean cultivars, Liaodou10 (L10), a soybean cyst nematode race 3 susceptible cultivar and Harbin xiaoheidou (HB) which resistant to SCN race 3 were selected for this study. Totally, four small RNA (sRNA) libraries were constructed from the two cultivars. The SCNinfected sRNA libraries from HB and L10 roots were named SHB, SL10 respectively. While sRNAs libraries from un-infected roots of HB and L10 were constructed as control, which designated as SHBC and SL10C. These four soybean root small RNA libraries were sequenced by Illumina Genome Analyze II, yielding a total of 69,861,808 sRNA raw reads. After removing the low quality reads, 11,831,532, 10,980,930, 10,494,778 and 11,238,552 clean reads were obtained for SHB, SHBC, SL10 and SL10C, respectively. These sRNAs were consisted by 11,317,750 unique sequences as shown in table 1, and these reads were searched against the public soybean genomic and expressed sequence tags using SOAP program, leading to 5,103,481 genome-matched reads [43].
Majority of sRNAs were sequenced only a few times although some sRNAs were highly abundant and present thousands of times in one sample. About 36% of the sRNAs were sequenced only once suggests soybean roots have large and diverse sRNAs population.
More than 90% of clean reads in this study were longer than 18 nt. Size distributions based both on total abundances and unique reads were assessed and shown in figure 1. For total abundance, more than 75% sRNAs were 19-24 nt in length, with 21, 22, 24 nt being the major length classes. For the proportion of unique reads, the 24 nt reads were prevailing in all four libraries and the 21, 22 nt reads were less abundant. Analysis revealed that there is a variation in abundance of sRNAs, the 21, 22 nt sRNAs showed higher redundancies than 24 nt sRNAs. The observation indicated the differently expression pattern of distinct categories of sRNAs in the same sample. This result was consistent with those of Citrus rifoliate, Arabidopsis, Medicago truncatula [44,45,46].
The common/specific sRNAs in SHB vs SHBC and SL10 vs SL10C were analyzed for total and unique tags, the results highlight the fact that sRNAs common in different libraries had more abundance than those library specific sRNAs (Table S1).

Identification of known miRNAs
In this study, miRBase release 18.0, which contains 391 soybean mature miRNAs, was searched for soybean known miRNAs [47]. A total of 364 known soybean mature miRNAs were identified in at least one library, among which 283 were sequenced in all 4 four libraries and 53 miRNAs were detected only in one or two libraries, which was summarized in table 2. There were 27 known soybean miRNAs, containing 7 conserved miRNAs and 20 fabaceae or soybean specific miRNAs, were absent in our data set (Table S2). This is probably due to low expression levels or stage-/tissue-specific RNA transcription. Two of reported miRNA*, gma-miR1507c* and gma-miR1516* were also identified.
Gma-miR1507ab was the most abundant miRNA in all the 4 samples in our study, followed by members of MIR156. The top 30 abundant miRNAs which occupied more than 85% of expressed miRNA tags were generally unanimous between samples. Most of the top 30 abundant miRNAs were conserved miRNAs such as, members of MIR156, MIR166 and MIR168, while three legume specific miRNAs, miR1507, miR1509 and miR1510 were highly abundant in all four libraries. Further analysis revealed that conserved miRNAs had relative high abundant in general, this observation suggested that conserved miRNAs may be essential for controlling basic cellular and developmental in plants. More than half of the miRNAs showed very low levels of expression, with fewer than 100 reads.
Within the tags that mapped to known miRNAs over 55% have length 21 nt, followed by 20 nt and 22 nt reads. Among the miRNA reads, 24 nt length was less than 1%. This result was consistent with current understanding. The canonical miRNAs are 21 nt length, while canonical siRNAs are 24 nt length [48]. About 94% of these miRNA reads start with a uracil contrast to 30% of all small RNAs. These distribution were examined in all four samples individually, there were no obvious difference among them.
We observed there are diverse isoforms of mature miRNAs from the same precursor and the 39 ends of miRNA showed stronger heterogeneity than the 59 ends. As exemplified by miR1509a in Figure S1, the highlighted 39 ends isoform were highly accumulated. Some most abundance tags mapped to hairpin were not previous annotated miRNAs. Gma-miR156b is an extreme example, the former annotated miRNA were not detected in all 4 samples, indicating the tagt0006055 should be regarded as the final functional molecular.

Prediction and validation of novel soybean miRNAs
Small RNAs perfectly matched to soybean genome sequences were searched against Rfam, ncRNA (non-coding RNA) of NCBI, EST database to eliminate other ncRNAs, such as rRNA, tRNA, snRNA, snoRNA and degradation products of protein-coding transcripts. Small RNAs were aligned to repeat associated RNA to find matched tags in the data set. With the analysis, 3,571,130 tags named unann were obtained for novel miRNA detection (Table S3). Characteristic hairpin structures of those reads were explored using Mireap. Then precursors were manually checked  Novel miRNA candidates were further assigned to miRNA families using sequence similarity to other known miRNA in miRBase database. Among the 21 novel miRNAs, 18 miRNAs have not been found in eukaryotes, except soy_1, soy_2, soy_3. All the three miRNAs were members of MIR166, which conserved in both monocots and dicots. The 18 miRNAs without homolog in database, along with the known soybean miRNAs were used for multiple alignments using ClustalW and the miRNA families were assigned based on the dendogram tree [49]. None of the miRNAs was assigned to families based on sequence similarity, while the miRNAs will get families with further study.
To validate the predicted new miRNAs, stem-loop RT-PCR were performed to examine whether the miRNAs were expressed in the soybean roots. Eight miRNAs which had more abundance were selected to confirm the analysis. Primers used in this experiment are showed in S5. The PCR products were about 60bp in length, and all the 8 miRNAs were found to be expressed in the soybean roots ( Figure 2).

SCN-infection associated miRNAs
Although the absolute expression level of miRNA is useful, identification of differential expression profile at the whole genome level in response to endogenous cues or stress is often desirable to detect miRNA function in particular cell processes. The sequencing frequencies for miRNAs in the four libraries were used as an index for estimating the relative abundance, expression level between SCN-infected soybean root and un-infected root were compared based on the ''transcripts per million''(TPM) of miRNAs ( Figure 3).
We used fold-change values $ 2 or #0.5 combined with Pvalue ,0.01 as the threshold, a total of 101 miRNAs belong to 40 families were identified to be significantly differentially expressed in at least one cultivar, the results were showed in Figure 4. Most of the differentially expressed miRNAs were down regulated during the SCN-infection and only 6 miRNAs were up-regulated. Thus, down-regulation of miRNAs appeared to be important in SCN infection. MiR171c and miR319 were highly induced by  MiR159b, miR398 and miR2119 showed the highest alternation with 12, 8.5 and 8 fold changes in HB while miR862b increased to 25 fold after SCN infection was the highest changed miRNA in L10. There were 16 miRNA families which were co-regulated in both two genotypes, among them, miR156, miR162, miR166a, miR167, miR319, miR397, miR398, miR408 were conserved miRNAs between plants species, two miRNA families, miR2119 and miR3522 were conserved in fabaceae, while miR1520, miR4365, miR4387, miR4413, miR4996 and miR5671 were soybean specific miRNAs. The results indicated that both conserved miRNAs and soybean specific miRNAs participated in nematode defense. HB, the SCN race 3 resistant cultivar, specific response miRNAs contained 15 members, while there were 9 susceptible cultivar L10 specific response miRNAs. Notably, some members of one miRNA family had varied frequencies and the members showed differential expressed pattern during SCN-infection. MiR1507 had three members, gma-miR1507ab frequency was more than 40 times of gma-miR1507c. MiR1507ab was down-regulated after SCN-infection while miR1507c was induced by nematode infection. We also found that members of MIR390 and MIR171 had different expression pattern showed in Table S4. Differential expressed pattern indicating the function divergence in miRNA families.
To confirm the Solexa sequencing results, we performed qRT-PCR to further validate several miRNAs that were detected as differential expressed. We found that most of the qRT-PCR results were consistent with the deep sequencing data (figure 5). However, there were differences between the sequencing results and qRT-pCR data, the log 2 ratio of gma-miR1512 in L10 was 20.1 with sequencing data, but the qRT-PCR ratio was 0.36, this could been due to sequence bias induced by the small RNA library construction or to the different normalization methods in this two strategy.

Targets prediction of soybean miRNAs
MiRNAs regulate expression of specific gene via hybridization to mRNA transcripts to promote RNA degradation, inhibit translation or both [50]. We predicted the potential targets of miRNAs using the webserver, psRNATarget [51]. For the  miRBase registered soybean miRNAs, a total of 293 target genes were identified showed in Table S5. Combined with previous studies, we investigated that most target of conserved miRNAs were transcript factors which involved in plant growth and stress responses. Some of miRNA families, such as MIR156, MIR160, MIR164, MIR393, MIR396, MIR1510 were conserved in their targets. The targets prediction results of MIR1507, MIR390 and MIR171 were checked manually, gma-1507a and gma-1507b had the same target genes, while potential targets of gma-1507c were totally different. This was consistent with the expression pattern between miRNAs. Members of MIR390 and MIR171 come to the same results.
Among the 21 novel identified miRNAs, 16 miRNAs belonged to 14 families had predicted targets, the results were showed in table 4. Soy_5, soy_8, soy_13, soy_17, soy_20 target only one locus, while other miRNA candidates had multiple distinct targets genes. Soy_1, soy_2 and soy_3 were members of MIR166, targets of the three miRNA candidates were HD-ZIP transcript factor which were conserved with other miRNA members. Function of new miRNA tragets were diverse and were not enriched in transcript factors compared to conserved miRNAs. Calmodulin, serves as an intracellular Ca 2+ receptor and mediates the Ca 2+ regulation of cyclic nucleotide and glycogen metabolism, secretion, motility and Ca 2+ transport, also a dynamic component of the mitotic apparatus. Calmodulin was detected as soy_9 potential targets indicating the miRNA may participate in cellular regulation [52].

Discussion
High-throughput sequencing is an efficient strategy to study miRNAs at the whole genome level and has dramatically expanded the number of miRNA families known to exit in plants. As to soybean, 391 miRNAs have been well annotated by deep sequencing and registered in miRBase database [47]. In present study, we sequenced small RNA libraries of two genotype soybean root treated by soybean cyst nematode, totally got 69,788,759 reads and 11,317,750 unique reads. The large amount of small RNA data revealed a diverse and complex small RNA population. Moreover, the large data set allows us to identify miRNAs with low abundance. In total, 364 known miRNAs and 21 novel miRNAs were detected, 171 of them had frequency less than 20. Combined with the previous report, 40% the soybean miRNAs are belonged to soybean or lineage specific, while roughly 70% known Arabidopsis miRNA families lack apparent homologs in outside of the brassicasseae [53]. This result indicates that miRNAs discovered in soybean has not reached saturation. Novel miRNA detected in this study showed relatively low expression levels, which was consistent with previous studies [54,55,56]. Thus much larger dataset are needed to identify new miRNAs than previous studies.
In our study, gma-miR1507ab was the most abundant miRNAs in all four samples followed by miR1509, miR1510 and miR2118. Although miR1509a and miR2019 were annotated other than 22 nt, based our data, previous register should be corrected. In M. truncatula, peanuts and common beans, miR1507, miR1509, and miR2118 were highly abundant 22-nt miRNAs [57]. The conservation and abundant expression of these 22 nt in legumes indicated unique evolution ways of these lineage-specific miRNAs. More interestingly, we found that all the four miRNA families were regulators of legume NB-LRR-coding genes. Furthermore, potential targets of miR5039, miR5041, miR5374, miR5376, miR5668 were plant NB-LRR transcripts. Until now, the only function of plant NB-LRR proteins is recognize pathogen effectors and activate defense responses [58]. The results implicating that the specialized group of miRNAs in legume target NB-LRR genes are crucial regulator to some legume specific plant defense to pathogens.
MiR319 was shown to be drought-responsive in Arabidopsis, rice, and sugarcane, while we firstly report miR319 was involved in biotic stress [61,62,63]. Gma-miR319 was the single up-regulated conserved miRNA family in both two cultivars. In contrast to deeply conserved miRNA families, most young miRNAs were weakly expressed, more divergent, and tend to lack targets. Thus, function information about non-conserved miRNA is rare. The 20 legume-or soybean specific miRNAs related to SCN infection were identified in current study, suggesting a number of nonconserved miRNAs function as regulators of genes in plant genomes, and play an important role in plant pathogen defense. The complex plant response to SCN infection spans the early migratory stages during penetration and migration through roots into the later sedentary stages of syncytium induction and feeding which have inspired numberous studies to identify genes involved in SCN. Here we report the miRNAs related to SCN infection. In previous study, 43 SCN-responsive miRNAs were identified. More than 80% of miRNAs were suppressed by SCN infection in both cultivars. Similarity, 10 of 11 miRNAs were down-regulated in the stem galls induced by the rust fungus Cronartium quercuum f. sp. Fusiforme [18]. While in Arabidopsis-H. schachtii interaction study, all the 14 differentially expressed miRNAs were statistically significantly down-regulated at the 4 dpi (day post-inoculation) while at 7 dpi, 5 miRNAs were up-regulated and 7 were down-regulated [59]. Indicating the tested miRNAs have complex expression profiles at different stage of cyst nematode infection which is consistent with gene expression pattern during SCN-infection [9,10]. Targets of SCN-responsive miRNAs include transcription factors, general stress-responsive genes, disease resistance responsive proteins, PR proteins. All these genes expression pattern were checked in previous studies, about 50% miRNAs showed negative correlations between miRNA accumulation and target gene mRNA abundance [60]. These discrepancies suggest that other regulatory mechanisms could be involved in regulating these genes.
MiR319 was the single up-regulated conserved miRNA family in both two cultivars, which was shown to be drought-responsive in Arabidopsis, rice, and sugarcane [61,62,63]. In tomato, miR319 was involved in Solanum lycopersicum, and showed opposite expression patterns in leaves [64]. The diverse functions of miR319 were due to the target gene TCP that presumably control cell divisions in stress response. Interestingly, during the same stress, several miRNAs had differential express pattern in the two cultivars. MiR169 was up-regulated in the sensitive cultivar, but down-regulated in the tolerant. MiR169 targeted nuclear factor Y subunit, which can regulate expression level of some stressresponsive genes [65]. MiR169 guided regulation of these transcription factors appeared to be important in variety of abiotic stress, including drought, cold, salinity, UV-B radiation, and heat stress [61,62,63,65,66]. Furthermore, MtHAP2-1, a new transcription factor of the CCAAT-binding family identified in M. truncatula, is regulated by miR169. MtHAP2-1 is essential for nodule meristematic persistence in M. truncatula, and that miR169 confers spatial and temporal accuracy of nodules development [67]. Contrarily, gma-miR390b was induced by SCN in HB and suppressed in L10. MiR390 was involved in auxin signaling and regulated root growth [68]. Gma-miR390c were down-regulated in HB but without significant changes in L10. Expression patterns of MIR390 members are different even in single cultivar which indicating function diversity of this miRNA family. With other 15 miRNAs up-or down-regulated only in one cultivar, the genotypespecific regulation of miRNAs might be part of the reason why the two soybean cultivars had differential tolerance to soybean cyst nematode.
In summary, our study have identified miRNAs form soybean roots, analyzed their expression pattern during soybean cyst nematode infection and predicted the putative targets of these miRNAs. Further function characterization of SCN-responsive miRNAs and targets will contribute substantial new knowledge to our understanding of the complex soybean-SCN pathosystem. Figure S1 Diversification of mature miRNA from precursors. Detected diverse isoforms of soybean miRNA, the most abundant sRNA are underline in red, sRNAs in red frame are annotated mature miRNA in miRBase.

(TIF)
Table S1 Statistics of common sRNAs from control and SCN-infected. Statistics of common sRNAs from control and SCN-infected sample in both total and unique reads. (XLSX)

Author Contributions
Conceived and designed the experiments: XL WD YD. Wrote the paper: XL. Carried out the small RNA cloning and library sequencing: XL SZ. Performed the quantitative RT-PCR analysis: XL XW DL. Performed the data management, processing and analysis: XL. Prepared the plant materials: XW DL. Critically revised the article: WD YD. Obtained the funding: WD YD.