BACKGROUND
small, non-coding rnas have been characterized in plants as important factors involved in gene expression regulation in developmental processes  <cit> , as well as adaption to biotic and abiotic stress conditions  <cit> .
in general, small rnas are grouped into two major classes: micrornas  and short-interfering rnas .
these two classes of small rnas cannot be discriminated by either their chemical composition or mechanism of action  <cit> .
however, sirnas and mirnas can be distinguished by their origin, evolutionary conservation and the types of genes that they silence  <cit> .
in this way, mirnas are well differentiated due to some particular characteristics.
these characteristics include the following: derived from genomic loci distinct from other recognized genes, processed from transcripts that can form local rna hairpin structures, and usually, mirnas sequences are nearly always conserved in related organisms  <cit> .
in plants, mirna genes are transcribed by rna polymerase ii enzymes  generating primary mirna .
the pri-mirna forms an imperfect fold-back structure, which is processed into a stem-loop precursor  by nuclear rnaseiii-like enzymes called dicer-like proteins   <cit> .
the resulting pre-mirna contains a mirna:mirna* intermediate duplex, formed by a self-complementary fold-back structure.
a mature mirna sequence can range from  <dig> to  <dig> nucleotides  in length and act as a regulatory molecule in post-transcriptional gene silencing by base pairing with target mrnas.
this leads to mrna cleavage or translational repression, depending on the degree of complementarity between the mirna and its target transcript  <cit> .
the same mature mirna can also present several variants of their sequence in length.
these populations of mirna variants are called isomirnas, which are isoforms of micrornas  <cit> .
they are caused by an imprecise or alternative cleavage of dicer during pre-mirna processing  <cit> .
isomirnas have been recently identified in both plants and animals  <cit> .
the first plant mirnas were described in arabidopsis thaliana  <cit>  and later in other species.
currently, mirnas have been reported in  <dig> plants species, and all of their sequences have been deposited in a publicly-available mirna database, mirbase   <cit> .
several mirnas have been identified in plants, and they are characterized in a wide variety of metabolic and biological processes in plants with important functions in development  <cit> , phytohormone signaling  <cit> , flowering and sex determination  <cit>  and responses to biotic and abiotic stresses  <cit> .
in soybean  merrill), the major legume crop worldwide, subramanian et al.
in  <dig>  <cit>  identified  <dig> novel mirna families for the first time.
in this study, the role of mirnas in soybean-rhizobial symbiosis was investigated  <cit> .
during that same year, zhang et al.
<cit>  used a comparative genome-based in silico screening of soybean est databases and quantitative pcr to provide evidence for  <dig> mirnas belonging to  <dig> families.
a second study involving mirnas and soybean root nodules was performed by wang and colleagues  <cit> .
they identified  <dig> mirnas belonging to  <dig> mirna families.
the identification of nine novel mirnas in wild soybean  was also reported by chen et al.
<cit> .
another study looked at four different soybean tissues  and identified  <dig> novel soybean mirnas  <cit> .
recently, song and coworkers  <cit>  identified  <dig> new mirnas and their related target genes from developing soybean seeds.
although these studies resulted in a large number of mirnas identified in soybean, none of them looked at micrornas with respect to biotic and abiotic stresses.
drought is the major abiotic stress factor to negatively affect soybean productivity around the world.
the impact of limited water during the flower formation can cause shorter flowering periods  <cit> , and water stress during the later phases of soybean reproductive development has been reported to accelerate senescence, which decreases the duration of the seed-filling period  <cit> .
with regards to biotic stress, asian soybean rust  is a foliar disease caused by the fungus phakopsora pachyrhizi sydow & sydow.
this pathogen presents a rapid aerial spread and a high capacity to colonize leaf tissue and, to a lesser extent, stem and pods  <cit> .
asr is one of the most severe diseases on the soybean culture, which causes damage between 10% and 90% in the different regions where it has been identified  <cit> .
this disease is the main threat in soybean-producing countries.
currently, there are  <dig> mirnas identified in glycine max ; however, none of these mirnas were associated with water deficit or asr stress conditions.
we consider that the identification of these mirnas is important to understanding small rna-mediated gene regulation in soybean roots under water deficit stress and in leaves during rust infection.
in this context, our goal was to identify new mirnas and to discover those that may be regulated by water deficit and soybean rust stress.
using high-throughput sequencing, we constructed four libraries of small rnas from the roots of drought-sensitive and tolerant seedlings in response to control or water deficit conditions.
we also constructed four libraries from leaves of rust-susceptible and resistant seedlings with mock and infected conditions.
a set of eight small rnas libraries was analyzed from soybean plants.
a total of  <dig> mirnas were detected in solexa sequencing.
we discovered  <dig> novel mirnas families and also detected several isomirnas in soybean.
in our rt-qpcr analysis, we verified that the expression profile of several mirnas varied during abiotic and biotic stresses.
this study has important implications for gene regulation under water deficit and pathogen-infection conditions and also contributes significantly to increase the number of identified mirnas in soybean.
methods
plant materials and treatments
water deficit assay
for water deficit treatment, we used the soybean  merrill) cultivars 'embrapa 48' as a drought-tolerant standard and 'br 16' as a sensitive standard  <cit> .
plants were grown in a greenhouse at embrapa-soybean  using a hydroponic system compound for plastic containers  and an aerated ph  <dig> -balanced nutrient solution.
seeds were pre-germinated on moist filter paper in the dark at 25°c ± 1°c and in 65% ± 5% relative humidity.
plantlets were then placed in polystyrene supports so the roots of the seedlings were completely immersed in the nutrient solution.
each seedling tray was maintained in a greenhouse at 25°c ± 2°c and in 60% ± 5% relative humidity under natural daylight  =  <dig>  ×  <dig> μmoles m- <dig> s- <dig>  equivalent to  <dig>  ×  <dig> lux) for a  <dig> h day.
after  <dig> days, seedlings with the first trifoliate leaf fully developed   <cit>  were submitted to different water-deficit treatments according to martins et al.
<cit> .
the nutrient solution was removed from each plastic container where the roots were kept in the tray in the dark without nutrient solution or water for  <dig> minutes ,  <dig> minutes  and  <dig> minutes .
at the end of each water-deficit period, the roots of the seedlings were immediately frozen in liquid nitrogen and stored at -80°c until rna extraction.
the experimental design was a factorial  with three replicates.
each replicate was composed of five plantlets that were sampled in bulk.
four libraries of small rnas were constructed for the water deficit-stress assays from the following root tissues: 1) roots of drought-sensitive seedlings submitted to  <dig> minutes of stress ); 2) roots of drought-sensitive seedlings submitted to  <dig> minutes and  <dig> minutes of stress ); 3) roots from drought-tolerant seedlings submitted to  <dig> minutes of stress ); and 4) roots of drought-tolerant seedlings submitted to  <dig> minutes and  <dig> minutes of stress ).
asian soybean rust assay
the asr reaction was evaluated in soybean plants in a greenhouse at embrapa-soybean  using a field population of phakopsora pachyrhizi collected from soybean fields in the state of mato grosso, which were maintained for over  <dig> generations on the susceptible cv.
brsms-bacuri.
asr identification was confirmed by its-sequencing analysis as described by silva et al.
<cit> , and it revealed a similarity to the mut zimbabwe isolate.
the soybean plants were grown in a pot-based system.
the 'embrapa 48' genotype was used as a susceptible host plant, which develops a susceptible lesion  after phakopsora pachyrhizi infection.
the 'pi561356' genotype was used as the resistant host, which carries an asr resistance gene mapped onto linkage group g  and develops a reddish-brown  lesion with few or no spores.
urediniospores were collected from infected brsms-bacuri plants in a separate greenhouse by tapping infected leaves over a plastic tray.
the urediniospores were then diluted in distilled water with  <dig> % tween- <dig> to a final concentration of  <dig> ×  <dig> spores/ml.
this spore suspension was sprayed onto three plants per pot at the v <dig> to v <dig> growth stages  <cit> .
a solution without the spores was used for the mock inoculations.
following the asr or mock inoculations, water-misting bags were placed over all plantlets for one day to aid the infection process and to prevent cross-contamination of the mock-infected plants.
the third trifoliolate leaves of six plants were collected  <dig> hours after inoculation  for rna extraction.
the experiment followed a completely randomized design with the three replicates as blocks and a full factorial treatment structure consisting of three treatment factors: hai , genotype , and inoculation type .
for the rust-stress assay, we constructed the other four libraries of small rnas from leaves which were compounded by: 1) leaves of rust-susceptible seedlings with mock inoculation ); 2) leaves of rust-susceptible seedlings with rust-spore inoculation ); 3) leaves of rust-resistant seedlings with mock inoculation ); and 4) leaves of rust-resistant seedlings with rust-spore inoculation ).
rna extraction and sequencing
total rna was isolated from fresh leaves and root materials using trizol , and the rna quality was evaluated by electrophoresis on a 1% agarose gel.
the amount of the rna was verified using a quibit fluorometer and quant-it rna assay kit according to the manufacturer's instructions .
total rna  was sent to fasteris life sciences sa  for processing and sequencing using solexa technology on the illumina genome analyzer gaii.
the libraries were constructed from the eight bar-coded samples  sequenced in a total of two channels.
quality scores were generated from illumina's data analysis pipeline, which are similar to sage phred scores with a maximum value of  <dig>  quality scores are based on the relative confidence of base calls using elements of cluster generation and image quality.
briefly, the processing by illumina for the mirna analyses consisted of the following successive steps: acrylamide gel purification of the rna bands corresponding to the size range 20- <dig> nt, ligation of the 3' and 5' adapters to the rna in two separate subsequent steps each followed by acrylamide gel purification,  cdna synthesis followed by acrylamide gel purification, and a final step of pcr amplification to generate cdna colonies template library for illumina sequencing.
after removing the adapter sequences, the sequences were trimmed into different read lengths from  <dig> to  <dig> nt for further analysis.
prediction of mirnas
the reads were grouped into unique sequences, and the read counts were calculated for each library.
the sequences that presented low read counts  were discarded from the final list of unique sequences, which are referred to as a tag.
the sequences were mapped into the soybean genome  assembly using the soap program  <cit> , which returns information concerning the alignment position, chromosome number and strand.
no mismatches were allowed in the alignments.
the tag alignment position's upstream and downstream genomic sequences  were extracted from the genome assembly using homemade perl scripts.
these genomic regions were then aligned against the reverse complement of its respective tag  using the smith-waterman algorithm  <cit> .
to ensure that these pre-mirna sequences could be precisely processed into mature mirna, the candidates were examined according the following criteria  <cit> : i) the mirna and anti-sense mirna should derive from the opposite stem-arms and must be entirely within the arm of the hairpin; ii) base-pairing between the mirna and anti-sense mirna were restricted to four or fewer mismatches; and iii) the frequency of asymmetric bulges was restricted to less than one and the size should be less than two bases.
the genomic regions that were not possible to align the tag and rc-tag were discarded.
finally, the genomic regions that were limited between the alignment positions of the tag and rc-tag were considered as pre-microrna candidates.
from all the pre-microrna candidate sequences that we selected, only the ones with no more than five matches to the soybean genome were selected for analyzing the secondary structure using the rna-folding program mfold  <cit> .
if a perfect stem-loop structure was formed, the small rna sequence was at one arm of the stem, and the respective anti-sense sequence was at the opposite arm; then, the small rna was consisted as a novel soybean mirna.
mirna validation and expression analysis by rt-qpcr
to validate predicted new mirnas, rt-qpcr in respect to eleven mirnas was performed to examine their expression across the eight different libraries.
from those, six were new mirnas belonging to conserved soybean mirnas families ; one new mirna pertencing to conserved mirnas families in other plants species ; and four were mirnas belonging to novel mirnas families .
the forward mirnas primers were designed based on the full mirnas sequence, and the reverse primer was the universal reverse primer for mirna  <cit> .
the stem-loop primer, used for mirna cdna synthesis, was designed according to cheng et al.
<cit> .
the stem-loop sequence consisted of  <dig> conserved and six variable nucleotides that were specific to the 3' end of the mirna sequence .
the rt-qpcr was performed in an abi  <dig> real-time pcr system  using sybr green i  to detect double-stranded cdna synthesis.
reactions were completed in a volume of  <dig> μl containing  <dig> μl of diluted cdna , 1x sybr green i ,  <dig>  mm dntp, 1x pcr buffer,  <dig> mm mgcl <dig>   <dig>  u platinum taq dna polymerase  and  <dig> nm of each reverse and forward primer.
the universal reverse primer  was used in all rt-qpcr reactions.
samples were analyzed in biological triplicate in a 96-well plate, and a no-template control was included.
we used mir156b , mir172ab  and mir1520d  as reference genes, which has been demonstrated as optimal normalizers for water deficit and rust-stress analysis in glycine max  <cit> .
the conditions were set as the following: an initial polymerase activation step for  <dig> minutes at 94°c,  <dig> cycles for  <dig> seconds at 94°c for denaturation,  <dig> seconds at 60°c for annealing and  <dig> seconds at 72°c for elongation.
a melting curve analysis was programmed at the end of the pcr run over the range 65- <dig>  increasing the temperature stepwise by  <dig> °c.
threshold and baselines were manually determined using the abi  <dig> real-time pcr system sds software v <dig> .
to calculate the relative expression of the mirnas, we used the 2-ΔΔct method.
student's t-test was performed to compare pair-wise differences in expression.
the parameters of two-tailed distribution and two samples assuming unequal variances were established.
the means were considered significantly different when p <  <dig> .
prediction of mirna targets
target prediction for mirnas is straightforward because it is assumed that most of them match their targets with almost perfect complementarity  <cit> .
the putative target genes for all mirnas identified were searched for by using the web-based computer psrna target server   <cit>  which can identify putative targets that may be regulated at post-transcriptional or at translational levels.
mature mirna sequences were used as queries to search for potential target mrnas in the glycine max database .
the total scoring for an alignment was calculated based on the mirna length, and the sequences were considered to be mirna targets if the total score were less than  <dig>  points .
results from these analyses were individually inspected on the phytozome, where the loci and protein annotation were obtained.
in order to look for evidences of the predicted targets of the novel identified mirna, we searched for the mirna targets sites in the soybean degradome libraries published by song et al.
<cit>  available under ncbi-geo accession nμ.
gse <dig>  finally, all putative targets regulated by soybean new mirnas which were analyzed by rt-qpcr were subjected to agrigo database to investigate the gene ontology  <cit> .
RESULTS
to identify mirnas from soybean under water deficit and rust stresses, we generated eight libraries of small rnas species.
from these libraries, a total of  <dig> mirnas ranged from  <dig> to  <dig> nt-long sequence sizes were identified .
all pre-mirna sequence candidates that were selected by the parameters stipulated during the mirna prediction and those that had no more than five matches on the soybean genome were folded using the mfold program.
all mirna sequences with the respective precursor sequence originating at a hairpin structure were submitted to the mirbase to determine if they were a new or known mirna.
we separated the results of these mirnas according the following classes: novel mirnas belonging to mirnas families never detected before ; new mirnas belonging to conserved mirna families in other plants species detected for the first time in soybean ; mirnas belonging to conserved mirna soybean families ; different isoforms of new and known mirnas ; and known mirnas already deposited into the mirbase database  .
identification of novel mirnas from soybean
a total of  <dig> new mirnas belonging to  <dig> novel families  were identified by solexa sequencing in libraries from water deficit and rust infections of glycine max.
these families were provisory nominated seq <dig> to seq <dig> .
the precursor mirna sequences varied from  <dig> to  <dig> nt in length.
precursors of these novel mirnas were identified, and they formed proper secondary hairpin structures, with mfes ranging from - <dig>  to - <dig>  kcal/mol .
the most abundant mature mirnas were  <dig> nt in length.
we also evaluated the genomic location of the new mirnas .
of the  <dig> new mirnas genes identified in soybean, around 86%  were located in intergenic regions and the rest were situated inside genes.
the mature mirnas sequences were localized inside the stem-loop sequence with almost half in each arm:  <dig> mirnas were localized in the 3' arm and  <dig> mirnas were in the 5' arm.
more than 63% of the pre-mirna sequences were in the same sense direction  as the soybean genome annotation.
for all  <dig> novel families identified, four were compounded by mirnas provided from two loci, and we detected only one mirna member for the rest.
sense and anti-sense mirnas were detected only in one family, the seq <dig>  and both were nominated according the arm localization .
most of the new mature mirna sequences presented a uracil  as their first nucleotide, which is in agreement with previous results for soybean root sequences  <cit> .
a the number refers to a new family and the letter refers to the new gene in that family.
b cds: codon sequence.
identification of homologues mirnas of other plant species
to determine whether any of the mirnas identified in our libraries were conserved among other plant species, we searched mirbase for homologues.
besides the novel families identified, we also detected  <dig> mirnas belonging to six conserved families in other plants species .
the families mir <dig>  mir <dig>  mir <dig>  mir <dig>  mir <dig> and mir <dig> were detected for the first time in soybean.
for families mir <dig> and mir <dig>  only a single locus was identified, and for mir <dig>  three genes were found.
in two families, mir <dig> and mir <dig>  we detected sense and antisense mirnas .
mir <dig> was only conserved in arabidopsis lyrata and arabidopsis thaliana.
mir <dig> was found in more different plants species than the other families.
it was found in  <dig> species: arabidopsis thaliana, populus trichocarpa, pinus taeda, vitis vinifera, arachis hypogaea, arabidopsis lyrata, citrus sinensis, oryza sativa, saccharum officinarum, zea mays, physcomitrella patens, selaginella moellendorffii, triticum aestivum, sorghum bicolor, brachypodium distachyon, ricinus communis and aquilegia coerulea .
we observed two families  to be conserved between glycine max and glycine soja; however, we expect that more mirna families could be conserved between these species considering that they are closely related.
this low number is probably due to glycine soja showing fewer mirnas identified to date.
arabidopsis lyrata , arabidopsis thaliana , brassica napus , ricinus communis , medicago truncatula , phaseolus vulgaris , arachis hypogaea , glycine soja , aquilegia coerulea , seleginella moellendorffii , physcomitrella patens , pinus taeda , picea abies , populus trichocarpa , citrus sinensis , vitis vinifera , solanum lycopersicum , brachypodium distachyon , hordeum vulgare , oryza sativa , saccharum officinarum , selaginella moellendorffii , sorghum bicolor , triticum aestivum , and zea mays .
identification of conserved soybean mirnas
to identify conserved soybean mirnas, all  <dig> sequences were searched using blastn against the soybean mirnas in mirbase.
we identified  <dig> families of conserved soybean mirnas in our libraries.
only  <dig> mirna soybean genes that were registered in the mirbase were observed .
from the remaining  <dig> mirna genes,  <dig> were mirnas antisense  to the mirnas presents in mirbase , and  <dig> were new members detected from new loci of known families .
of the  <dig> mirnas identified from the opposite strand of previously known mirnas, six were in the 5' arm and six in the 3'arm.
for the  <dig> new members of conserved soybean families,  <dig> mirnas were  <dig> nt in length.
the family with the largest number of new mirna genes  was mir <dig> .
interestingly, in family mir <dig>  we found three new members with sense and antisense mirnas.
also, in mir <dig>  two new genes with sequences originated from both the 3'and 5'arms were identified.
one new gene was detected in mir <dig>  mir <dig>  mir <dig> and mir <dig> with mature sequences originated from both the 3'and 5'arms .
similar to the observation for the novel soybean mirnas , the new genes in these conserved soybean families were compounded for a majority of mature mirnas with a uracil as the first nucleotide in the 5' end.
a the group number refers to:  the new mirnas identified in the conserved soybean mirnas families;  mirnas originated from the opposite arm of mirnas previously identified; and  mirnas registered in the mirbase that were detected in our libraries.
b * mirnas registered in the mirbase database that were not detected in our libraries.
identification of mirnas isoforms
isoforms of micrornas  are a population of known mirna variants.
they are caused by an imprecise or alternative cleavage of dicer during pre-mirna processing  <cit> .
we detected numerous mirnas with additional nucleotides in the 5'or 3' terminus compared to the recorded mature mirnas.
as isomirnas were previously reported in soybean high-throughput sequencing  <cit> , we found  <dig> isomirnas in our libraries .
these isoforms were observed in  <dig> conserved mirna families and in four novel families.
these mirna isoforms occurred in both strands from the 5' or 3' arm.
the conserved mir1507a and mir1507b were found with the most isomirnas detected .
mir1507a showed a variation of three nucleotides in the 5'end and six nucleotides in the 3'end, and mir1507b showed a variation of three and five nucleotides in the 5'and 3' terminal region respectively .
from the novel mirnas identified, the mir-seq <dig> was the read with the most isoforms detected in our sequencing.
this mirna presented a total of  <dig> different sequences with  <dig> varying nucleotides in both the 5'and 3' ends from six fixed nucleotides .
all isoforms and their respective nominated mature mirnas can be found in additional file  <dig>
a the group number refers to:  the mirnas previously identified in other plant species as described in table 2;  the new mirnas identified in the families of conserved mirnas in soybean;  mirnas originated from the opposite arm of mirnas previously identified; and  mirnas registered in the mirbase database that were detected in our libraries.
b sequence conserved between all isoforms and the number of nucleotide variations in each end.
c total number of isoforms  including the typical member for that gene.
validation of mirnas validation and expression profile by rt-qpcr
the stem-loop rt-qpcr was used to validate and measure the expression of the respective mirnas: mir166a-5p, mir166f, mir169f-3p, mir397ab, mir482bd-3p, mir1513c, mir4415b, mir-seq <dig>  mir-seq <dig>  mir-seq <dig> and mir-seq15ab, detected by solexa sequencing.
these mirnas were validated in all genotypes analyzed during dehydration and rust stress.
the relative expressions of these mirnas in the same eight conditions are shown in figure  <dig>
expression patterns of mirnas during water deficit
to identify water deficit-responsive mirnas, we compared the expression profiles of the  <dig> mirnas in both genotypes before and after stress .
a set of five different mirnas  presented the same behavior during the water deficit stress.
these mirnas were commonly up-regulated during the stress condition in the sensitive genotype, and the opposite occur in the tolerant genotype, where they were down-regulated during the water deficit.
mir-seq <dig> and mir-seq <dig> demonstrated a similar expression across the four conditions.
water deficit significantly increased mir-seq <dig> and mir-seq <dig> expression in the roots compared to the control condition in the sensitive genotype, but both mirnas did not vary in the tolerant plants.
mir166f had its level increased in the sensitive genotype and decreased in the tolerant during the stress compared to the control situation.
interestingly, both genotypes presented the same level during the control condition.
in the sensitive plants, mir-482bd-3p showed a strong decrease when submitted to water deficit, being this low level equally observed in the tolerant genotype during the control condition and decreasing when subjected to stress.
mir4415b presented an effective rise in its expression level during the water deficit in the sensitive plants, and its high level was also observed in the tolerant genotype independent of the condition.
both sensitive and tolerant genotype exhibited the same expression pattern for mir-seq <dig> and its level was increased during the stress compared to the control situation.
expression patterns of mirnas during soybean rust stress
the rt-qpcr analyses of four libraries from the rust assays are shown in figure 1b.
the differential expression analyses revealed that mir166a-5p, mir166f, mir169-3p, mir397ab and mir-seq <dig> were dow-regulated in the susceptible genotype during pathogen infection, and equally expressed in the resistant plants.
the level of mir482bd-3p did not vary significantly between the two different conditions in the susceptible.
however in the resistant genotype, its level is higher during the control condition and decrease with the pathogen attack.
mir1513c presented unchangeable expression in the control and stressed condition for both genotypes, but when we compared the two genotypes; the resistant was down-regulated compared to the susceptible.
a strong decrease was observed for mir4415b in the rust infection when compared with the control in the susceptible plants, and its level is higher in the resistant genotype showing no expression alteration between the conditions.
mir-seq <dig> was down-regulated with respect to the soybean rust infection in both genotypes.
significant difference was observed in mir-seq <dig> expression between the mock and infected plants from the susceptible genotypes.
this mirna presented a low expression level after rust inoculation, and its level decreased in the resistant genotypes remaining similar in the both conditions.
mir-seq15ab expression level was significantly decreased in the rust compared to the mock treatment in the susceptible genotype, the opposite occurs in the resistant genotype, when the control showed a lower level of expression compared to the stressed condition.
target prediction of the soybean mirnas
mirnas suppress gene expression by inhibiting translation, promoting mrna decay or both  <cit> .
target gene identification is challenging due to many factors including the following: binding to their target mrnas by partial complementarity over a short sequence, suppression of an individual target genes is often small, and targeting rules are not completely understood.
we predicted the potential mirnas targets in the psrna database using all identified mirnas as queries.
the results of the analysis were divided into two tables, showing the targets predicted for the novel  and for the conserved mirnas families .
a the data from phytozome version  <dig> .
b pairing obtained in psrnatarget server: "|" indicates a watson-crick base pairing; ":" is a g:u base pairing, and "-" indicates a mismatch.
among the  <dig> novel identified mirnas families, only  <dig> families had targets predicted .
the mirnas families mir-seq <dig>  mir-seq <dig>  mir-seq <dig>  mir-seq <dig>  mir-seq <dig>  mir-seq <dig> and mir-seq <dig> had multiple distinct targets.
mir-seq <dig>  mir-seq <dig> and mir-seq <dig> targeted only one locus.
although, mir-seq <dig>  mir-seq <dig>  mir-seq <dig> and mir-seq <dig> presented several loci as targets, all of them are coding for the same proteins.
fructose-bisphosphate aldolase, lrr -containing proteins, translation elongation factor were predicted to be potential targets of the novel mir-seq <dig> which was investigated by rt-qpcr.
the search for a target of the novel mir-seq <dig>  also analyzed by rt-qpcr, showed a match to glycine max peroxidase precursors mrnas as potential targets.
the oxidoreductase and a transcription regulator factor were predicted to be targeted by mir-seq13; and for the mir-seq <dig> only a translation initiator factor was predicted as a target.
after a comparative analysis of our novel identified mirnas and the degradome libraries of developing soybean seeds it was possible to identify specific sequences in the degradome that corresponds to the downstream sequence of the predicted mirna recognition site.
we identified target sequences to six among the  <dig> novel soybean mirnas .
the list of the  <dig> identified genes is composed by a glucosyl transferase, serine carboxypeptidase, fructose biphosphate aldolase, three leucine-rich repeat protein, two peroxidases and two atp dependent rna helicases .
although many soybean conserved mirnas targets have been predicted and validated by previous studies  <cit> , we also investigated the possible targets for the  <dig> known families of mirnas detected in our sequencing.
of these, only  <dig> families had predicted targets and they are listed in the additional file  <dig>  the conserved mirna families showed multiples targets, however families mir <dig>  mir <dig>  mir <dig>  mir <dig>  mir <dig> and mir <dig> were highly conserved about their targets.
for example, all members from the mir <dig> family  targeted sbp -domain protein.
ap apetala  <dig> transcription factors were targeted by mir <dig> family.
the same occur with mir <dig>  mir <dig>  mir <dig> and mir <dig> families that targeted various genes families as grf  transcription factor, multicopper oxidases, lrr -containing proteins and f-box domain proteins respectively.
these results were already observed across several plant species   <cit> .
gene ontology analysis
the targets of those mirnas which the expression was analyzed by rt-qpcr were also investigated in respect to their gene ontology   <cit> .
among the  <dig> mirna genes, six presented target predictions, which were: mir397ab, mir1513c, mir-seq <dig>  mir-seq <dig>  mir-seq <dig> and mir-seq15ab.
the putative soybean mirnas targets presented diverse functions, however the most representative group was the proteins involved in oxidoreductase activity followed by the proteins involved in the catabolic process .
the result demonstrates that more than 76% of the target proteins are involved in oxidoreductase activity is consistent with the fact that some of the mirnas libraries are originated from stressed plants.
a consequence of many environmental stresses - including water deficit and pathogen attack - is a oxidative stress, i. e. the accumulation of reactive oxygen species , which damage cellular structures  <cit> .
as mirnas mir <dig>  mir-seq <dig> and mir-seq <dig> were predicted to match proteins with oxidative activity, they may act in some level of regulation during water deficit or asr stress.
discussion
the use of deep-sequencing technology was efficient to identify  <dig> mirnas of glycine max.
these mirnas were identified from eight different libraries from precursors with stem-loop secondary structures that also map to the soybean genome .
they were detected from water deficit and rust libraries and were characterized as following: detected for the first time, already detected in some plant species, conserved in soybean, or a variant of a known mirna .
from these analyses, we found  <dig> novel families that had not been detected before, six families that had already been detected in coniferophytes, embryophytes and magnoliophytes , and  <dig> conserved soybean families.
in terms of conserved soybean mirnas, we only detected  <dig> known mirnas in our sequencing.
this small number of known mirna genes detected in our libraries could be due to the two filters used in our processing.
these filters may have missed some known, conserved soybean mirnas because they discarded reads with low frequency and those with more than five matches in the genome.
we detected  <dig> mirnas with additional nucleotides in the 3' or 5' terminus compared to the recorded mature mirna.
these mirna variants  were very common in our population of small, detected rnas.
out of the isomirnas, we observed  <dig> pairs of sense and antisense mirnas.
the duplex presents the antisense strand paired to the corresponding mirna with two nucleotides 3' overhangs .
this shows that the sense and antisense mirnas originated from dcl <dig> processing and supports their validation as true mirnas  <cit> .
in addition, we validated the conserved mirnas in our libraries based on homology to known mirnas in mirbase.
the phylogenetic conservation of mirna sequences is one rule proposed by ambros et al.
<cit>  to characterize mirnas.
in this study, we established new mirnas in soybean that were already detected in other plants species.
however, as opposed to some studies that only blast the candidate to the known mirna mature sequence, our identifications were determined by precursor sequence folding and verification of the genuine hairpin structures.
the complexity of the plant response to biotic and abiotic stresses involves many genes and biochemical and molecular mechanisms, and adaptation to these stresses is achieved through regulating gene expression at the transcriptional and post-transcriptional levels.
with regard to post-transcriptional regulation, mirnas are associated with water deficit response in others plants, but this was the first time that differential expression of these small rnas were observed in soybean during water deficit.
in order to validate  <dig> of the novel mirnas detected in sequencing by the rt-qpcr method, we constructed primers stem-loop and analyzed their expression during abiotic and biotic stresses .
we observed that several mirnas were up-regulated during the water deficit in the sensitive genotype .
however, during the same stress, these mirnas had a different expression in the tolerant genotype.
this distinct mirnas behavior between the two contrasting genotypes under the same conditions could be involved with the drought-tolerance that is observed in the tolerant genotype.
one of these mirnas with this expression pattern is the new mir-seq <dig>  interestingly, mir-seq <dig> was predicted to target peroxidase protein.
as known, stress conditions can produce excess concentrations of reactive oxygen species , resulting in oxidative damage at the cellular level  <cit> .
the increase of this mirna in the sensitive genotype, when subjected to water deficit, could be one of the factors associated with vulnerability of these sensitive plants.
whereas in tolerant genotype during the two conditions, the expression levels of mir-seq <dig> are lower than in the sensitive cultivar during stress.
this situation could indicate that the unchangeable mir-seq <dig> levels in the tolerant genotype may be related to its drought-tolerance capacity.
another interesting point is the expression of a novel mirna mir-seq <dig> that showed increased expression levels during the water deficit stress for both genotypes.
this result allows us to associate this mirna with water deficit stress mechanism independently of the genotype background.
our computational approach showed that one of the loci targeted by mir-seq <dig> corresponds to a fructose-bipfosphato-aldolase enzyme which is a constituent of both the glycolytic/gluconeogenic pathway and the pentose phosphate cycle in plants  <cit> .
therefore increase and/or activation of aldolase appear to be implicated in the plant growth mainly through promotion of the glycolytic pathway function to synthesize atp  <cit> .
since, mir-seq <dig> expression was increased during the stress condition in both genotypes and assuming that it can inhibit or degrade aldolases, it could be associated to metabolism decreasing during water deficit in roots.
plants possess several adaptive traits to support pathogen attacks.
in glycine max, asr is responsible for significant losses in soybean growth areas.
nevertheless, no study investigating mirnas and asr disease had been preformed to date.
to determine if mirnas act as key factors during rust infection or for resistance maintenance, we performed expression analyses with the same  <dig> mirnas during mock and infected conditions in two different genotypes .
in general the mirnas under the fungus infection were down-regulated in the susceptible genotype .
for example, mir-seq <dig>  mir-seq <dig> and mir-seq <dig> which had predicted peroxidases, oxidoreductases and translational initiation factor respectively as targets proteins, were down regulated when infected with asr.
the peroxidases enzymes help to metabolize h2o <dig> in higher plants, and these proteins, as also others proteins with oxidoreductase activity, have already been reported to be up-regulated after pathogen infection and especially after asr  <cit> , indicating a possible involvement of mir-seq <dig> and mir-seq <dig> with the responses to asr infection.
considering, that a translational initiator factor was predicted to be targeted by mir-seq <dig>  we could speculate about the participation of this mirna in the protein synthesis machinery.
in the resistant plants, most of the mirnas analyzed by rt-qpcr  did not vary across the mock and rust infection.
surprisingly, mir-seq <dig> was the unique mirna that was down-regulated during the fungi infection for both genotypes analyzed in our study.
we already mentioned that the mir-seq <dig> had predicted protein target related to metabolism and thus its possible association with water stress.
however mir-seq <dig> also had predicted lrrs -domain target which are known to be present in disease resistance proteins  <cit> .
this suggested a good candidate for the investigation of the mirnas' regulatory mechanisms during asr stress.
although we investigated the expression patterns of some mirnas detected in our sequencing and predicted the target genes that it regulates, additional experimental approaches must be addressed to confirm these hypotheses.
CONCLUSIONS
the present study detected a large number of small rna sequences that were characterized as novel and as already known soybean mirnas.
we grouped some of these unique sequences into  <dig> novel soybean mirnas and further classified several of new members in known families or as new loci in the soybean genome.
validation of new mirna with quantitative rt-qpcr revealed that solexa sequencing is a powerful tool for mirna discovery.
many mirna expression patterns were up- or down-regulated by water deficit and rust stresses, which is an important discovery.
future investigations should use supplementary experimental approaches to verify the targets and to understand the complex gene regulatory network of these mirnas.
this work will contribute to improve systems to support soybean crop production and to mitigating crop losses during biotic or abiotic stresses.
authors' contributions
frk and rm conceived and designed the study.
frk performed the sequence analyses to identify the mirnas and secondary structures and to predict the target genes conceived, executed the rt-qpcr, performed the data management and processing, and wrote the draft manuscript.
rm was the supervisor of this study, provided critical revision, obtained financial support and performed data interpretation.
lfvo contributed to the data assembly, prediction and identification de new mirnas.
lm and ma contributed to the analysis of the mirna secondary structures and processing of the data.
fr, jm, jfb and rsm performed the plant experiments and rna extractions.
aln, fcmg and rva provided the studied material, critically revised the article for important intellectual content and obtained funding.
mfc, gagp and lcn created the perl scripts to identify the micrornas.
mfc participated in writing the methods section.
all authors read and approved the final version of manuscript.
supplementary material
additional file 1
predicted precursor structures of all mirnas identified.
the mature mirnas  and pre-mirna sequences with chromosome and locus information.
the pre-mirna length  and its directional information  or anti-sense  compared to the soybean genome sequence) is provided.
the fold-back structure with respect to the free energy value  was predicted using the mfold program.
click here for file
additional file 2
identified targets of known conserved plant mirnas families.
a the data from phytozome version  <dig> .
b pairing obtained in psrnatarget server: "|" indicates a watson-crick base pairing; ":" is a g:u base pairing, and "-"indicates a mismatch.
click here for file
additional file 3
the soybean transcript loci which were identified as new-mirna families target by degradome sequencing.
the mirna target site is indicated in red and underlined while the degradome sequence is highlighted.
click here for file
acknowledgements
frk was sponsored by a phd grant , ma by a pdj  and rm by a productivity and research grant  from the national council for scientific and technological development .
this work was financially supported by the biotecsur, genosoja consortium  and genoprot .

 small, non-coding rnas have been characterized in plants as important factors involved in gene expression regulation in developmental processes, as well as adaption to biotic and abiotic stress conditions. in general, small rnas are grouped into two major classes: micrornas and short-interfering rnas. these two classes of small rnas cannot be discriminated by either their chemical composition or mechanism of action. however, sirnas and mirnas can be distinguished by their origin, evolutionary conservation and the types of genes that they silence. in this way, mirnas are well differentiated due to some particular characteristics. these characteristics include the following: derived from genomic loci distinct from other recognized genes, processed from transcripts that can form local rna hairpin structures, and usually, mirnas sequences are nearly always conserved in related organisms. in plants, mirna genes are transcribed by rna polymerase ii enzymes generating primary mirna. the pri-mirna forms an imperfect fold-back structure, which is processed into a stem-loop precursor by nuclear rnaseiii-like enzymes called dicer-like proteins. the resulting pre-mirna contains a mirna:mirna* intermediate duplex, formed by a self-complementary fold-back structure. a mature mirna sequence can range from to nucleotides in length and act as a regulatory molecule in post-transcriptional gene silencing by base pairing with target mrnas. this leads to mrna cleavage or translational repression, depending on the degree of complementarity between the mirna and its target transcript. the same mature mirna can also present several variants of their sequence in length. these populations of mirna variants are called isomirnas, which are isoforms of micrornas. they are caused by an imprecise or alternative cleavage of dicer during pre-mirna processing. isomirnas have been recently identified in both plants and animals. the first plant mirnas were described in arabidopsis thaliana and later in other species. currently, mirnas have been reported in plants species, and all of their sequences have been deposited in a publicly-available mirna database, mirbase. several mirnas have been identified in plants, and they are characterized in a wide variety of metabolic and biological processes in plants with important functions in development, phytohormone signaling, flowering and sex determination and responses to biotic and abiotic stresses. in soybean merrill), the major legume crop worldwide, subramanian et al. in identified novel mirna families for the first time. in this study, the role of mirnas in soybean-rhizobial symbiosis was investigated. during that same year, zhang et al. used a comparative genome-based in silico screening of soybean est databases and quantitative pcr to provide evidence for mirnas belonging to families. a second study involving mirnas and soybean root nodules was performed by wang and colleagues. they identified mirnas belonging to mirna families. the identification of nine novel mirnas in wild soybean was also reported by chen et al.. another study looked at four different soybean tissues and identified novel soybean mirnas. recently, song and coworkers identified new mirnas and their related target genes from developing soybean seeds. although these studies resulted in a large number of mirnas identified in soybean, none of them looked at micrornas with respect to biotic and abiotic stresses. drought is the major abiotic stress factor to negatively affect soybean productivity around the world. the impact of limited water during the flower formation can cause shorter flowering periods, and water stress during the later phases of soybean reproductive development has been reported to accelerate senescence, which decreases the duration of the seed-filling period. with regards to biotic stress, asian soybean rust is a foliar disease caused by the fungus phakopsora pachyrhizi sydow & sydow. this pathogen presents a rapid aerial spread and a high capacity to colonize leaf tissue and, to a lesser extent, stem and pods. asr is one of the most severe diseases on the soybean culture, which causes damage between 10% and 90% in the different regions where it has been identified. this disease is the main threat in soybean-producing countries. currently, there are mirnas identified in glycine max; however, none of these mirnas were associated with water deficit or asr stress conditions. we consider that the identification of these mirnas is important to understanding small rna-mediated gene regulation in soybean roots under water deficit stress and in leaves during rust infection. in this context, our goal was to identify new mirnas and to discover those that may be regulated by water deficit and soybean rust stress. using high-throughput sequencing, we constructed four libraries of small rnas from the roots of drought-sensitive and tolerant seedlings in response to control or water deficit conditions. we also constructed four libraries from leaves of rust-susceptible and resistant seedlings with mock and infected conditions. a set of eight small rnas libraries was analyzed from soybean plants. a total of mirnas were detected in solexa sequencing. we discovered novel mirnas families and also detected several isomirnas in soybean. in our rt-qpcr analysis, we verified that the expression profile of several mirnas varied during abiotic and biotic stresses. this study has important implications for gene regulation under water deficit and pathogen-infection conditions and also contributes significantly to increase the number of identified mirnas in soybean. methods plant materials and treatments water deficit assay for water deficit treatment, we used the soybean merrill) cultivars'embrapa 48' as a drought-tolerant standard and'br 16' as a sensitive standard. plants were grown in a greenhouse at embrapa-soybean using a hydroponic system compound for plastic containers and an aerated ph -balanced nutrient solution. seeds were pre-germinated on moist filter paper in the dark at 25°c ± 1°c and in 65% ± 5% relative humidity. plantlets were then placed in polystyrene supports so the roots of the seedlings were completely immersed in the nutrient solution. each seedling tray was maintained in a greenhouse at 25°c ± 2°c and in 60% ± 5% relative humidity under natural daylight = × μmoles m- s- equivalent to × lux) for a h day. after days, seedlings with the first trifoliate leaf fully developed were submitted to different water-deficit treatments according to martins et al.. the nutrient solution was removed from each plastic container where the roots were kept in the tray in the dark without nutrient solution or water for minutes, minutes and minutes. at the end of each water-deficit period, the roots of the seedlings were immediately frozen in liquid nitrogen and stored at -80°c until rna extraction. the experimental design was a factorial with three replicates. each replicate was composed of five plantlets that were sampled in bulk. four libraries of small rnas were constructed for the water deficit-stress assays from the following root tissues: 1) roots of drought-sensitive seedlings submitted to minutes of stress ); 2) roots of drought-sensitive seedlings submitted to minutes and minutes of stress ); 3) roots from drought-tolerant seedlings submitted to minutes of stress ); and 4) roots of drought-tolerant seedlings submitted to minutes and minutes of stress ). asian soybean rust assay the asr reaction was evaluated in soybean plants in a greenhouse at embrapa-soybean using a field population of phakopsora pachyrhizi collected from soybean fields in the state of mato grosso, which were maintained for over generations on the susceptible cv. brsms-bacuri. asr identification was confirmed by its-sequencing analysis as described by silva et al., and it revealed a similarity to the mut zimbabwe isolate. the soybean plants were grown in a pot-based system. the'embrapa 48' genotype was used as a susceptible host plant, which develops a susceptible lesion after phakopsora pachyrhizi infection. the'pi561356' genotype was used as the resistant host, which carries an asr resistance gene mapped onto linkage group g and develops a reddish-brown lesion with few or no spores. urediniospores were collected from infected brsms-bacuri plants in a separate greenhouse by tapping infected leaves over a plastic tray. the urediniospores were then diluted in distilled water with % tween- to a final concentration of × spores/ml. this spore suspension was sprayed onto three plants per pot at the v to v growth stages. a solution without the spores was used for the mock inoculations. following the asr or mock inoculations, water-misting bags were placed over all plantlets for one day to aid the infection process and to prevent cross-contamination of the mock-infected plants. the third trifoliolate leaves of six plants were collected hours after inoculation for rna extraction. the experiment followed a completely randomized design with the three replicates as blocks and a full factorial treatment structure consisting of three treatment factors: hai, genotype, and inoculation type. for the rust-stress assay, we constructed the other four libraries of small rnas from leaves which were compounded by: 1) leaves of rust-susceptible seedlings with mock inoculation ); 2) leaves of rust-susceptible seedlings with rust-spore inoculation ); 3) leaves of rust-resistant seedlings with mock inoculation ); and 4) leaves of rust-resistant seedlings with rust-spore inoculation ). rna extraction and sequencing total rna was isolated from fresh leaves and root materials using trizol, and the rna quality was evaluated by electrophoresis on a 1% agarose gel. the amount of the rna was verified using a quibit fluorometer and quant-it rna assay kit according to the manufacturer's instructions. total rna was sent to fasteris life sciences sa for processing and sequencing using solexa technology on the illumina genome analyzer gaii. the libraries were constructed from the eight bar-coded samples sequenced in a total of two channels. quality scores were generated from illumina's data analysis pipeline, which are similar to sage phred scores with a maximum value of quality scores are based on the relative confidence of base calls using elements of cluster generation and image quality. briefly, the processing by illumina for the mirna analyses consisted of the following successive steps: acrylamide gel purification of the rna bands corresponding to the size range 20- nt, ligation of the 3' and 5' adapters to the rna in two separate subsequent steps each followed by acrylamide gel purification, cdna synthesis followed by acrylamide gel purification, and a final step of pcr amplification to generate cdna colonies template library for illumina sequencing. after removing the adapter sequences, the sequences were trimmed into different read lengths from to nt for further analysis. prediction of mirnas the reads were grouped into unique sequences, and the read counts were calculated for each library. the sequences that presented low read counts were discarded from the final list of unique sequences, which are referred to as a tag. the sequences were mapped into the soybean genome assembly using the soap program, which returns information concerning the alignment position, chromosome number and strand. no mismatches were allowed in the alignments. the tag alignment position's upstream and downstream genomic sequences were extracted from the genome assembly using homemade perl scripts. these genomic regions were then aligned against the reverse complement of its respective tag using the smith-waterman algorithm. to ensure that these pre-mirna sequences could be precisely processed into mature mirna, the candidates were examined according the following criteria: i) the mirna and anti-sense mirna should derive from the opposite stem-arms and must be entirely within the arm of the hairpin; ii) base-pairing between the mirna and anti-sense mirna were restricted to four or fewer mismatches; and iii) the frequency of asymmetric bulges was restricted to less than one and the size should be less than two bases. the genomic regions that were not possible to align the tag and rc-tag were discarded. finally, the genomic regions that were limited between the alignment positions of the tag and rc-tag were considered as pre-microrna candidates. from all the pre-microrna candidate sequences that we selected, only the ones with no more than five matches to the soybean genome were selected for analyzing the secondary structure using the rna-folding program mfold. if a perfect stem-loop structure was formed, the small rna sequence was at one arm of the stem, and the respective anti-sense sequence was at the opposite arm; then, the small rna was consisted as a novel soybean mirna. mirna validation and expression analysis by rt-qpcr to validate predicted new mirnas, rt-qpcr in respect to eleven mirnas was performed to examine their expression across the eight different libraries. from those, six were new mirnas belonging to conserved soybean mirnas families; one new mirna pertencing to conserved mirnas families in other plants species; and four were mirnas belonging to novel mirnas families. the forward mirnas primers were designed based on the full mirnas sequence, and the reverse primer was the universal reverse primer for mirna. the stem-loop primer, used for mirna cdna synthesis, was designed according to cheng et al.. the stem-loop sequence consisted of conserved and six variable nucleotides that were specific to the 3' end of the mirna sequence. the rt-qpcr was performed in an abi real-time pcr system using sybr green i to detect double-stranded cdna synthesis. reactions were completed in a volume of μl containing μl of diluted cdna, 1x sybr green i, mm dntp, 1x pcr buffer, mm mgcl u platinum taq dna polymerase and nm of each reverse and forward primer. the universal reverse primer was used in all rt-qpcr reactions. samples were analyzed in biological triplicate in a 96-well plate, and a no-template control was included. we used mir156b, mir172ab and mir1520d as reference genes, which has been demonstrated as optimal normalizers for water deficit and rust-stress analysis in glycine max. the conditions were set as the following: an initial polymerase activation step for minutes at 94°c, cycles for seconds at 94°c for denaturation, seconds at 60°c for annealing and seconds at 72°c for elongation. a melting curve analysis was programmed at the end of the pcr run over the range 65- increasing the temperature stepwise by °c. threshold and baselines were manually determined using the abi real-time pcr system sds software v. to calculate the relative expression of the mirnas, we used the 2-ΔΔct method. student's t-test was performed to compare pair-wise differences in expression. the parameters of two-tailed distribution and two samples assuming unequal variances were established. the means were considered significantly different when p <. prediction of mirna targets target prediction for mirnas is straightforward because it is assumed that most of them match their targets with almost perfect complementarity. the putative target genes for all mirnas identified were searched for by using the web-based computer psrna target server which can identify putative targets that may be regulated at post-transcriptional or at translational levels. mature mirna sequences were used as queries to search for potential target mrnas in the glycine max database. the total scoring for an alignment was calculated based on the mirna length, and the sequences were considered to be mirna targets if the total score were less than points. results from these analyses were individually inspected on the phytozome, where the loci and protein annotation were obtained. in order to look for evidences of the predicted targets of the novel identified mirna, we searched for the mirna targets sites in the soybean degradome libraries published by song et al. available under ncbi-geo accession nμ. gse finally, all putative targets regulated by soybean new mirnas which were analyzed by rt-qpcr were subjected to agrigo database to investigate the gene ontology.  to identify mirnas from soybean under water deficit and rust stresses, we generated eight libraries of small rnas species. from these libraries, a total of mirnas ranged from to nt-long sequence sizes were identified. all pre-mirna sequence candidates that were selected by the parameters stipulated during the mirna prediction and those that had no more than five matches on the soybean genome were folded using the mfold program. all mirna sequences with the respective precursor sequence originating at a hairpin structure were submitted to the mirbase to determine if they were a new or known mirna. we separated the results of these mirnas according the following classes: novel mirnas belonging to mirnas families never detected before; new mirnas belonging to conserved mirna families in other plants species detected for the first time in soybean; mirnas belonging to conserved mirna soybean families; different isoforms of new and known mirnas; and known mirnas already deposited into the mirbase database. identification of novel mirnas from soybean a total of new mirnas belonging to novel families were identified by solexa sequencing in libraries from water deficit and rust infections of glycine max. these families were provisory nominated seq to seq. the precursor mirna sequences varied from to nt in length. precursors of these novel mirnas were identified, and they formed proper secondary hairpin structures, with mfes ranging from - to - kcal/mol. the most abundant mature mirnas were nt in length. we also evaluated the genomic location of the new mirnas. of the new mirnas genes identified in soybean, around 86% were located in intergenic regions and the rest were situated inside genes. the mature mirnas sequences were localized inside the stem-loop sequence with almost half in each arm: mirnas were localized in the 3' arm and mirnas were in the 5' arm. more than 63% of the pre-mirna sequences were in the same sense direction as the soybean genome annotation. for all novel families identified, four were compounded by mirnas provided from two loci, and we detected only one mirna member for the rest. sense and anti-sense mirnas were detected only in one family, the seq and both were nominated according the arm localization. most of the new mature mirna sequences presented a uracil as their first nucleotide, which is in agreement with previous results for soybean root sequences. a the number refers to a new family and the letter refers to the new gene in that family. b cds: codon sequence. identification of homologues mirnas of other plant species to determine whether any of the mirnas identified in our libraries were conserved among other plant species, we searched mirbase for homologues. besides the novel families identified, we also detected mirnas belonging to six conserved families in other plants species. the families mir mir mir mir mir and mir were detected for the first time in soybean. for families mir and mir only a single locus was identified, and for mir three genes were found. in two families, mir and mir we detected sense and antisense mirnas. mir was only conserved in arabidopsis lyrata and arabidopsis thaliana. mir was found in more different plants species than the other families. it was found in species: arabidopsis thaliana, populus trichocarpa, pinus taeda, vitis vinifera, arachis hypogaea, arabidopsis lyrata, citrus sinensis, oryza sativa, saccharum officinarum, zea mays, physcomitrella patens, selaginella moellendorffii, triticum aestivum, sorghum bicolor, brachypodium distachyon, ricinus communis and aquilegia coerulea. we observed two families to be conserved between glycine max and glycine soja; however, we expect that more mirna families could be conserved between these species considering that they are closely related. this low number is probably due to glycine soja showing fewer mirnas identified to date. arabidopsis lyrata, arabidopsis thaliana, brassica napus, ricinus communis, medicago truncatula, phaseolus vulgaris, arachis hypogaea, glycine soja, aquilegia coerulea, seleginella moellendorffii, physcomitrella patens, pinus taeda, picea abies, populus trichocarpa, citrus sinensis, vitis vinifera, solanum lycopersicum, brachypodium distachyon, hordeum vulgare, oryza sativa, saccharum officinarum, selaginella moellendorffii, sorghum bicolor, triticum aestivum, and zea mays. identification of conserved soybean mirnas to identify conserved soybean mirnas, all sequences were searched using blastn against the soybean mirnas in mirbase. we identified families of conserved soybean mirnas in our libraries. only mirna soybean genes that were registered in the mirbase were observed. from the remaining mirna genes, were mirnas antisense to the mirnas presents in mirbase, and were new members detected from new loci of known families. of the mirnas identified from the opposite strand of previously known mirnas, six were in the 5' arm and six in the 3'arm. for the new members of conserved soybean families, mirnas were nt in length. the family with the largest number of new mirna genes was mir. interestingly, in family mir we found three new members with sense and antisense mirnas. also, in mir two new genes with sequences originated from both the 3'and 5'arms were identified. one new gene was detected in mir mir mir and mir with mature sequences originated from both the 3'and 5'arms. similar to the observation for the novel soybean mirnas, the new genes in these conserved soybean families were compounded for a majority of mature mirnas with a uracil as the first nucleotide in the 5' end. a the group number refers to: the new mirnas identified in the conserved soybean mirnas families; mirnas originated from the opposite arm of mirnas previously identified; and mirnas registered in the mirbase that were detected in our libraries. b * mirnas registered in the mirbase database that were not detected in our libraries. identification of mirnas isoforms isoforms of micrornas are a population of known mirna variants. they are caused by an imprecise or alternative cleavage of dicer during pre-mirna processing. we detected numerous mirnas with additional nucleotides in the 5'or 3' terminus compared to the recorded mature mirnas. as isomirnas were previously reported in soybean high-throughput sequencing, we found isomirnas in our libraries. these isoforms were observed in conserved mirna families and in four novel families. these mirna isoforms occurred in both strands from the 5' or 3' arm. the conserved mir1507a and mir1507b were found with the most isomirnas detected. mir1507a showed a variation of three nucleotides in the 5'end and six nucleotides in the 3'end, and mir1507b showed a variation of three and five nucleotides in the 5'and 3' terminal region respectively. from the novel mirnas identified, the mir-seq was the read with the most isoforms detected in our sequencing. this mirna presented a total of different sequences with varying nucleotides in both the 5'and 3' ends from six fixed nucleotides. all isoforms and their respective nominated mature mirnas can be found in additional file a the group number refers to: the mirnas previously identified in other plant species as described in table 2; the new mirnas identified in the families of conserved mirnas in soybean; mirnas originated from the opposite arm of mirnas previously identified; and mirnas registered in the mirbase database that were detected in our libraries. b sequence conserved between all isoforms and the number of nucleotide variations in each end. c total number of isoforms including the typical member for that gene. validation of mirnas validation and expression profile by rt-qpcr the stem-loop rt-qpcr was used to validate and measure the expression of the respective mirnas: mir166a-5p, mir166f, mir169f-3p, mir397ab, mir482bd-3p, mir1513c, mir4415b, mir-seq mir-seq mir-seq and mir-seq15ab, detected by solexa sequencing. these mirnas were validated in all genotypes analyzed during dehydration and rust stress. the relative expressions of these mirnas in the same eight conditions are shown in figure expression patterns of mirnas during water deficit to identify water deficit-responsive mirnas, we compared the expression profiles of the mirnas in both genotypes before and after stress. a set of five different mirnas presented the same behavior during the water deficit stress. these mirnas were commonly up-regulated during the stress condition in the sensitive genotype, and the opposite occur in the tolerant genotype, where they were down-regulated during the water deficit. mir-seq and mir-seq demonstrated a similar expression across the four conditions. water deficit significantly increased mir-seq and mir-seq expression in the roots compared to the control condition in the sensitive genotype, but both mirnas did not vary in the tolerant plants. mir166f had its level increased in the sensitive genotype and decreased in the tolerant during the stress compared to the control situation. interestingly, both genotypes presented the same level during the control condition. in the sensitive plants, mir-482bd-3p showed a strong decrease when submitted to water deficit, being this low level equally observed in the tolerant genotype during the control condition and decreasing when subjected to stress. mir4415b presented an effective rise in its expression level during the water deficit in the sensitive plants, and its high level was also observed in the tolerant genotype independent of the condition. both sensitive and tolerant genotype exhibited the same expression pattern for mir-seq and its level was increased during the stress compared to the control situation. expression patterns of mirnas during soybean rust stress the rt-qpcr analyses of four libraries from the rust assays are shown in figure 1b. the differential expression analyses revealed that mir166a-5p, mir166f, mir169-3p, mir397ab and mir-seq were dow-regulated in the susceptible genotype during pathogen infection, and equally expressed in the resistant plants. the level of mir482bd-3p did not vary significantly between the two different conditions in the susceptible. however in the resistant genotype, its level is higher during the control condition and decrease with the pathogen attack. mir1513c presented unchangeable expression in the control and stressed condition for both genotypes, but when we compared the two genotypes; the resistant was down-regulated compared to the susceptible. a strong decrease was observed for mir4415b in the rust infection when compared with the control in the susceptible plants, and its level is higher in the resistant genotype showing no expression alteration between the conditions. mir-seq was down-regulated with respect to the soybean rust infection in both genotypes. significant difference was observed in mir-seq expression between the mock and infected plants from the susceptible genotypes. this mirna presented a low expression level after rust inoculation, and its level decreased in the resistant genotypes remaining similar in the both conditions. mir-seq15ab expression level was significantly decreased in the rust compared to the mock treatment in the susceptible genotype, the opposite occurs in the resistant genotype, when the control showed a lower level of expression compared to the stressed condition. target prediction of the soybean mirnas mirnas suppress gene expression by inhibiting translation, promoting mrna decay or both. target gene identification is challenging due to many factors including the following: binding to their target mrnas by partial complementarity over a short sequence, suppression of an individual target genes is often small, and targeting rules are not completely understood. we predicted the potential mirnas targets in the psrna database using all identified mirnas as queries. the results of the analysis were divided into two tables, showing the targets predicted for the novel and for the conserved mirnas families. a the data from phytozome version. b pairing obtained in psrnatarget server: "|" indicates a watson-crick base pairing; ":" is a g:u base pairing, and "-" indicates a mismatch. among the novel identified mirnas families, only families had targets predicted. the mirnas families mir-seq mir-seq mir-seq mir-seq mir-seq mir-seq and mir-seq had multiple distinct targets. mir-seq mir-seq and mir-seq targeted only one locus. although, mir-seq mir-seq mir-seq and mir-seq presented several loci as targets, all of them are coding for the same proteins. fructose-bisphosphate aldolase, lrr -containing proteins, translation elongation factor were predicted to be potential targets of the novel mir-seq which was investigated by rt-qpcr. the search for a target of the novel mir-seq also analyzed by rt-qpcr, showed a match to glycine max peroxidase precursors mrnas as potential targets. the oxidoreductase and a transcription regulator factor were predicted to be targeted by mir-seq13; and for the mir-seq only a translation initiator factor was predicted as a target. after a comparative analysis of our novel identified mirnas and the degradome libraries of developing soybean seeds it was possible to identify specific sequences in the degradome that corresponds to the downstream sequence of the predicted mirna recognition site. we identified target sequences to six among the novel soybean mirnas. the list of the identified genes is composed by a glucosyl transferase, serine carboxypeptidase, fructose biphosphate aldolase, three leucine-rich repeat protein, two peroxidases and two atp dependent rna helicases. although many soybean conserved mirnas targets have been predicted and validated by previous studies, we also investigated the possible targets for the known families of mirnas detected in our sequencing. of these, only families had predicted targets and they are listed in the additional file the conserved mirna families showed multiples targets, however families mir mir mir mir mir and mir were highly conserved about their targets. for example, all members from the mir family targeted sbp -domain protein. ap apetala transcription factors were targeted by mir family. the same occur with mir mir mir and mir families that targeted various genes families as grf transcription factor, multicopper oxidases, lrr -containing proteins and f-box domain proteins respectively. these results were already observed across several plant species. gene ontology analysis the targets of those mirnas which the expression was analyzed by rt-qpcr were also investigated in respect to their gene ontology. among the mirna genes, six presented target predictions, which were: mir397ab, mir1513c, mir-seq mir-seq mir-seq and mir-seq15ab. the putative soybean mirnas targets presented diverse functions, however the most representative group was the proteins involved in oxidoreductase activity followed by the proteins involved in the catabolic process. the result demonstrates that more than 76% of the target proteins are involved in oxidoreductase activity is consistent with the fact that some of the mirnas libraries are originated from stressed plants. a consequence of many environmental stresses - including water deficit and pathogen attack - is a oxidative stress, i. e. the accumulation of reactive oxygen species, which damage cellular structures. as mirnas mir mir-seq and mir-seq were predicted to match proteins with oxidative activity, they may act in some level of regulation during water deficit or asr stress. discussion the use of deep-sequencing technology was efficient to identify mirnas of glycine max. these mirnas were identified from eight different libraries from precursors with stem-loop secondary structures that also map to the soybean genome. they were detected from water deficit and rust libraries and were characterized as following: detected for the first time, already detected in some plant species, conserved in soybean, or a variant of a known mirna. from these analyses, we found novel families that had not been detected before, six families that had already been detected in coniferophytes, embryophytes and magnoliophytes, and conserved soybean families. in terms of conserved soybean mirnas, we only detected known mirnas in our sequencing. this small number of known mirna genes detected in our libraries could be due to the two filters used in our processing. these filters may have missed some known, conserved soybean mirnas because they discarded reads with low frequency and those with more than five matches in the genome. we detected mirnas with additional nucleotides in the 3' or 5' terminus compared to the recorded mature mirna. these mirna variants were very common in our population of small, detected rnas. out of the isomirnas, we observed pairs of sense and antisense mirnas. the duplex presents the antisense strand paired to the corresponding mirna with two nucleotides 3' overhangs. this shows that the sense and antisense mirnas originated from dcl processing and supports their validation as true mirnas. in addition, we validated the conserved mirnas in our libraries based on homology to known mirnas in mirbase. the phylogenetic conservation of mirna sequences is one rule proposed by ambros et al. to characterize mirnas. in this study, we established new mirnas in soybean that were already detected in other plants species. however, as opposed to some studies that only blast the candidate to the known mirna mature sequence, our identifications were determined by precursor sequence folding and verification of the genuine hairpin structures. the complexity of the plant response to biotic and abiotic stresses involves many genes and biochemical and molecular mechanisms, and adaptation to these stresses is achieved through regulating gene expression at the transcriptional and post-transcriptional levels. with regard to post-transcriptional regulation, mirnas are associated with water deficit response in others plants, but this was the first time that differential expression of these small rnas were observed in soybean during water deficit. in order to validate of the novel mirnas detected in sequencing by the rt-qpcr method, we constructed primers stem-loop and analyzed their expression during abiotic and biotic stresses. we observed that several mirnas were up-regulated during the water deficit in the sensitive genotype. however, during the same stress, these mirnas had a different expression in the tolerant genotype. this distinct mirnas behavior between the two contrasting genotypes under the same conditions could be involved with the drought-tolerance that is observed in the tolerant genotype. one of these mirnas with this expression pattern is the new mir-seq interestingly, mir-seq was predicted to target peroxidase protein. as known, stress conditions can produce excess concentrations of reactive oxygen species, resulting in oxidative damage at the cellular level. the increase of this mirna in the sensitive genotype, when subjected to water deficit, could be one of the factors associated with vulnerability of these sensitive plants. whereas in tolerant genotype during the two conditions, the expression levels of mir-seq are lower than in the sensitive cultivar during stress. this situation could indicate that the unchangeable mir-seq levels in the tolerant genotype may be related to its drought-tolerance capacity. another interesting point is the expression of a novel mirna mir-seq that showed increased expression levels during the water deficit stress for both genotypes. this result allows us to associate this mirna with water deficit stress mechanism independently of the genotype background. our computational approach showed that one of the loci targeted by mir-seq corresponds to a fructose-bipfosphato-aldolase enzyme which is a constituent of both the glycolytic/gluconeogenic pathway and the pentose phosphate cycle in plants. therefore increase and/or activation of aldolase appear to be implicated in the plant growth mainly through promotion of the glycolytic pathway function to synthesize atp. since, mir-seq expression was increased during the stress condition in both genotypes and assuming that it can inhibit or degrade aldolases, it could be associated to metabolism decreasing during water deficit in roots. plants possess several adaptive traits to support pathogen attacks. in glycine max, asr is responsible for significant losses in soybean growth areas. nevertheless, no study investigating mirnas and asr disease had been preformed to date. to determine if mirnas act as key factors during rust infection or for resistance maintenance, we performed expression analyses with the same mirnas during mock and infected conditions in two different genotypes. in general the mirnas under the fungus infection were down-regulated in the susceptible genotype. for example, mir-seq mir-seq and mir-seq which had predicted peroxidases, oxidoreductases and translational initiation factor respectively as targets proteins, were down regulated when infected with asr. the peroxidases enzymes help to metabolize h2o in higher plants, and these proteins, as also others proteins with oxidoreductase activity, have already been reported to be up-regulated after pathogen infection and especially after asr, indicating a possible involvement of mir-seq and mir-seq with the responses to asr infection. considering, that a translational initiator factor was predicted to be targeted by mir-seq we could speculate about the participation of this mirna in the protein synthesis machinery. in the resistant plants, most of the mirnas analyzed by rt-qpcr did not vary across the mock and rust infection. surprisingly, mir-seq was the unique mirna that was down-regulated during the fungi infection for both genotypes analyzed in our study. we already mentioned that the mir-seq had predicted protein target related to metabolism and thus its possible association with water stress. however mir-seq also had predicted lrrs -domain target which are known to be present in disease resistance proteins. this suggested a good candidate for the investigation of the mirnas' regulatory mechanisms during asr stress. although we investigated the expression patterns of some mirnas detected in our sequencing and predicted the target genes that it regulates, additional experimental approaches must be addressed to confirm these hypotheses.  the present study detected a large number of small rna sequences that were characterized as novel and as already known soybean mirnas. we grouped some of these unique sequences into novel soybean mirnas and further classified several of new members in known families or as new loci in the soybean genome. validation of new mirna with quantitative rt-qpcr revealed that solexa sequencing is a powerful tool for mirna discovery. many mirna expression patterns were up- or down-regulated by water deficit and rust stresses, which is an important discovery. future investigations should use supplementary experimental approaches to verify the targets and to understand the complex gene regulatory network of these mirnas. this work will contribute to improve systems to support soybean crop production and to mitigating crop losses during biotic or abiotic stresses. authors' contributions frk and rm conceived and designed the study. frk performed the sequence analyses to identify the mirnas and secondary structures and to predict the target genes conceived, executed the rt-qpcr, performed the data management and processing, and wrote the draft manuscript. rm was the supervisor of this study, provided critical revision, obtained financial support and performed data interpretation. lfvo contributed to the data assembly, prediction and identification de new mirnas. lm and ma contributed to the analysis of the mirna secondary structures and processing of the data. fr, jm, jfb and rsm performed the plant experiments and rna extractions. aln, fcmg and rva provided the studied material, critically revised the article for important intellectual content and obtained funding. mfc, gagp and lcn created the perl scripts to identify the micrornas. mfc participated in writing the methods section. all authors read and approved the final version of manuscript. supplementary material additional file 1 predicted precursor structures of all mirnas identified. the mature mirnas and pre-mirna sequences with chromosome and locus information. the pre-mirna length and its directional information or anti-sense compared to the soybean genome sequence) is provided. the fold-back structure with respect to the free energy value was predicted using the mfold program. click here for file additional file 2 identified targets of known conserved plant mirnas families. a the data from phytozome version. b pairing obtained in psrnatarget server: "|" indicates a watson-crick base pairing; ":" is a g:u base pairing, and "-"indicates a mismatch. click here for file additional file 3 the soybean transcript loci which were identified as new-mirna families target by degradome sequencing. the mirna target site is indicated in red and underlined while the degradome sequence is highlighted. click here for file acknowledgements frk was sponsored by a phd grant, ma by a pdj and rm by a productivity and research grant from the national council for scientific and technological development. this work was financially supported by the biotecsur, genosoja consortium and genoprot.
most mirnas were down regulated during the fungus infection in the susceptible genotype; however, in the resistant genotype, most mirnas did not vary during rust attack. among the mirnas analyzed, all showed different expression profiles during biotic and abiotic stresses to soybean. in order to identify new mirnas in soybean and to verify those that are possibly water deficit and rust-stress regulated, eight libraries of small rnas were constructed and submitted to solexa sequencing. a prediction of the putative targets was carried out for conserved and novel mirnas families. validation of our results with quantitative rt-qpcr revealed that solexa sequencing is a powerful tool for mirna discovery. from this total, we identified families of novel mirnas that had not been reported before, six families of conserved mirnas that exist in other plants species, and families previously reported in soybean. small rnas are key regulators of gene expression that guide both transcriptional and post-transcriptional silencing mechanisms in eukaryotes. the identification of differentially expressed plant mirnas provides molecular evidence for the possible involvement of mirnas in the process of water deficit- and rust-stress responses. however, for the tolerant genotype, most of the mirnas were down regulated. the majority of mirnas were up-regulated during water deficit stress in the sensitive plants. to validate novel mirnas, we performed rt-qpcr across the eight different libraries. current studies have demonstrated that micrornas act in several plant pathways associated with tissue proliferation, differentiation, and development and in response to abiotic and biotic stresses. the libraries were developed from drought-sensitive and tolerant seedlings and rust-susceptible and resistant soybeans with or without stressors.
 small rnas are key regulators of gene expression that guide both transcriptional and post-transcriptional silencing mechanisms in eukaryotes. current studies have demonstrated that micrornas act in several plant pathways associated with tissue proliferation, differentiation, and development and in response to abiotic and biotic stresses. in order to identify new mirnas in soybean and to verify those that are possibly water deficit and rust-stress regulated, eight libraries of small rnas were constructed and submitted to solexa sequencing.  the libraries were developed from drought-sensitive and tolerant seedlings and rust-susceptible and resistant soybeans with or without stressors. sequencing the library and subsequent analyses detected mirnas. from this total, we identified families of novel mirnas that had not been reported before, six families of conserved mirnas that exist in other plants species, and families previously reported in soybean. we also observed the presence of several isomirnas during our analyses. to validate novel mirnas, we performed rt-qpcr across the eight different libraries. among the mirnas analyzed, all showed different expression profiles during biotic and abiotic stresses to soybean. the majority of mirnas were up-regulated during water deficit stress in the sensitive plants. however, for the tolerant genotype, most of the mirnas were down regulated. the pattern of mirnas expression was also different for the distinct genotypes submitted to the pathogen stress. most mirnas were down regulated during the fungus infection in the susceptible genotype; however, in the resistant genotype, most mirnas did not vary during rust attack. a prediction of the putative targets was carried out for conserved and novel mirnas families.  validation of our results with quantitative rt-qpcr revealed that solexa sequencing is a powerful tool for mirna discovery. the identification of differentially expressed plant mirnas provides molecular evidence for the possible involvement of mirnas in the process of water deficit- and rust-stress responses. 