BACKGROUND
micrornas  are a class of endogenous, single-stranded, small regulatory rna molecules, which broadly exist in diverse eukaryotes.
in plants, mirnas range from  <dig> to 24 nt in length and regulate the expression of target genes mainly at post-transcriptional levels  <cit> .
plant mirnas are derived from their congate mirna genes that are mostly located in intergenic regions, sometimes in intragenic regions, of the plant genome .
generally, a mir gene is transcribed by polymerase ii into a capped and polyadenylated primary mirna   <cit> .
subsequently, dcl protein carries out the cleavage of pri-mirna into the stem-loop precursor , which is further processed by dcl <dig> to generate mature mirna and mirna*.
mature mirna is then loaded into an argonaute protein to form the mirna-induced silencing complex .
finally, miriscs can bind to the transcripts of target genes through perfect or near-perfect pairing with its mature mirna, and guild mrna cleavage or translational inhibition in plants  <cit> .
thus, plant mirnas are thought to mediate almost all plant cellular and metabolic processes via regulating posttranscriptional gene silencing  <cit> .
mir166s are highly conserved in plants, and  <dig> mir <dig> species have been identified in  <dig> plant species according to the mirbase database .
it is generally accepted that conserved mirnas might play crucial roles in regulating fundamentally important biological processes  <cit> .
so, many efforts have been made to elucidate functional roles and regulatory mechanisms of mir <dig> family in different plant species.
to date, mir166s have been found to be involved in the modulation of various developmental processes via negatively mediating their targets, including sam development, organ polarity, seed development, vascular patterning of shoot, root development, and nutrition ion uptake .
additionally, some evidences indicated that mir <dig> family might play crucial roles in response to abiotic and biotic stresses.
for example, mir <dig> up-regulation upon salinity stress and concomitant depression of its targets were observed in andigena potato, suggesting an important network node for modulating gene expression responsive for growth adjustments  <cit> .
similarly, mir <dig> was induced by phytophthora sojae infection in soybean, indicating that it may be implicated in soybean basal defense  <cit> .
plants often harbor a number of mir166s derived from multigene family with the individual loci.
for example, there are  <dig> mir166s in soybean,  <dig> in populus trichocarpa,  <dig> in arabidopsis thaliana,  <dig> in maize, rice, and physcomitrella patens according to mirbase database .
the mir166s in multigene family were found to be highly conserved to target hd-zip iii family genes such as revoluta , phabulosa , phavoluta , corona  and athb <dig> in a broad range of plant species, indicating that they might exhibit high degree of functional redundancy .
nevertheless, emerging evidences indicate that functional specialization exists in mir <dig> family.
firstly, mir166s in different species can be predicted and/or demonstrated to target non-hd-zipiii genes.
for example, rice dof daily fluctuations  <dig>  can be targeted by mir <dig> in rice, thereby regulating nutrient ion uptake and accumulation in rice  <cit> ; while mir166m rather than other mir <dig> members in physcomitrella patens was predicted to target  <dig> non-hd-zipiii transcripts such as tc <dig>  fc <dig> and tc <dig>  <cit> .
these different set of targets suggest functional diversity of mir <dig> family in plant species.
secondly, spatio-temporal expression patterns of mir <dig> genes also reflect their functional specialization  <cit> .
a good example is that only five members of mir <dig> gene family  in arabidopsis were found to be expressed in embryo, which might provide a positional signal from the basal–peripheral region of early embryos, whereas the other four mir166s can not be detected in any stage of embryogenesis  <cit> .
similarly, three mir <dig> genes  were observed to be expressed specifically in the root endodermis and post-embryonic meristem in arabidopsis  <cit> .
additionally, phylogenetic analysis revealed diversification of mir <dig> gene family among arabidopsis, rice and physcomitrella patens, which laid a foundation for further exploration of the evolutionary and functional divergence of plant mirnas  <cit> .
evidently, the existence of multiple copies of mir <dig> genes contributes to the functional diversity and specificity of mir <dig> gene family members in plants.
thus, addressing the evolutionary conservation and diversification of mir <dig> gene family becomes an important step towards understanding their fine-tuned regulatory roles in plant developmental processes and/or resistance to stresses.
soybean harbors the largest mir <dig> family with  <dig> members in plant.
although the co-evolution and characteristics of mirna genes were globally investigated in soybean , the evolutionary conservation and functional diversification of mir <dig> family members in soybean remain poorly understood.
in this study, five novel mir166s in soybean were identified by data mining approach.
subsequently, the genomic organization and evolution of all the  <dig> mir <dig> genes were investigated; and as a result, eight segmental duplications and four tandem duplications were found to have possibly occurred during evolution of mir <dig> gene family in soybean.
furthermore, promoter analysis and target prediction results pointed to functional diversification of soybean mir <dig> gene family.
finally, the expression patterns of the five newly identified mir166s were examined during seed development and in response to abiotic stresses.
these findings laid a foundation for elucidating functional diversification of soybean mir <dig> gene family, especially during seed development or in response to abiotic stresses.
RESULTS
identification of novel mir166s in soybean
to identify novel mir166s in soybean, pre-mir165/ <dig> sequences from soybean, arabidopsis and rice were used as queries to conduct blast search against the soybean genomic database .
subsequently, the matched genomic sequences were analyzed following a series of screening criteria for encoding mirna sequence, and five novel pre-mir166s were identified in soybean.
based on soybean mir <dig> nomenclature  in mirbase , the five newly identified mir166s were orderly named as mir166v-z.
the characteristics of the five newly identified mir166s were summarized in tables  <dig> and  <dig>  all the five newly identified mature mir <dig> sequences were identical to their paralogs in soybean.
also, they were located on the 3’arm of the secondary stem-loop hairpin structure of the pre-mirnas .
the nucleotide lengths of the five newly identified pre-mirnas ranged from  <dig> to 244 nt, which fall right within the normal range of plant mirna precursors   <cit> .
negative minimal folding free energy  is an important criterion for distinguishing the mirna from other types of rna  <cit> .
in the study, all the potential mirna precursors exhibited low mfes ranging from - <dig>  to - <dig>  kcal/mol .
furthermore, to normalize the potential effect of sequence length on mfes, the minimal folding free energy index  was calculated to differentiate mirnas from other rnas.
as shown in table  <dig>  the pre-mirnas had a high mfei  with an average of about  <dig> , suggesting that they are very likely to be mirna precursors.
thus, the size of soybean mir <dig> family was enlarged from  <dig> to  <dig> members .table  <dig> characteristics of the five newly identified mir166s in soybean
alength of mirna; arm = location of mature mirnas on secondary stem-loop structures of pre- mirna sequences; bnumber of mismatch between predicted and homologous mirna
alength of pre-mirna; bnumber of mismatches between predicted mirna and mirna*
fig.
<dig> stem-loop structures of the five newly identified pre-mir166s in soybean.
the panel a-e sequentially correspond to pre-mir166v-z.
the mature mirna portion is highlighted in green bar
evolutionary relationship of pre-mir166s in soybean
to examine the evolutionary relationship among the  <dig> mir166s in soybean, phylogenetic analysis was conducted using the sequences of their precursors and mature mirnas, respectively.
the phylogram of mir <dig> precursors revealed that  <dig> soybean pre-mir166s were grouped into three classes representing similarities and divergence.
as shown in fig.
2a,  <dig> pre-mir166s formed the largest class, and the rest of pre-mir166s were clustered into two separate branches away from the majority of pre-mir166s.
the five newly identified pre-mir166s  were unevenly distributed in the phylogenetic tree.
as indicated in fig.
2a, pre-mir166w was grouped together with pre-mir166b, n ,u, m, while pre-mir166v, pre-mir166x and pre-mir166y formed a branch with mir166e in class i. in contrast, pre-mir166z formed a separate branch with pre-mir166k in class iii, suggesting that the newly identified pre-mir166s might have evolutionary diversification.
it is interesting to note that ten pairs of mir <dig> precursors with high sequence identity , were individually grouped into discrete clades in the phylogenetic tree , including three pairs containing the newly identified pre-mir166s  and seven previously known pre-mir <dig> pairs .
these observations suggest that these pairs of pre-mir166s are evolutionarily closer to each other compared with the rest members in mir <dig> family.fig.
<dig> phylogenetic analysis and alignment of mir <dig> family in soybean.
a phylogenetic analysis of mir <dig> precursors.
b phylogenetic analysis and alignment of mature mir166s.
the five newly identified mir166s in soybean are highlighted with a star, and the identity between paralogous pair is listed
meanwhile, the topological tree of mature mir166s strongly supported the conserved property of mir <dig> family.
as shown in fig.
2b,  <dig> mature mir166s were grouped into two clades, and all the mature mir166s except mir166l showed identical or nearly identical sequence.
to further understand the conservation and divergence of soybean mir166s,  <dig> unique mature mir <dig> sequences were isolated from mirbase  and designated as umir166- <dig> to  <dig>  respectively .
the phylogenetic analysis indicated that the  <dig> umir166s were clustered into four classes representing conservation and divergence .
six umir166s presentative of  <dig> soybean mir166s were grouped in the class iv , while the representative of gma-mir166l  formed a branch with  <dig> umir166s in the class i. the observation indicated that soybean mir166s are highly conserved among plant species.
we further performed kalmogorov-smirnov analysis to assess the percentage identity of soybean mir166s and their precursors, respectively.
as shown in figure s <dig> , ~ <dig>  fraction of mature mir166s had more than 80% sequence identity, whereas ~ <dig>  fraction of pre-mir166s only showed ~20% sequence identity, indicating that mature mir166s are more conserved than their precursors.
it was also observed that each of the pre-mir <dig> paralogous pairs generated identical mature sequence except the pair of mir166h and mir166j, suggesting that they might target the same set of genes in soybean.
chromosomal distributions and duplications of mir166s in soybean
to determine the distribution of mir166s on different chromosomes in soybean, a chromosome map was constructed.
as shown in fig.
<dig>  the  <dig> pre-mir166s are unevenly located on  <dig> different chromosomes in soybean.
chromosome  <dig> has the largest number with four pre-mir166s, followed by chromosome  <dig> and  <dig> with three pre-mir166s, while two mir166s were distributed on chromosome  <dig>   <dig>   <dig>   <dig>   <dig>  and one mir <dig> on chromosome  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig>  further examination revealed that  <dig> pre-mir166s were found in the intergenic regions of the genome.
intriguingly,  <dig> intragenic mir166s were mapped to the region of intron, exon or/and utr in host genes .
for example, pre-mir166p and pre-mir166l were situated in the utr region of the unag kinase gene and anti-sense strand of intron region of the chaperonin gene, respectively.
these results implied that interactions might exist between intragenic mir166s and host genes.fig.
<dig> chromosomal localization and duplication of mir <dig> family genes in soybean.
each colored box represents a chromosome.
the approximate distribution of each soybean mir <dig> gene is marked on the circle with a short black line.
tandem duplication clusters are indicated with star.
colored lines indicate the linkage group with segmental duplicated mir <dig> genes, and segmental duplication regions were determined using the plant genome duplication database
it has been demonstrated previously that approximately 20% of plant mirnas are clustered, and generally contain conserved mirnas of the same family  <cit> .
to estimate cluster numbers in soybean mir <dig> family, a maximal distance of 3 kb between two consecutive mir166s was used as a stringent criterion.
as a result, four mir <dig> clusters on the same dna strand were identified in soybean , and each cluster contained two copies of mir <dig> family members .
the distance between two pre-mir166s within each cluster ranged from  <dig> to 158 nt , suggesting that they might be polycistronic precursors.
to confirm their polycistronic nature, data mining was conducted against soybean expressed sequence tags  database since ests could provide evidence for mir gene expression in different tissues or processes.
as shown in table s <dig> , the corresponding est simultaneously contained the sequences of two pre-mir166s that belong to the same cluster.
these observations indicated that the corresponding mir genes can indeed be transcribed into polycistronic precursors carrying two tandem copies of mir166s.
thus, we tentatively named these genes as mir166e/q, mir166v/o, mir166x/b and mir166y/n, and three of the five newly identified pre-mir166s belonged to polycistronic precursors.
it has been reported that homologous mir gene clusters are associated with dosage effect  <cit> .
thus, existence of homologous mir166s from a polycistronic mir gene implied that they might exert dosage effect on the regulation of target gene expression.
based on coordinates of pre-mir166s or their neighbor genes , we further investigated whether traceable genome duplications contributed to the expansion of mir <dig> gene family in soybean.
eight sets of mir166s were mapped on eight distinct duplicate blocks , including three newly identified mir166-containing blocks  and five previously known mir166-containing blocks .
the analysis suggested that these pairs of mir166s on the same block were possibly derived from segmental duplication events during evolution.
to investigate the selective evolutionary pressure on mir <dig> gene divergence after duplication, the non-synonymous/synonymous substitution ratio  was retrieved for the eight duplicated pairs of mir <dig> genes.
as shown in additional file 7: table s <dig> the ka/ks value of all the duplicated gene pairs are less than  <dig>  suggesting that these genes might have undergone a purifying selection with limited functional divergence after duplication  <cit> .
analyzing cis-regulatory motifs in mir <dig> promoters
to provide clues for understanding the regulatory mechanism of mir <dig> expression, bioinformatics analysis was conducted to identify putative cis-regulatory elements in promoter region, which was defined based on the estimated transcription start site .
as shown in additional file 8: table s <dig> most putative promoters were located within 800 bp upstream region of pre-mirna foldbacks, consistent with previous report that the distance between tss and precursor is less than 1 kb for more than 85% mirnas in arabidopsis  <cit> .
as expected, the core promoter element tata-box was observed in all the mir <dig> promoter regions, while  <dig> out of  <dig> mir <dig> genes contained caat-box, a common cis-acting element in promoter and enhancer regions.
furthermore, specific promoter elements were analyzed for understanding mir <dig> functional roles.
as indicated in additional file 9: table s <dig> the specific promoter elements were classified into eight classes, including light response, hormone response, seed development, leaf development, biosynthesis and metabolism, cell cycle and circadian, biotic or abiotic defense, and others.
the cis-regulatory elements were non-uniformly distributed in different mir <dig> promoters.
light-responsive elements were most enriched in the promoters of mir <dig> family genes , followed by two endosperm-specific regulatory motifs  and four stress-responsive elements such as tc-rich , hse , are  and mbs , suggesting that these cis-elements are fundamental to the expression of mir <dig> family genes.
in contrast, some motifs exist only in three or fewer members of mir <dig> gene family , implying that these motifs might be involved in controlling the specificity of mir <dig> expression.
it has been well demonstrated that mir <dig> family can control the expression of hd-zipiii genes via directly targeting their transcripts for degradation  <cit> .
interestingly, mir166g, mir166r and mir166z contained hd-zipiii binding site in their promoter regions , implying a possible hd-zipiii-mediated feedback regulatory loop for mir <dig> expression in soybean.fig.
<dig> promoter analysis of mir <dig> genes in soybean.
phylogenetic relationship of mir <dig> promoters is shown in left panel.
the  <dig> classes of cis-elements  are responsive to hormones, abiotic stresses, seed development, hd-zipiii binding, biosynthesis and metabolism, cell cycle and circadian, respectively.
if mir <dig> harbors the cis-element in promoter region, the box is labeled with the number  <dig> and highlighted in blue color, otherwise the box is labeled with the number 0
to investigate the evolutionary relationship between the promoters of mir <dig> family genes, phylogenetic analysis was also conducted.
as shown in fig.
<dig>  nine pairs of mir <dig> promoters can individually form a discrete clade in phylogenetic tree  and have similar motif composition in their promoter regions.
expression analysis of mir166s in soybean
to investigate the expression patterns of mir <dig> genes, we mined the publicly available transcript profiling data of soybean tissues at the phytozome database  and found that  <dig> mir <dig> genes showed tissue-specific expression patterns.
based on their expression patterns, the  <dig> mir <dig> genes can be categorized into three groups .
group  <dig> comprises mir166a, mir166j, mir166h and mir166z, which were mainly expressed in root, root hair, nodule and/or shoot apical meristem.
group  <dig> contains mir166g and mir166i with high expression in flowers, shoot apical meristem, leaf and/or nodules, while group  <dig> consists of mir166n/y, mir166e/q and mir166p with high transcript abundance either in seed, pod, flower or shoot apical meristem.
the maximum fragments per kilobase of transcript per million mapped reads  for mir166n/y was higher  compared to the reads for the rest of mir <dig> gene family , suggesting that mir166n/y might perform a major function among mir <dig> family genes in soybean.fig.
<dig> expression analysis of soybean mir <dig> genes in various tissues.
the transcript profiling data of soybean tissues were extracted from the publicly-available phytozome database  for heatmap generation.
the color scale above the heat map indicates gene expression levels, low transcript abundance indicated by green color and high transcript abundance indicated by red color.
maximum fpkm value for each mir <dig> is shown
additionally, blast search against ests in ncbi database revealed that  <dig> pre-mir166s can perfectly match with one or more ests including nine gene-mapped pre-mir166s and three gene-unmapped pre-mir166s .
among them, three ests generated from mir166e/q, mir166b/x, mir166o/v were derived from immature seeds, and one est for mir166j from germinating shoot.
additionally, five ests generating pre-mir166a, i, n, r, y were originated from the tissues under abiotic stresses.
these observations not only confirmed that at least  <dig> pre-mir166s , out of the  <dig> soybean mir166s, are indeed transcribed, but also implied that these mir166s might probably be involved in seed development, germination and/or in response to abiotic stresses.
targets of mir166s in soybean
to obtain a further understanding of biological functions of mir <dig> family in soybean, the transcripts of glycine max from jgi genomic project were utilized as a reference system to predict the targets of soybean mir166s by employing the psrnatarget tool.
using three representative mature mir166s as query sequences,  <dig> target genes were identified.
as shown in fig.
6a,  <dig> mir166s  were predicted to target ten members of the hd-zipiii family genes , whereas mir166h,k,u,z,l might target the transcripts other than hd-zip iii family genes, such as phytoene synthase , ammonium transporter 2-like  and/or two-component response regulator-like aprr <dig>  this results support the previous report that mir <dig> family displays the functional diversity and specificity in plants  <cit> .
among the predicted targets, cleavage of athb14-like transcript  was chosen to be validated using rlm-race approach.
as shown in fig.
6b-c, the product of athb14-like was precisely cleaved at the 10th or 12th position of complementarity from the 5’ end of mir166s, indicating that mir166s can be indeed involved in the regulation of biological processes via targeting hd-zip iii family in soybean.fig.
<dig> predicted targets of mir <dig> family and validation of mir <dig> target gene using rlm-race pcr.
a the alignments of fragments of target mrnas that have complementarity to mir166s are shown.
mir <dig> sequences are displayed to highlight complementarity to target mrnas.
the targets are classified into seven groups .
b rlm-race pcr for athb14-like.
lane  <dig> and m show cleaved product of desired size  and dna ladder, respectively.
c mapping of athb14-like mrna cleavage sites by rlm-race.
the arrows indicate the cleavage sites and the numbers show the frequency of clones sequenced
to understand the evolutionary relationship between mir166s and their target genes, phylogenetic analysis were performed using  <dig> unique target sequences  and  <dig> target genes from diverse plant species .
as shown in fig.
<dig>   <dig> out of  <dig> utss were clustered with hd-zip family in the phylogenetic tree.
it is reasonable since hd-zip iii family is highly conserved targets of mir166s among diverse plant species.
further observation pointed that uts <dig> and uts <dig> were very close to hd-zip iii family genes in soybean including gma-hb14-like, gma-rev-like and gma-hb <dig>  indicating highly conserved binding site of mi166s in soybean hd-zip iii family.
in contrast, two non-hd-zip iii target genes in soybean  had seven closely related utss, which were also close to sbi-calcium channel alpha- <dig> subunit , ptr-hd-zip iii  <dig>  phv-plastoglobulin- <dig> , ptr-atp synthase b chain atpf.
similarly, gma-appr <dig> was grouped with  <dig> utss  and some target genes including csi-resistance protein-like protein , bna-prk, ath-photosystem i reaction center subunit xi .
the observations suggested that non-hd-zip iii target genes are diverse among different plant species.fig.
<dig> phylogenetic analysis of unique target sequences  present in plant species and binding sites of mir <dig> in target genes.
target genes of mir166s in soybean are highlighted with a star.
the details of utss were listed in additional file 10: table s7
to provide clues for their functions during soybean development, the expression patterns of  <dig> annotated target genes of mir166s were extracted from genome-wide transcript profiling data of soybean tissues from the phytozome database.
as shown in fig.
<dig>  the members of the hd-zipiii family were highly expressed in stem, shoot apical meristem and/or pods, whereas low transcript accumulation was observed in roots and leaves.
additionally, the transcript accumulation of aprr <dig>  psy and amt2-like was highest in flower, root and stem, respectively.fig.
<dig> expression patterns of mir <dig> target genes in different soybean tissues.
the transcript profiling data of  <dig> annotated target genes in soybean tissues were extracted from the publicly-available phytozome database
expression pattern of the five newly identified mir <dig> genes during seed development or in response to abiotic stresses
several studies have documented that mir <dig> plays important roles during seed development and in response to abiotic stresses .
to provide clues for their functional roles in soybean, we took the five newly identified pre-mir166s as examples to examine their expression patterns during seed development and in response to cold, drought and salinity stresses.
expression patterns of the five newly identified mir166s were investigated using seeds collected at  <dig>   <dig>   <dig>  65 days after flowering .
as shown in fig.
9a, all the five newly identified mir <dig> genes showed gradual increase in transcript accumulation as seeds developed, and reached the peak at  <dig> daf.
their fold changes were up to  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig> , respectively, compared with the ones at  <dig> daf, suggesting that the five newly identified mir166s might play an important role at late stages of seed development.fig.
<dig> expression analysis of the five newly identified mir <dig> genes in response to seed development and abiotic stresses.
a expression patterns of the  <dig> mir <dig> genes during seed development .
b expression patterns of the  <dig> mir <dig> genes in seedlings exposed to cold stress for  <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and 72 h. c expression patterns of the  <dig> mir <dig> genes in seedlings exposed to drought stress for  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> d. d expression patterns of the  <dig> mir <dig> genes in seedlings exposed to salinity stress for  <dig>   <dig>   <dig>   <dig>   <dig> and 72 h. error bars indicate se of two biological and three technical replicates.
values were normalized against the subi <dig> gene
ten-day-old soybean seedlings were exposed to cold stress at 4 °c for  <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> h, and expression of the five newly identified mir166s were monitored.
as indicated in fig.
9b, the expression of mir166v and mir166z were dramatically increased by  <dig> - and  <dig> -fold, respectively, at 3 h after cold treatment, implying that mir166v and mir166z might respond to cold stress at early stage.
mir166w showed gradual increase in transcript accumulation as the stress prolonged, while mir166x and mir166y exhibited slight decrease at late stage of cold stress.
when soybean seedlings were exposed to drought stress, the expression of mir166v and mir166y was promptly decreased by  <dig>  and  <dig>  folds, respectively, at 2 days after treatment , indicating that they might be relatively sensitive to drought stress.
also, a distinct decrease in transcript accumulation after 10 days of drought treatment was observed for mir166w, mir166x, mir166y and mir166z with  <dig> ,  <dig> ,  <dig>  and  <dig>  fold changes, respectively as compared to control.
when young seedlings were subjected to salt stress, the expressions of all the  <dig> mir <dig> genes were promptly increased at 3 h after salt treatment .
especially, the expression of mir166v and mir166z reached their maximum level at 3 h salt stress with  <dig>  and  <dig>  fold changes, respectively as compared to non-salinity control, suggesting that they might be responsive to salt stress at early stage.
mir166w and mir166y showed gradual increase in transcript accumulation as the stress prolonged, and reached their maximum level at 72 h salt stress with  <dig>  and  <dig>  fold changes, respectively as compared to control.
expression pattern of mir <dig> target genes during seed development or in response to abiotic stresses
to explore the regulatory function of mir166s in soybean, the expression patterns of  <dig> target genes with annotation in phytozome were also analyzed by qrt-pcr during seed development and in response to abiotic stresses.
as shown in fig.
10a, the expressions of all the hd-zip iii genes except glyma05g <dig> were gradually increased as seeds developed, and reached the peak at  <dig> daf.
their fold changes ranged from  <dig>  to  <dig> , compared with the ones at  <dig> daf.
on the contrary, psy was gradually decreased in transcript accumulation as seeds developed, while amt2-like was dramatically decreased by  <dig>  folds at  <dig> daf.fig.
<dig> expression analysis of the  <dig> target genes in response to seed development and abiotic stresses.
expression patterns of target genes during seed development , under cold stress , drought stress  and salinity stress .
stress treatments and seed development stage were used as same as the ones in fig.
<dig>  psy, phytoene synthase; amt2-like, ammonium transporter 2-like.
error bars indicate se of two biological and three technical replicates.
values were normalized against the subi <dig> gene
as indicated in fig.
10b, the expressions of all the hd-zip iii genes and aprr <dig> were obviously decreased by cold stress, and reached their minimal levels at 24 h stress with  <dig> – <dig>  fold changes as compared to control.
in contrast, the transcript accumulations of psy and amt2-like were increased by  <dig>  and  <dig>  folds at  <dig> and 48 h after treatment, respectively.
four athb14-like, rev-like  and psy displayed gradual decrease in transcript accumulation as drought stress prolonged , and their maximal fold changes ranged from  <dig>  to  <dig> .
although atm2-like displayed similar expression level to the control at early stage , the transcript accumulation was dramatically increased from 6 days after treatment until 10 days with  <dig>  fold changes, suggesting that atm2-like might respond to drought stress at middle and late stages.
the decreased transcript accumulation was also observed for the rest of target genes on one or more treatment points as compared to control.
as shown in fig.
10d, all the hd-zip iii genes and aprr <dig> were promptly increased to their peaks in transcript accumulation at 3 h after salt stress with  <dig> – <dig>  fold changes, respectively as compared to control, indicating that they might be responsive to salt stress at early stage.
the expressions of psy and amt2-like were gradually increased, and reached their peaks at 72 h after salt treatment with  <dig> – <dig>  fold changes, respectively as compared to control.
discussion
mir <dig> acts as a highly conserved mirna to be implicated in a wide range of cellular and physiological processes in plants .
based on the data from mirbase, most of plant species harbor multiple mir166s.
it has been accepted that the existence of multiple copies of mir genes within plant species provides one possibility for functional redundancy and specificity  <cit> .
thus, the diversity and evolutionary conservation of mir <dig> family genes are expected to be elucidated for further understanding their fine-tuned regulatory roles in plant developmental processes and/or resistance to stresses.
generally, mirnas are annotated using comparative criteria for both their expression and biogenesis, including size and sequence of mature mirna, phylogenetic conservation, secondary structure of mirna precursor and increased precursor accumulation with decreased dicer function in vivo  <cit> .
in this study, the five newly identified mirnas exhibited the same size and sequence as their paralogs in soybean .
the lengths of the potential mirna precursors range from  <dig> to 244 nt , which is acceptable since the length of plant mirna precursors has been found to be in the range of 55–930 nt   <cit> .
furthermore, their precursors also satisfied the criteria with aspect to both phylogenetic conservation and secondary structure.
it has been proposed that low mfe and high mfei are reliable criteria to distinguish mirnas from all coding and other non-coding rnas, as these parameters, to certain extent, reflect the stability of the prefect or near-perfect secondary hairpin structure of mirna precursor  <cit> .
in the study, mfe values of the five newly identified pre-mir166s ranged from - <dig>  to - <dig>  kcal/mol, which were much lower than folding free energies of trna  and rrna .
at the same time, their individual mfei value was also higher than that of trna, rrna and mrna.
taken together, the five newly identified mirnas in soybean are very likely to be true mirnas.
thus, the study enlarged the size of soybean mir <dig> family from  <dig> to  <dig> members.
it has been proposed that conserved mirna families have been both conserved and diversified during mirna evolution  <cit> .
in this study, our data indicate that soybean mir166s showed evolutionary conservation.
as shown in fig.
2a, ten pre-mir <dig> paralogous pairs with high sequence identity  were individually grouped into a discrete clade in the phylogenetic tree.
further analysis revealed that eight out of ten paralogous pairs were located at corresponding loci on the same duplication block , indicating that they might have arisen from segmental duplications without major sequence change.
additionally, smaller than  <dig> of ka/ks values suggest that these mir <dig> genes likely had undergone through a purifying selection with limited functional divergence after duplication.
it is also worth noting that five paralogous pairs were clustered into discrete clades in phylogenetic tree based on their promoter sequences , supporting the view that gene duplication could happen simultaneously at the coding region and the promoter  <cit> .
evidently, the origins of these soybean mir166s were in line with soybean evolutionary history that soybean had undergone whole-genome duplication event ~ <dig>  million years ago  <cit> .
apart from the evolutionary conservation, soybean mir166s also exhibited evolutionary diversification.
in the study, five pre-mir166s  were separated far from the majority of pre-mir166s , and the identity between the five pre-mir166s and the rest of the pre-mir166s was low  .
these observations suggested that they had diverged from the majority of mir166s during soybean evolution.
mature mir166s are ultimately responsible for functional diversification since they act as functional unit to regulate the expression of their targets.
in the study, all the mir166s except mir166l were highly conserved in mature sequence, and identical mature sequence exists in each of the pre-mir <dig> paralogous pairs with exception of pre-mir166h and pre-mir166j .
consequently, most of soybean mir166s shared the same putative targets  in soybean.
however, mir166z, k, h, u, l are likely to have functional diversity since they are predicted to target genes other than those of the hd-zip iii family, i.e., aprr <dig>  phytoene synthase and ammonium transporter 2-like .
thus, our findings demonstrate that mature mir166s exhibit both functional conservation and diversification.
it has been well-known that mir <dig> plays important roles during seed development.
for example, mir <dig> can repress the seed maturation program during vegetative development via controlling the expression of hd-zip iii family genes such as phb  <cit> ; while it might also be involved in somatic embryo formation via acting as a positional signal from the basal–peripheral region of early embryos  <cit> .
in this study, all the five newly identified mir166s showed gradual increase during seed development, and reached their expression peaks at  <dig> daf, a late stage of soybean seed development .
these expression patterns are consistent with previous report that mir <dig> in brassica napus was strongly accumulated at late stage of seed development  <cit> .
also, it was proposed that the accumulation of mature mir <dig> reached its major peak at late cotyledonary embryo stage in larch, suggesting that it might be associated with cotyledon formation  <cit> .
the expression patterns of the five newly identified mir166s suggested that they might be one of the major contributors to the network controlling seed development, especially seed maturation program.
it has been documented that plant mir166s might be involved in response to diverse abiotic stresses, and main evidences were derived from the altered accumulation pattern of mature mir <dig> under diverse abiotic stresses.
for instance, the accumulation of mature mir <dig> was decreased by drought stress in sorghum bicolor  <cit> , whereas drought stress led to increase of mirn <dig> accumulation in root and leaf of wheat  <cit> .
another kind of examples are that mir <dig> in potato was significantly up-regulated by salinity stress  <cit> , and depression of mir <dig> by small tandem target mimic can lead to improvement of salt tolerance in arabidopsis  <cit> .
in the study, the presence of stress-responsive elements in the promoter regions of mir <dig> genes pointed to their possible roles under various stresses .
furthermore, the altered expression patterns of the five newly identified mir <dig> genes suggest their prominent roles in response to cold, drought and salinity .
for example, mir166v and mir166w harbored the ltr cis-acting element in promoter regions, which can be involved in low-temperature responsiveness .
indeed, mir166v was promptly increased by cold stress with  <dig>  fold changes, while mir166w showed gradual increase in transcript accumulation as cold stress prolonged , implying that mir166v and mir166w can possibly function as a modulator of cold-induced signaling pathways.
these results are consistent with previous reports that mir166e and mir166r were increased by chilling stress in vegetable soybean  <cit> , and that the expression levels of mir166u were obviously increased under cold condition in nitrogen-fixing nodules of soybean  <cit> .
it has been proposed that spatial-temporal or specific expression of mirna isoforms is a crucial determinant of isoform specificity  <cit> .
in the study, the differential expression patterns of  <dig> mir <dig> genes under same stress indicated their functional specificity .
as shown in fig.
7c, mir166v was dramatically decreased at 2 days after drought treatment, whereas mir166w was increased at  <dig>   <dig>  8 days after drought treatment.
thus, the stress-specific properties of mir <dig> genes provided important clues for elucidating their functional specificity.
although five novel mir166s in soybean were differentially regulated during seed development or in response to stresses, their expression patterns were not negatively correlated with the ones of mir <dig> target genes .
for example, the expressions of the hd-zip iii genes were gradually increased as seeds developed.
this is consistent with the previous report that hd-zip iii family in arabidopsis can be positively involved in regulation of seed mature program  <cit> .
however, the five new identified mir166s also showed gradual increase in transcript accumulations.
these results suggested that mir166-mediated silencing of target genes was under sophisticated regulation.
mir <dig> gene family has large number of members in plant species, and at least  <dig> mir <dig> members exist in soybean.
so, we speculated that mir <dig> gene family possibly contributes to the accumulation pattern of target genes in a family member-specific manner.
it has been widely accepted that members of a certain evolutionary branch have potential to share similar functions  <cit> .
in the study, pre-mir166x and pre-mir166y were the most closely related paralogs in the phylogenetic tree , and they also showed  <dig> % sequence identity .
as expected, similar expression patterns were observed for mir166x and mir166y in response to cold, drought and salinity stresses .
additionally, the most closely related pre-mir <dig> paralogs  also shared similar tissue-specific expression pattern in soybean .
thus, it was hypothesized that these most closely related pre-mir <dig> paralogs in soybean perform similar cellular functions during developmental processes or in response to various stresses.
CONCLUSIONS
this study enlarged the size of soybean mir <dig> family from  <dig> to  <dig> members.
a comprehensive analysis indicated that the  <dig> soybean mir166s exhibited evolutionary conservation and diversification.
furthermore, the expression patterns of five newly identified mir166s during seed development and in response to abiotic stresses provided important clues for elucidating the functional specificity and redundancy of mir <dig> family in soybean.
future study will aim at investigating the effect of each mir <dig> gene on seed development and stress tolerance, identifying target genes and dissecting the regulatory mechanisms for each mir <dig>  and exploring combination of multiple mir <dig> genes during seed development and in response to specific stress.

 micrornas are a class of endogenous, single-stranded, small regulatory rna molecules, which broadly exist in diverse eukaryotes. in plants, mirnas range from to 24 nt in length and regulate the expression of target genes mainly at post-transcriptional levels. plant mirnas are derived from their congate mirna genes that are mostly located in intergenic regions, sometimes in intragenic regions, of the plant genome. generally, a mir gene is transcribed by polymerase ii into a capped and polyadenylated primary mirna. subsequently, dcl protein carries out the cleavage of pri-mirna into the stem-loop precursor, which is further processed by dcl to generate mature mirna and mirna*. mature mirna is then loaded into an argonaute protein to form the mirna-induced silencing complex. finally, miriscs can bind to the transcripts of target genes through perfect or near-perfect pairing with its mature mirna, and guild mrna cleavage or translational inhibition in plants. thus, plant mirnas are thought to mediate almost all plant cellular and metabolic processes via regulating posttranscriptional gene silencing. mir166s are highly conserved in plants, and mir species have been identified in plant species according to the mirbase database. it is generally accepted that conserved mirnas might play crucial roles in regulating fundamentally important biological processes. so, many efforts have been made to elucidate functional roles and regulatory mechanisms of mir family in different plant species. to date, mir166s have been found to be involved in the modulation of various developmental processes via negatively mediating their targets, including sam development, organ polarity, seed development, vascular patterning of shoot, root development, and nutrition ion uptake. additionally, some evidences indicated that mir family might play crucial roles in response to abiotic and biotic stresses. for example, mir up-regulation upon salinity stress and concomitant depression of its targets were observed in andigena potato, suggesting an important network node for modulating gene expression responsive for growth adjustments. similarly, mir was induced by phytophthora sojae infection in soybean, indicating that it may be implicated in soybean basal defense. plants often harbor a number of mir166s derived from multigene family with the individual loci. for example, there are mir166s in soybean, in populus trichocarpa, in arabidopsis thaliana, in maize, rice, and physcomitrella patens according to mirbase database. the mir166s in multigene family were found to be highly conserved to target hd-zip iii family genes such as revoluta, phabulosa, phavoluta, corona and athb in a broad range of plant species, indicating that they might exhibit high degree of functional redundancy. nevertheless, emerging evidences indicate that functional specialization exists in mir family. firstly, mir166s in different species can be predicted and/or demonstrated to target non-hd-zipiii genes. for example, rice dof daily fluctuations can be targeted by mir in rice, thereby regulating nutrient ion uptake and accumulation in rice; while mir166m rather than other mir members in physcomitrella patens was predicted to target non-hd-zipiii transcripts such as tc fc and tc. these different set of targets suggest functional diversity of mir family in plant species. secondly, spatio-temporal expression patterns of mir genes also reflect their functional specialization. a good example is that only five members of mir gene family in arabidopsis were found to be expressed in embryo, which might provide a positional signal from the basal–peripheral region of early embryos, whereas the other four mir166s can not be detected in any stage of embryogenesis. similarly, three mir genes were observed to be expressed specifically in the root endodermis and post-embryonic meristem in arabidopsis. additionally, phylogenetic analysis revealed diversification of mir gene family among arabidopsis, rice and physcomitrella patens, which laid a foundation for further exploration of the evolutionary and functional divergence of plant mirnas. evidently, the existence of multiple copies of mir genes contributes to the functional diversity and specificity of mir gene family members in plants. thus, addressing the evolutionary conservation and diversification of mir gene family becomes an important step towards understanding their fine-tuned regulatory roles in plant developmental processes and/or resistance to stresses. soybean harbors the largest mir family with members in plant. although the co-evolution and characteristics of mirna genes were globally investigated in soybean, the evolutionary conservation and functional diversification of mir family members in soybean remain poorly understood. in this study, five novel mir166s in soybean were identified by data mining approach. subsequently, the genomic organization and evolution of all the mir genes were investigated; and as a result, eight segmental duplications and four tandem duplications were found to have possibly occurred during evolution of mir gene family in soybean. furthermore, promoter analysis and target prediction results pointed to functional diversification of soybean mir gene family. finally, the expression patterns of the five newly identified mir166s were examined during seed development and in response to abiotic stresses. these findings laid a foundation for elucidating functional diversification of soybean mir gene family, especially during seed development or in response to abiotic stresses.  identification of novel mir166s in soybean to identify novel mir166s in soybean, pre-mir165/ sequences from soybean, arabidopsis and rice were used as queries to conduct blast search against the soybean genomic database. subsequently, the matched genomic sequences were analyzed following a series of screening criteria for encoding mirna sequence, and five novel pre-mir166s were identified in soybean. based on soybean mir nomenclature in mirbase, the five newly identified mir166s were orderly named as mir166v-z. the characteristics of the five newly identified mir166s were summarized in tables  and all the five newly identified mature mir sequences were identical to their paralogs in soybean. also, they were located on the 3’arm of the secondary stem-loop hairpin structure of the pre-mirnas. the nucleotide lengths of the five newly identified pre-mirnas ranged from to 244 nt, which fall right within the normal range of plant mirna precursors. negative minimal folding free energy is an important criterion for distinguishing the mirna from other types of rna. in the study, all the potential mirna precursors exhibited low mfes ranging from - to -  kcal/mol. furthermore, to normalize the potential effect of sequence length on mfes, the minimal folding free energy index was calculated to differentiate mirnas from other rnas. as shown in table  the pre-mirnas had a high mfei with an average of about, suggesting that they are very likely to be mirna precursors. thus, the size of soybean mir family was enlarged from to members.table characteristics of the five newly identified mir166s in soybean alength of mirna; arm = location of mature mirnas on secondary stem-loop structures of pre- mirna sequences; bnumber of mismatch between predicted and homologous mirna alength of pre-mirna; bnumber of mismatches between predicted mirna and mirna* fig. stem-loop structures of the five newly identified pre-mir166s in soybean. the panel a-e sequentially correspond to pre-mir166v-z. the mature mirna portion is highlighted in green bar evolutionary relationship of pre-mir166s in soybean to examine the evolutionary relationship among the mir166s in soybean, phylogenetic analysis was conducted using the sequences of their precursors and mature mirnas, respectively. the phylogram of mir precursors revealed that soybean pre-mir166s were grouped into three classes representing similarities and divergence. as shown in fig. 2a, pre-mir166s formed the largest class, and the rest of pre-mir166s were clustered into two separate branches away from the majority of pre-mir166s. the five newly identified pre-mir166s were unevenly distributed in the phylogenetic tree. as indicated in fig. 2a, pre-mir166w was grouped together with pre-mir166b, n,u, m, while pre-mir166v, pre-mir166x and pre-mir166y formed a branch with mir166e in class i. in contrast, pre-mir166z formed a separate branch with pre-mir166k in class iii, suggesting that the newly identified pre-mir166s might have evolutionary diversification. it is interesting to note that ten pairs of mir precursors with high sequence identity, were individually grouped into discrete clades in the phylogenetic tree, including three pairs containing the newly identified pre-mir166s and seven previously known pre-mir pairs. these observations suggest that these pairs of pre-mir166s are evolutionarily closer to each other compared with the rest members in mir family.fig. phylogenetic analysis and alignment of mir family in soybean. a phylogenetic analysis of mir precursors. b phylogenetic analysis and alignment of mature mir166s. the five newly identified mir166s in soybean are highlighted with a star, and the identity between paralogous pair is listed meanwhile, the topological tree of mature mir166s strongly supported the conserved property of mir family. as shown in fig. 2b, mature mir166s were grouped into two clades, and all the mature mir166s except mir166l showed identical or nearly identical sequence. to further understand the conservation and divergence of soybean mir166s, unique mature mir sequences were isolated from mirbase and designated as umir166- to respectively. the phylogenetic analysis indicated that the umir166s were clustered into four classes representing conservation and divergence. six umir166s presentative of soybean mir166s were grouped in the class iv, while the representative of gma-mir166l formed a branch with umir166s in the class i. the observation indicated that soybean mir166s are highly conserved among plant species. we further performed kalmogorov-smirnov analysis to assess the percentage identity of soybean mir166s and their precursors, respectively. as shown in figure s, ~ fraction of mature mir166s had more than 80% sequence identity, whereas ~ fraction of pre-mir166s only showed ~20% sequence identity, indicating that mature mir166s are more conserved than their precursors. it was also observed that each of the pre-mir paralogous pairs generated identical mature sequence except the pair of mir166h and mir166j, suggesting that they might target the same set of genes in soybean. chromosomal distributions and duplications of mir166s in soybean to determine the distribution of mir166s on different chromosomes in soybean, a chromosome map was constructed. as shown in fig.  the pre-mir166s are unevenly located on different chromosomes in soybean. chromosome has the largest number with four pre-mir166s, followed by chromosome and with three pre-mir166s, while two mir166s were distributed on chromosome and one mir on chromosome and further examination revealed that pre-mir166s were found in the intergenic regions of the genome. intriguingly, intragenic mir166s were mapped to the region of intron, exon or/and utr in host genes. for example, pre-mir166p and pre-mir166l were situated in the utr region of the unag kinase gene and anti-sense strand of intron region of the chaperonin gene, respectively. these results implied that interactions might exist between intragenic mir166s and host genes.fig. chromosomal localization and duplication of mir family genes in soybean. each colored box represents a chromosome. the approximate distribution of each soybean mir gene is marked on the circle with a short black line. tandem duplication clusters are indicated with star. colored lines indicate the linkage group with segmental duplicated mir genes, and segmental duplication regions were determined using the plant genome duplication database it has been demonstrated previously that approximately 20% of plant mirnas are clustered, and generally contain conserved mirnas of the same family. to estimate cluster numbers in soybean mir family, a maximal distance of 3 kb between two consecutive mir166s was used as a stringent criterion. as a result, four mir clusters on the same dna strand were identified in soybean, and each cluster contained two copies of mir family members. the distance between two pre-mir166s within each cluster ranged from to 158 nt, suggesting that they might be polycistronic precursors. to confirm their polycistronic nature, data mining was conducted against soybean expressed sequence tags database since ests could provide evidence for mir gene expression in different tissues or processes. as shown in table s, the corresponding est simultaneously contained the sequences of two pre-mir166s that belong to the same cluster. these observations indicated that the corresponding mir genes can indeed be transcribed into polycistronic precursors carrying two tandem copies of mir166s. thus, we tentatively named these genes as mir166e/q, mir166v/o, mir166x/b and mir166y/n, and three of the five newly identified pre-mir166s belonged to polycistronic precursors. it has been reported that homologous mir gene clusters are associated with dosage effect. thus, existence of homologous mir166s from a polycistronic mir gene implied that they might exert dosage effect on the regulation of target gene expression. based on coordinates of pre-mir166s or their neighbor genes, we further investigated whether traceable genome duplications contributed to the expansion of mir gene family in soybean. eight sets of mir166s were mapped on eight distinct duplicate blocks, including three newly identified mir166-containing blocks and five previously known mir166-containing blocks. the analysis suggested that these pairs of mir166s on the same block were possibly derived from segmental duplication events during evolution. to investigate the selective evolutionary pressure on mir gene divergence after duplication, the non-synonymous/synonymous substitution ratio was retrieved for the eight duplicated pairs of mir genes. as shown in additional file 7: table s the ka/ks value of all the duplicated gene pairs are less than suggesting that these genes might have undergone a purifying selection with limited functional divergence after duplication. analyzing cis-regulatory motifs in mir promoters to provide clues for understanding the regulatory mechanism of mir expression, bioinformatics analysis was conducted to identify putative cis-regulatory elements in promoter region, which was defined based on the estimated transcription start site. as shown in additional file 8: table s most putative promoters were located within 800 bp upstream region of pre-mirna foldbacks, consistent with previous report that the distance between tss and precursor is less than 1 kb for more than 85% mirnas in arabidopsis. as expected, the core promoter element tata-box was observed in all the mir promoter regions, while out of mir genes contained caat-box, a common cis-acting element in promoter and enhancer regions. furthermore, specific promoter elements were analyzed for understanding mir functional roles. as indicated in additional file 9: table s the specific promoter elements were classified into eight classes, including light response, hormone response, seed development, leaf development, biosynthesis and metabolism, cell cycle and circadian, biotic or abiotic defense, and others. the cis-regulatory elements were non-uniformly distributed in different mir promoters. light-responsive elements were most enriched in the promoters of mir family genes, followed by two endosperm-specific regulatory motifs and four stress-responsive elements such as tc-rich, hse, are and mbs, suggesting that these cis-elements are fundamental to the expression of mir family genes. in contrast, some motifs exist only in three or fewer members of mir gene family, implying that these motifs might be involved in controlling the specificity of mir expression. it has been well demonstrated that mir family can control the expression of hd-zipiii genes via directly targeting their transcripts for degradation. interestingly, mir166g, mir166r and mir166z contained hd-zipiii binding site in their promoter regions, implying a possible hd-zipiii-mediated feedback regulatory loop for mir expression in soybean.fig. promoter analysis of mir genes in soybean. phylogenetic relationship of mir promoters is shown in left panel. the classes of cis-elements are responsive to hormones, abiotic stresses, seed development, hd-zipiii binding, biosynthesis and metabolism, cell cycle and circadian, respectively. if mir harbors the cis-element in promoter region, the box is labeled with the number and highlighted in blue color, otherwise the box is labeled with the number 0 to investigate the evolutionary relationship between the promoters of mir family genes, phylogenetic analysis was also conducted. as shown in fig.  nine pairs of mir promoters can individually form a discrete clade in phylogenetic tree and have similar motif composition in their promoter regions. expression analysis of mir166s in soybean to investigate the expression patterns of mir genes, we mined the publicly available transcript profiling data of soybean tissues at the phytozome database and found that mir genes showed tissue-specific expression patterns. based on their expression patterns, the mir genes can be categorized into three groups. group comprises mir166a, mir166j, mir166h and mir166z, which were mainly expressed in root, root hair, nodule and/or shoot apical meristem. group contains mir166g and mir166i with high expression in flowers, shoot apical meristem, leaf and/or nodules, while group consists of mir166n/y, mir166e/q and mir166p with high transcript abundance either in seed, pod, flower or shoot apical meristem. the maximum fragments per kilobase of transcript per million mapped reads for mir166n/y was higher compared to the reads for the rest of mir gene family, suggesting that mir166n/y might perform a major function among mir family genes in soybean.fig. expression analysis of soybean mir genes in various tissues. the transcript profiling data of soybean tissues were extracted from the publicly-available phytozome database for heatmap generation. the color scale above the heat map indicates gene expression levels, low transcript abundance indicated by green color and high transcript abundance indicated by red color. maximum fpkm value for each mir is shown additionally, blast search against ests in ncbi database revealed that pre-mir166s can perfectly match with one or more ests including nine gene-mapped pre-mir166s and three gene-unmapped pre-mir166s. among them, three ests generated from mir166e/q, mir166b/x, mir166o/v were derived from immature seeds, and one est for mir166j from germinating shoot. additionally, five ests generating pre-mir166a, i, n, r, y were originated from the tissues under abiotic stresses. these observations not only confirmed that at least pre-mir166s, out of the soybean mir166s, are indeed transcribed, but also implied that these mir166s might probably be involved in seed development, germination and/or in response to abiotic stresses. targets of mir166s in soybean to obtain a further understanding of biological functions of mir family in soybean, the transcripts of glycine max from jgi genomic project were utilized as a reference system to predict the targets of soybean mir166s by employing the psrnatarget tool. using three representative mature mir166s as query sequences, target genes were identified. as shown in fig. 6a, mir166s were predicted to target ten members of the hd-zipiii family genes, whereas mir166h,k,u,z,l might target the transcripts other than hd-zip iii family genes, such as phytoene synthase, ammonium transporter 2-like and/or two-component response regulator-like aprr this results support the previous report that mir family displays the functional diversity and specificity in plants. among the predicted targets, cleavage of athb14-like transcript was chosen to be validated using rlm-race approach. as shown in fig. 6b-c, the product of athb14-like was precisely cleaved at the 10th or 12th position of complementarity from the 5’ end of mir166s, indicating that mir166s can be indeed involved in the regulation of biological processes via targeting hd-zip iii family in soybean.fig. predicted targets of mir family and validation of mir target gene using rlm-race pcr. a the alignments of fragments of target mrnas that have complementarity to mir166s are shown. mir sequences are displayed to highlight complementarity to target mrnas. the targets are classified into seven groups. b rlm-race pcr for athb14-like. lane and m show cleaved product of desired size and dna ladder, respectively. c mapping of athb14-like mrna cleavage sites by rlm-race. the arrows indicate the cleavage sites and the numbers show the frequency of clones sequenced to understand the evolutionary relationship between mir166s and their target genes, phylogenetic analysis were performed using unique target sequences and target genes from diverse plant species. as shown in fig.  out of utss were clustered with hd-zip family in the phylogenetic tree. it is reasonable since hd-zip iii family is highly conserved targets of mir166s among diverse plant species. further observation pointed that uts and uts were very close to hd-zip iii family genes in soybean including gma-hb14-like, gma-rev-like and gma-hb indicating highly conserved binding site of mi166s in soybean hd-zip iii family. in contrast, two non-hd-zip iii target genes in soybean had seven closely related utss, which were also close to sbi-calcium channel alpha- subunit, ptr-hd-zip iii phv-plastoglobulin-, ptr-atp synthase b chain atpf. similarly, gma-appr was grouped with utss and some target genes including csi-resistance protein-like protein, bna-prk, ath-photosystem i reaction center subunit xi. the observations suggested that non-hd-zip iii target genes are diverse among different plant species.fig. phylogenetic analysis of unique target sequences present in plant species and binding sites of mir in target genes. target genes of mir166s in soybean are highlighted with a star. the details of utss were listed in additional file 10: table s7 to provide clues for their functions during soybean development, the expression patterns of annotated target genes of mir166s were extracted from genome-wide transcript profiling data of soybean tissues from the phytozome database. as shown in fig.  the members of the hd-zipiii family were highly expressed in stem, shoot apical meristem and/or pods, whereas low transcript accumulation was observed in roots and leaves. additionally, the transcript accumulation of aprr psy and amt2-like was highest in flower, root and stem, respectively.fig. expression patterns of mir target genes in different soybean tissues. the transcript profiling data of annotated target genes in soybean tissues were extracted from the publicly-available phytozome database expression pattern of the five newly identified mir genes during seed development or in response to abiotic stresses several studies have documented that mir plays important roles during seed development and in response to abiotic stresses. to provide clues for their functional roles in soybean, we took the five newly identified pre-mir166s as examples to examine their expression patterns during seed development and in response to cold, drought and salinity stresses. expression patterns of the five newly identified mir166s were investigated using seeds collected at 65 days after flowering. as shown in fig. 9a, all the five newly identified mir genes showed gradual increase in transcript accumulation as seeds developed, and reached the peak at daf. their fold changes were up to,,, and, respectively, compared with the ones at daf, suggesting that the five newly identified mir166s might play an important role at late stages of seed development.fig. expression analysis of the five newly identified mir genes in response to seed development and abiotic stresses. a expression patterns of the mir genes during seed development. b expression patterns of the mir genes in seedlings exposed to cold stress for and 72 h. c expression patterns of the mir genes in seedlings exposed to drought stress for and d. d expression patterns of the mir genes in seedlings exposed to salinity stress for and 72 h. error bars indicate se of two biological and three technical replicates. values were normalized against the subi gene ten-day-old soybean seedlings were exposed to cold stress at 4 °c for and h, and expression of the five newly identified mir166s were monitored. as indicated in fig. 9b, the expression of mir166v and mir166z were dramatically increased by - and -fold, respectively, at 3 h after cold treatment, implying that mir166v and mir166z might respond to cold stress at early stage. mir166w showed gradual increase in transcript accumulation as the stress prolonged, while mir166x and mir166y exhibited slight decrease at late stage of cold stress. when soybean seedlings were exposed to drought stress, the expression of mir166v and mir166y was promptly decreased by and folds, respectively, at 2 days after treatment, indicating that they might be relatively sensitive to drought stress. also, a distinct decrease in transcript accumulation after 10 days of drought treatment was observed for mir166w, mir166x, mir166y and mir166z with,, and fold changes, respectively as compared to control. when young seedlings were subjected to salt stress, the expressions of all the mir genes were promptly increased at 3 h after salt treatment. especially, the expression of mir166v and mir166z reached their maximum level at 3 h salt stress with and fold changes, respectively as compared to non-salinity control, suggesting that they might be responsive to salt stress at early stage. mir166w and mir166y showed gradual increase in transcript accumulation as the stress prolonged, and reached their maximum level at 72 h salt stress with and fold changes, respectively as compared to control. expression pattern of mir target genes during seed development or in response to abiotic stresses to explore the regulatory function of mir166s in soybean, the expression patterns of target genes with annotation in phytozome were also analyzed by qrt-pcr during seed development and in response to abiotic stresses. as shown in fig. 10a, the expressions of all the hd-zip iii genes except glyma05g were gradually increased as seeds developed, and reached the peak at daf. their fold changes ranged from to, compared with the ones at daf. on the contrary, psy was gradually decreased in transcript accumulation as seeds developed, while amt2-like was dramatically decreased by folds at daf.fig. expression analysis of the target genes in response to seed development and abiotic stresses. expression patterns of target genes during seed development, under cold stress, drought stress and salinity stress. stress treatments and seed development stage were used as same as the ones in fig.  psy, phytoene synthase; amt2-like, ammonium transporter 2-like. error bars indicate se of two biological and three technical replicates. values were normalized against the subi gene as indicated in fig. 10b, the expressions of all the hd-zip iii genes and aprr were obviously decreased by cold stress, and reached their minimal levels at 24 h stress with – fold changes as compared to control. in contrast, the transcript accumulations of psy and amt2-like were increased by and folds at and 48 h after treatment, respectively. four athb14-like, rev-like and psy displayed gradual decrease in transcript accumulation as drought stress prolonged, and their maximal fold changes ranged from to. although atm2-like displayed similar expression level to the control at early stage, the transcript accumulation was dramatically increased from 6 days after treatment until 10 days with fold changes, suggesting that atm2-like might respond to drought stress at middle and late stages. the decreased transcript accumulation was also observed for the rest of target genes on one or more treatment points as compared to control. as shown in fig. 10d, all the hd-zip iii genes and aprr were promptly increased to their peaks in transcript accumulation at 3 h after salt stress with – fold changes, respectively as compared to control, indicating that they might be responsive to salt stress at early stage. the expressions of psy and amt2-like were gradually increased, and reached their peaks at 72 h after salt treatment with – fold changes, respectively as compared to control. discussion mir acts as a highly conserved mirna to be implicated in a wide range of cellular and physiological processes in plants. based on the data from mirbase, most of plant species harbor multiple mir166s. it has been accepted that the existence of multiple copies of mir genes within plant species provides one possibility for functional redundancy and specificity. thus, the diversity and evolutionary conservation of mir family genes are expected to be elucidated for further understanding their fine-tuned regulatory roles in plant developmental processes and/or resistance to stresses. generally, mirnas are annotated using comparative criteria for both their expression and biogenesis, including size and sequence of mature mirna, phylogenetic conservation, secondary structure of mirna precursor and increased precursor accumulation with decreased dicer function in vivo. in this study, the five newly identified mirnas exhibited the same size and sequence as their paralogs in soybean. the lengths of the potential mirna precursors range from to 244 nt, which is acceptable since the length of plant mirna precursors has been found to be in the range of 55–930 nt. furthermore, their precursors also satisfied the criteria with aspect to both phylogenetic conservation and secondary structure. it has been proposed that low mfe and high mfei are reliable criteria to distinguish mirnas from all coding and other non-coding rnas, as these parameters, to certain extent, reflect the stability of the prefect or near-perfect secondary hairpin structure of mirna precursor. in the study, mfe values of the five newly identified pre-mir166s ranged from - to -  kcal/mol, which were much lower than folding free energies of trna and rrna. at the same time, their individual mfei value was also higher than that of trna, rrna and mrna. taken together, the five newly identified mirnas in soybean are very likely to be true mirnas. thus, the study enlarged the size of soybean mir family from to members. it has been proposed that conserved mirna families have been both conserved and diversified during mirna evolution. in this study, our data indicate that soybean mir166s showed evolutionary conservation. as shown in fig. 2a, ten pre-mir paralogous pairs with high sequence identity were individually grouped into a discrete clade in the phylogenetic tree. further analysis revealed that eight out of ten paralogous pairs were located at corresponding loci on the same duplication block, indicating that they might have arisen from segmental duplications without major sequence change. additionally, smaller than of ka/ks values suggest that these mir genes likely had undergone through a purifying selection with limited functional divergence after duplication. it is also worth noting that five paralogous pairs were clustered into discrete clades in phylogenetic tree based on their promoter sequences, supporting the view that gene duplication could happen simultaneously at the coding region and the promoter. evidently, the origins of these soybean mir166s were in line with soybean evolutionary history that soybean had undergone whole-genome duplication event ~ million years ago. apart from the evolutionary conservation, soybean mir166s also exhibited evolutionary diversification. in the study, five pre-mir166s were separated far from the majority of pre-mir166s, and the identity between the five pre-mir166s and the rest of the pre-mir166s was low. these observations suggested that they had diverged from the majority of mir166s during soybean evolution. mature mir166s are ultimately responsible for functional diversification since they act as functional unit to regulate the expression of their targets. in the study, all the mir166s except mir166l were highly conserved in mature sequence, and identical mature sequence exists in each of the pre-mir paralogous pairs with exception of pre-mir166h and pre-mir166j. consequently, most of soybean mir166s shared the same putative targets in soybean. however, mir166z, k, h, u, l are likely to have functional diversity since they are predicted to target genes other than those of the hd-zip iii family, i.e., aprr phytoene synthase and ammonium transporter 2-like. thus, our findings demonstrate that mature mir166s exhibit both functional conservation and diversification. it has been well-known that mir plays important roles during seed development. for example, mir can repress the seed maturation program during vegetative development via controlling the expression of hd-zip iii family genes such as phb; while it might also be involved in somatic embryo formation via acting as a positional signal from the basal–peripheral region of early embryos. in this study, all the five newly identified mir166s showed gradual increase during seed development, and reached their expression peaks at daf, a late stage of soybean seed development. these expression patterns are consistent with previous report that mir in brassica napus was strongly accumulated at late stage of seed development. also, it was proposed that the accumulation of mature mir reached its major peak at late cotyledonary embryo stage in larch, suggesting that it might be associated with cotyledon formation. the expression patterns of the five newly identified mir166s suggested that they might be one of the major contributors to the network controlling seed development, especially seed maturation program. it has been documented that plant mir166s might be involved in response to diverse abiotic stresses, and main evidences were derived from the altered accumulation pattern of mature mir under diverse abiotic stresses. for instance, the accumulation of mature mir was decreased by drought stress in sorghum bicolor, whereas drought stress led to increase of mirn accumulation in root and leaf of wheat. another kind of examples are that mir in potato was significantly up-regulated by salinity stress, and depression of mir by small tandem target mimic can lead to improvement of salt tolerance in arabidopsis. in the study, the presence of stress-responsive elements in the promoter regions of mir genes pointed to their possible roles under various stresses. furthermore, the altered expression patterns of the five newly identified mir genes suggest their prominent roles in response to cold, drought and salinity. for example, mir166v and mir166w harbored the ltr cis-acting element in promoter regions, which can be involved in low-temperature responsiveness. indeed, mir166v was promptly increased by cold stress with fold changes, while mir166w showed gradual increase in transcript accumulation as cold stress prolonged, implying that mir166v and mir166w can possibly function as a modulator of cold-induced signaling pathways. these results are consistent with previous reports that mir166e and mir166r were increased by chilling stress in vegetable soybean, and that the expression levels of mir166u were obviously increased under cold condition in nitrogen-fixing nodules of soybean. it has been proposed that spatial-temporal or specific expression of mirna isoforms is a crucial determinant of isoform specificity. in the study, the differential expression patterns of mir genes under same stress indicated their functional specificity. as shown in fig. 7c, mir166v was dramatically decreased at 2 days after drought treatment, whereas mir166w was increased at 8 days after drought treatment. thus, the stress-specific properties of mir genes provided important clues for elucidating their functional specificity. although five novel mir166s in soybean were differentially regulated during seed development or in response to stresses, their expression patterns were not negatively correlated with the ones of mir target genes. for example, the expressions of the hd-zip iii genes were gradually increased as seeds developed. this is consistent with the previous report that hd-zip iii family in arabidopsis can be positively involved in regulation of seed mature program. however, the five new identified mir166s also showed gradual increase in transcript accumulations. these results suggested that mir166-mediated silencing of target genes was under sophisticated regulation. mir gene family has large number of members in plant species, and at least mir members exist in soybean. so, we speculated that mir gene family possibly contributes to the accumulation pattern of target genes in a family member-specific manner. it has been widely accepted that members of a certain evolutionary branch have potential to share similar functions. in the study, pre-mir166x and pre-mir166y were the most closely related paralogs in the phylogenetic tree, and they also showed % sequence identity. as expected, similar expression patterns were observed for mir166x and mir166y in response to cold, drought and salinity stresses. additionally, the most closely related pre-mir paralogs also shared similar tissue-specific expression pattern in soybean. thus, it was hypothesized that these most closely related pre-mir paralogs in soybean perform similar cellular functions during developmental processes or in response to various stresses.  this study enlarged the size of soybean mir family from to members. a comprehensive analysis indicated that the soybean mir166s exhibited evolutionary conservation and diversification. furthermore, the expression patterns of five newly identified mir166s during seed development and in response to abiotic stresses provided important clues for elucidating the functional specificity and redundancy of mir family in soybean. future study will aim at investigating the effect of each mir gene on seed development and stress tolerance, identifying target genes and dissecting the regulatory mechanisms for each mir and exploring combination of multiple mir genes during seed development and in response to specific stress. 
phylogenetic analyses of the mir166s and their precursors indicated that soybean mir family exhibited both evolutionary conservation and diversification, and ten pairs of mir precursors with high sequence identity were individually grouped into a discrete clade in the phylogenetic tree. the targets of soybean mir166s were predicted, and the cleavage of athb14-like transcript was experimentally validated by race pcr. this study enlarged the size of soybean mir family from to members, and the soybean mir166s exhibited evolutionary conservation and diversification. further, the expression patterns of the five newly identified mir166s and target genes were examined during seed development and in response to abiotic stresses, which provided important clues for dissecting their functions and isoform specificity. these findings have laid a foundation for elucidating functional conservation and diversification of mir family members, especially during seed development or under abiotic stresses. microrna is a highly conserved family of mirnas implicated in a wide range of cellular and physiological processes in plants. the analysis of genomic organization and evolution of mir gene family revealed that eight segmental duplications and four tandem duplications might occur during evolution of the mir family in soybean. the cis-elements in promoters of mir family genes and their putative targets pointed to their possible contributions to the functional conservation and diversification. mir family generally comprises multiple mir members in plants, which might exhibit functional redundancy and specificity. however, the evolutionary conservation and functional diversification of mir family members in soybean remain poorly understood. mir familysoybeanevolutionary conservation and diversificationpromoter analysisgene expression patternhttp://dx.doi.org/ /501100001809national natural science foundation of china31370340bian shaomin issue-copyright-statement© the author 2017 the soybean mir family consists of members according to the mirbase database. we identified five novel mir166s in soybean by data mining approach, thus enlarging the size of mir family from to members.
 microrna is a highly conserved family of mirnas implicated in a wide range of cellular and physiological processes in plants. mir family generally comprises multiple mir members in plants, which might exhibit functional redundancy and specificity. the soybean mir family consists of members according to the mirbase database. however, the evolutionary conservation and functional diversification of mir family members in soybean remain poorly understood.  we identified five novel mir166s in soybean by data mining approach, thus enlarging the size of mir family from to members. phylogenetic analyses of the mir166s and their precursors indicated that soybean mir family exhibited both evolutionary conservation and diversification, and ten pairs of mir precursors with high sequence identity were individually grouped into a discrete clade in the phylogenetic tree. the analysis of genomic organization and evolution of mir gene family revealed that eight segmental duplications and four tandem duplications might occur during evolution of the mir family in soybean. the cis-elements in promoters of mir family genes and their putative targets pointed to their possible contributions to the functional conservation and diversification. the targets of soybean mir166s were predicted, and the cleavage of athb14-like transcript was experimentally validated by race pcr. further, the expression patterns of the five newly identified mir166s and target genes were examined during seed development and in response to abiotic stresses, which provided important clues for dissecting their functions and isoform specificity.  this study enlarged the size of soybean mir family from to members, and the soybean mir166s exhibited evolutionary conservation and diversification. these findings have laid a foundation for elucidating functional conservation and diversification of mir family members, especially during seed development or under abiotic stresses. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords mir familysoybeanevolutionary conservation and diversificationpromoter analysisgene expression patternhttp://dx.doi.org/ /501100001809national natural science foundation of china31370340bian shaomin issue-copyright-statement© the author 2017 