Genome-Wide Identification of Different Dormant Medicago sativa L. MicroRNAs in Response to Fall Dormancy

Background MicroRNAs (miRNAs) are a class of regulatory small RNAs (sRNAs) that regulate gene post-transcriptional expression in plants and animals. High-throughput sequencing technology is capable of identifying small RNAs in plant species. Alfalfa (Medicago sativa L.) is one of the most widely cultivated perennial forage legumes worldwide, and fall dormancy is an adaptive characteristic related to the biomass production and winter survival in alfalfa. Here, we applied high-throughput sRNA sequencing to identify some miRNAs that were responsive to fall dormancy in standard variety (Maverick and CUF101) of alfalfa. Results Four sRNA libraries were generated and sequenced from alfalfa leaves in two typical varieties at distinct seasons. Through integrative analysis, we identified 51 novel miRNA candidates of 206 families. Additionally, we identified 28 miRNAs associated with fall dormancy in standard variety (Maverick and CUF101), including 20 known miRNAs and eight novel miRNAs. Both high-throughput sequencing and RT-qPCR confirmed that eight known miRNA members were up-regulated and six known miRNA members were down-regulated in response to fall dormancy in standard variety (Maverick and CUF101). Among the 51 novel miRNA candidates, five miRNAs were up-regulated and three miRNAs were down-regulated in response to fall dormancy in standard variety (Maverick and CUF101), and five of them were confirmed by Northern blot analysis. Conclusion We identified 20 known miRNAs and eight new miRNA candidates that were responsive to fall dormancy in standard variety (Maverick and CUF101) by high-throughput sequencing of small RNAs from Medicago sativa. Our data provide a useful resource for investigating miRNA-mediated regulatory mechanisms of fall dormancy in alfalfa, and these findings are important for our understanding of the roles played by miRNAs in the response of plants to abiotic stress in general and fall dormancy in alfalfa.


Introduction
There are two abundant classes of small RNAs (sRNAs) in plants; i.e., endogenous small interfering RNAs (siRNAs) and microRNAs (miRNAs) [1,2]. Both siRNAs and miRNAs regulate messenger RNA (mRNA) stability and translation [1], and the miRNAs regulate gene expression by targeting specific mRNAs [2]. MiRNAs were firstly discovered as developmental timing regulators in Caenorhabditis elegans [3]. There is evidence that miRNAs play a critical role in posttranscriptional gene regulation and various biological and metabolic processes [4]. For instance, miRNAs regulate the modulation of the processes associated with growth and development in plants, including leaf morphogenesis, floral organs and root development [5][6][7]. Many new conserved Hg-responsive miRNAs and miRNA candidates were identified from Medicago truncatula [8]. Some of the miRNAs were involved in the regulation of development and plant responses to heavy metal stress; they were up-or down-regulated by heavy metals (Hg, Cd, and Al) [9]. In addition, miRNAs are involved in responses to various abiotic and biotic stresses, such as drought, cold, salinity and nutrient starvation [10][11][12].
Leguminous plants account for one third of primary crop production in the world [13][14][15]. In legumes, alfalfa (Medicago sativa L.) has been gradually concerned as an important forage crop in the world, due to its high nutrition value, wide adaptation, yield potential [16] and as a kind of nutritious feed. However, alfalfa varieties developed from different geographic regions have distinct characteristics in autumn. Ample studies have shown that there were high correlations between fall growth and winter injury in alfalfa [17][18][19][20][21], and the sugar content of roots and crown buds in alfalfa were closely related with winter hardiness [22][23][24][25][26]. Recent studies also revealed that winter injury was negatively associated with carbohydrate content in the roots of alfalfa, especially the hexosan content [27]. Increased sucrose, stachyose and raffinose, and reduced glucose, fructose and starch levels were related to freezing tolerance [24]. These studies mainly concerned the effects of temperature on fall dormancy of alfalfa. However, it has been reported that photoperiod was a key factor regulating fall dormancy in alfalfa [28]. Based on the responses of different varieties to low temperatures and short day-length (photoperiod), alfalfa varieties were classified into three types; i.e., fall dormant varieties grow very slowly and even cease to grow with short and decumbent shoots in the autumn. In contrast, non-dormant varieties continue to grow and produce tall upright shoots in the autumn. Semi-dormant varieties show phenotypes between dormant and non-dormant alfalfa [16,29]. In many countries, such as the United States and Canada, FD class is commonly used as the first index of selecting alfalfa varieties because of its important role in various adaptations to particular regions associated with winter survival [30].
Recently, high-throughput sequencing analysis of sRNAs from Medicago truncatula leaves identified eight new miRNAs and sequencing sRNAs from specific tissues of legume plants by deep sequencing is expected to reveal more new miRNAs [31]. To identify the miRNAs in response to Hg, two small RNA libraries were generated from Hg-treated and Hg-free (control) Medicago truncatula seedlings. 201 individual miRNAs were identified representing 63 known Medicago truncatula miRNA families, including 12 new conserved and one non-conserved miRNAs that have not been described previously [8]. Zhou et al. used a bioinformatics approach for ESTs (Expressed Sequence Tags)-and GSS (Genomic Survey Sequences)-wide prediction of novel miRNAs in Medicago truncatula [9]. Lelandais et al. showed that spatial regulation of miRNAs may determine the specialization of regulatory RNA networks in plant differentiation processes in Medicago truncatula, such as root nodule formation [32].
Though the genome sequence of Medicago truncatula has been completed, which is an annual legume species with a small diploid genome and easy transformation and is a reference model plant to study functional genomics of alfalfa, the physiological, biochemical and molecular mechanisms causing FD are not clear. The genome sequence of alfalfa, which is an autotetraploid, allogamous and heterozygous species, is not available at present. Currently, we have sequenced the transcriptome of alfalfa using RNA-seq technology and identified some differentially expressed genes. Regulation of gene expression can be achieved at transcriptional, post-transcriptional and translational levels. In plants and animals, gene silencing occurs when endonuclease complexes are guided by small RNAs to target RNAs. For instance, a more effective approach to induce gene silencing in the alfalfa has been reported. Artificial microRNAs (amiRNAs) can be used to very specifically target genes for silencing because only a short sequence of 21 nucleotides of the gene of interest is used [33]. It is important to identify small RNAs and their target messenger RNAs, which are essential to understanding small RNA-mediated gene regulation of growth and stress responses and fall dormancy.
In comparison with miRNAs from Medicago truncatula, fewer miRNAs in Medicago sativa have been identified. In this study, we used previously known Medicago truncatula miRNAs to search for some miRNAs in response to fall dormancy of Medicago sativa. In this context, high-throughput sequencing has been used to identify alfalfa target miRNAs (including fall dormant and non-fall dormant types at two time points). To understand the role of miRNAs in response to fall dormancy of alfalfa, the objectives of this study were: (i) to sequence and annotate miRNAs of alfalfa varieties in different fall dormant types, and their expression in response to fall dormancy in standard variety (Maverick and CUF101), and (ii) to identify the target miRNAs related to fall dormancy in alfalfa and provide a scientific basis for miRNA function in fall dormancy of two standard alfalfa varieties (Maverick and CUF101).

Plant materials and growth condition
There are 11 fall dormancy classes (FDC) for Medicago sativa L. (alfalfa) at present, which are classified by the regrowth height after cutting [34] and can be broadly divided into dormant (FDC 1-4), semi-dormant (FDC 5-7) and nondormant categories (FDC 8-11) [29]. In this study, alfalfa standard varieties Maverick (FDC1) and CUF101 (FDC9) were obtained from the United States and planted on a sandy loam soil at the Experimental Station of Henan Agricultural University, Zhengzhou, China (34˚196N, 113˚356E). The region experiences a monsoon-influenced, four-season humid subtropical climate with 640.8 mm of annual precipitation. Over the past 62 years, it has had a mean temperature of 14.7˚C with extremes of 217.9˚C and 43˚C. Spring lasts from February 3 to May 4, summer from May 5 to August 6, autumn from August 7 to November 6, and winter from November 7 to February 2. Soil samples were collected and tested for amounts of nitrogen (N), phosphorous (P), and potassium (K). According to the results, the land was prepared by fertilizing 81 kg/ha of N as urea and 96 kg/ha with P 2 O 5 as calcium phosphate to a depth of 15 cm. The field was in an alfalfaalfalfa-alfalfa-corn (Zea mays)-winter wheat (Triticum aestivum L.) cropping system before the current study. On October 1, 2009, seed was sprinkled in 5 rows of alfalfa with 0.6 m row spacing plots. During drought, irrigation was provided in a timely manner. Weed control was performed by hand or hoeing, and insects were controlled as required.
Alfalfa leaves were collected from Maverick and CUF101, both on May 19 and September 23, 2011, which was 14 days after cutting. Samples of two different alfalfa varieties in May and September were named DM1 (Maverick, dormancy-May-1), NM2 (CUF101, non-dormancy-May-2), DS3 (Maverick, dormancy-September -3), NS4 (CUF101, non-dormancy-September-4), respectively. The two time points were selected because both the fall-dormant (Maverick) and nonfall dormant (CUF101) alfalfa grows normally in May, but in September, Maverick enters into fall dormancy while CUF101 still continues growing in our experimental station. Leaves were immediately frozen in liquid nitrogen and then stored at 280˚C.

RNA isolation and qualification
Total RNA was isolated using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. RNA degradation and contamination was monitored on 1% agarose gels. The RNA samples were measured with an ND 1000 spectrophotometer (Nanodrop) for contamination with either protein (A260 nm/A280 nm ratio) or reagent (A260 nm/A230 nm ratio). All the sample integrity numbers (RINs) were .7.5.

Library construction and deep sequencing
Total RNAs were extracted from DM1, NM2, DS3 and NS4 following the manufacturer's protocol. A total amount of 3 mg total RNA (>100 ng/ml) per sample was used as input material to construct libraries. Sequencing libraries were generated using a NEB Next Multiplex Small RNA Library Prep Set for Illumina (NEB). Following manufacturer's recommendations, index codes were added to attribute sequences to each sample. Briefly, NEB 39 SR Adaptor was directly and specifically ligated to the 39 end of miRNA. After the 39 ligation reaction, we hybridized the SR RT primer. This step was important to prevent adaptor-dimer formation. For each sample, small RNAs (18-30 nt) were ligated with 59 and 39 RNA adapter by T4 RNA ligase, after they were purified by electrophoretic separation on a 15% TBE-urea denaturing PAGE gel. The sRNAs were transcribed into cDNA using Superscript II reverse transcriptase (Invitrogen, USA). Products were purified and quantified using the Agilent high sensitivity DNA assay on the Agilent Bioanalyzer 2100 system. After the clustering of the index-coded samples, four libraries, three for non-dormancy samples, one for dormancy sample, were sequenced on an Illumina Hiseq 2000 platform at the Beijing Novo Gene Genomics Institute, China.

Analysis of sRNA sequencing data
Raw sequence reads were first processed through data cleaning by removing reads containing ploy-N, with 59 primer contaminants, without 39 primers or the insert tag, containing ploy A or T or G or C and low quality reads. At the same time, Q20, Q30, and GC contents of the raw data were calculated. Then, we chose a certain range of length from the clean reads to do all the downstream analyses. The available software miREvo and mirdeep2 were integrated to predict novel miRNA by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the small RNA tags unannotated in the former steps [35,36].

Mapping of sequence reads and expression analysis of miRNAs
Mapping of sequence reads back to the transcriptome was carried out using the alignread.pl script that comes with the Trinity assembler with the aligner bowtie (--bowtie) [37]. Sequences mapped to the genome using the Medicago truncatula (Mt3.0) used SOAP [38]. The expression levels of miRNA were estimated by transcript per million (TPM) through the following criteria Normalization formula: Normalized expression5mapped read count/Total reads610 6 [39]. Pvalues were adjusted using Q-value. Q-value,0.01 and |log 2 (fold change)|.1 was set as the threshold for significantly differential expression by default [40].

Prediction and annotation of candidate potential miRNAs
Target accessibility allowed maximum energy to unpair the target-site (UPE). Target predictions were performed using the psRNA Target http://bioinfo3.noble. org/psRNA-Target/ and psRobot_v1.2 [41], through alignment with the genome of Medicago truncatula and the transcriptome of Medicago sativa, and then we took the intersection of the two predicting tools at the genome level of Medicago truncatula and transcriptome of Medicago sativa.

Expression validation of different fall dormancy alfalfa miRNAs
This validation was done on five individual RNAs. We used SYBR Prime Script miRNA RT-PCR (TaKaRa, Code No.: RR716) for the RT reactions. All primers are listed in Table S1. The kit could add ploy (A) to the 39 end of miRNAs and start reversely transcription. With known sequence at its 59 end, the reverse transcription was led by special oligo-dT ligations. Real-time quantitative PCR was conducted on a light cycler Real-Time PCR Detection System. The reactions were incubated in a 96-well plate at 95˚C for 300 s, followed by 40 cycles of 95˚C for 10 s, 60˚C for 10 s, and 72˚C for 10 s. Each reaction included 1 mL of product from the diluted RT reactions, 0.8 mL of each forward primer, 0.8 mL of Uni-miR qPCR Primer, 10 mL of SYBR Premix Ex TaqII (TaKaRa, Code No.: RR716), and 7.4 mL of nuclease-free water. All reactions were run in three replicates for each sample. A relative quantitative method (2 2DDCt ) was used to evaluate quantitative variation [42].

Northern blot analysis
Northern blotting was performed as previously described [43][44][45]. Total RNA was isolated using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. About 30 mg of total RNA samples were run on 15% polyacrylamide denaturing (urea) gels, and then transferred to a Hybond-N+ nylon membrane (Amersham Biosciences) by electrophoresis using a semidry transfer cell (Bio-Rad). Hybridization was performed according to a standard protocol. Digoxin-labeled oligo-nucleotide probes complementary to the mature miR-novel_1, miR-novel_3, miR-novel_45, miR-novel_17 and miR-novel_76 were used in the hybridization. The probes are listed in Table 1.

Natural height and leaf area of alfalfa
Leaf areas were measured using A5KLW (A: leave area; K: coefficient of correction; L: Leave length; W: leave wide [46]). Twenty randomly selected plants were measured at each of the four conditions (alfalfa has a terrately compound leaf). Natural height (S Zhang, personal communication) and leaf areas are shown in Table 2. In Figure S1, DM1 (Maverick, May) and NM2 (CUF101, May) are shown to be growing well and similarly in May, while DS3 (Maverick, September) turned into a decumbent state and became much shorter than DM4 (CUF101, September) in September. From the plant height, Maverick might enter fall dormancy in September, while CUF101 did not show fall dormancy morphology at this time and continued to grow and produce tall upright shoots in the fall with a similar natural height and leaf area as in May and September [see Materials and methods].

Analysis of sequences from four libraries
To investigate the enrichment of miRNAs and identify the miRNAs in response to fall dormancy, we generated four libraries from alfalfa leaves sampled under four conditions. After sequencing via Illumina Hiseq 2000, we removed low-quality reads and corrupted adapter sequences (reads,18 nt). High-throughput sequencing generated 7,536,934 reads (1,693,394 unique) representing DM1, 6,572,924 reads (1,985,915 unique) for NM2, 6,548,094 reads (1,449,168 unique) for DS3, and 6,627,260 reads (1,185,319 unique) for NS4. In order to control the quality of the libraries, we inspected the error rate and GC content of the sequencing results. They are presented in Table 3. The error rate distribution and GC distribution were also configured (ref. Figures S2 and S3).
The expression of miRNAs in dormant alfalfa (DM1 and DS3) and nondormant alfalfa (NM2 and NS4) were examined with Solexa technology. The length distribution of reads showed that the majority of the reads, except for corrupted adapter sequences, were 20-25 nt in size ( Figure 1). In the four libraries, 21 nt and 24 nt RNAs were the four most abundant classes ( Figure 1). This was consistent with the distribution patterns of small RNAs from other plant species [32]. After initial processing, the high quality small RNA reads were mapped to the alfalfa transcriptome sequence (NCBI accession number: SRA057663) (S Zhang, personal communication). The number of total/unique sequences that matched

Identification of known miRNAs and novel miRNAs in alfalfa
To identify known miRNAs in dormant and non-dormant alfalfa, we aligned the sRNA sequences with known mature miRNAs from plants in the miRBase20.0. Modified software mirdeep2 [36] and srna-tools-cli were used to obtain the potential miRNAs and draw the secondary structures. We identified 866, 923, 777 and 797 unique sequences as known miRNAs in the four alfalfa libraries (Table 4), respectively. Novel miRNAs were predicted according to the characteristic hairpin structures of their precursors, which distinguishes them from other endogenous sRNAs [47,48]. A total of 51 novel hairpin miRNAs candidates (206 families) were identified (Table S2, Table S3), which were not found in other species. These miRNA candidates are likely to be new miRNAs or new members of known miRNA families in alfalfa.

Potential miRNAs related to fall dormancy in alfalfa
We took the intersection of the two predicting tools (psRNA Target and psRobot_v1.2) at the genome level of Medicago truncatula (Table S4) and the transcriptome of Medicago sativa (Table S5). We identified 566 potential miRNAs related to fall dormancy at the genome level of Medicago truncatula, and 583 potential miRNAs related to fall dormancy at the level of the transcriptome of Medicago sativa. Their predicted function annotation, target accession, target start, target end, expectation, UPE, target aligned fragment are listed in Table S4 at the genome level of Medicago truncatula and an Table S5 at the transcriptome level of Medicago sativa.

Functional annotation of differentially expressed miRNAs
Gene Ontology (GO) enrichment analysis was used on the target gene candidates of differentially expressed miRNAs [49]. GO functional enrichment analysis was performed to characterize the functional consequences of changes in gene expression associated with fall dormancy. According to the GO enrichment, the differential expression miRNA target genes in alfalfa related to defense responses, signal transduction, signaling, single organism signaling, cell communication, cellular responses to a stimulus, regulation of cellular processes, response to stress, biological regulation, regulation of biological processes, ADP binding, hydrolase activity, and pyrophosphatase activity (Figure 2). From Figure 2, we can see that fall dormancy is mainly a biological process in alfalfa. The main GO enriched factors are biological regulation, regulation of biological processes, response to stress, and regulation of cellular processes.
KEGG pathways were used to assess the statistical enrichment of the target gene candidates via KOBAS software. In three sets, the statistical enrichment of the gene candidates in KEGG pathways showed that the majority of target gene clusters included metabolic pathways, biosynthesis of secondary metabolites, porphyrin and chlorophyll metabolism, glycerolipid metabolism, cell cycle, glycerolipid metabolism, meiosis, plant hormone signal transduction, glycerophospholipid metabolism, ubiquitin-mediated proteolysis, etc. (Table S6). From the most enriched pathway terms, we can see that the main enriched factors included metabolic pathways, biosynthesis of secondary metabolites, cell cycle, and ubiquitin mediated proteolysis. From the analysis of GO enrichment and KEGG pathways, we inferred that the regulation mechanism of fall dormancy was biological regulation of responses to stress in the metabolic pathway.

Response of known miRNAs to fall dormancy in different types of dormant alfalfa
To identify fall dormancy-responsive miRNAs in alfalfa, the normalized expression of miRNAs in the DS3 library was compared to that in the DM1, NM2 and NS4 libraries (Tables S7, 8, and 9). In contrast to DS3 vs DM1, DS3 vs NM2 and DS3 vs NS4, there were 11 up-regulated miRNAs and nine down-regulated miRNAs found in three sets ( Figure S4). The differences in the expression profiles of these miRNAs indicated that each miRNA may play a different role in alfalfa, and they may explain the differences between dormant alfalfa and non-dormant alfalfa. Based on the results of high-throughput sequencing, we selected 14 miRNAs from these 20 known miRNAs randomly with changes in expression levels with Log 2 (Fold Change).51 (P value,0.01) [40] in response to fall dormancy treatment to validate the expression patterns with real-time quantitative PCR. As shown in Figure 3, high-throughput sequencing and RT-PCR produced similar results ( Figure 3). Figure S3 shows miR172a, miR52671, miR2676a, miR2593e, miR5286b, miR5281a, miR2634, miR5286a, and miR2629h to be down-regulated in different types of dormant alfalfa. Conversely, miR5228, miR5205b, miR5279, miR2674, miR2655b, miR164a, miR164d, miR156g-3p, miR156e, miR2643a, and miR5299 were up-regulated in different types of dormant alfalfa. We picked up these common miRNAs from three contrasts; i.e., DS3 vs DM1, DS3 vs NM2 and DS3 vs NS4 ( Figure S2), and gave their expression pattern and function annotations in Table 5. Because little is known about miRNA-mediated developmental regulation of fall dormancy in alfalfa, we can see that most function annotations of known miRNAs in response to fall dormancy in different types of dormant alfalfa are hypothetical proteins, as shown in Table 5.

Response of new miRNAs to fall dormancy in different types of dormant alfalfa
Among novel miRNAs, there were five candidate up-regulated miRNAs and three down-regulated miRNAs found in three sets (Table S10). It shows novel_17, novel_75, and novel_76 to be down-regulated in different dormant varieties. Conversely, novel_1, novel_3, novel_42, novel_45, and novel_85 were upregulated in different types of dormant alfalfa. Based on the results of highthroughput sequencing ( Figure 4A), we selected five novel miRNAs from these eight novel miRNA candidates randomly. The five novel miRNA levels were further confirmed by Northern blotting analysis ( Figure 4B). As shown in Figure 4, high-throughput sequencing and Northern blotting analysis of the five Figure 2. GO functional enrichment of genes differentially expressed in alfalfa. Gene classification based on gene ontology (GO) for the differential expression miRNA target genes in DS3 vs DM1, DS3 vs NM2 DS3, and DS3 vs NS4. The number of genes in GO terms was analyzed using GO Slim Assignment. Biological processes and molecular function were used for GO analysis. (Hydrolase activity*5hydrolase activity, acting on acid anhydrides, in phosphorus-containing anhydrides;hydrolase activity**5hydrolase activity, acting on acid anhydrides; BP, Biological Process; MF, Molecular Function). doi:10.1371/journal.pone.0114612.g002 novel miRNAs produced similar results ( Figure 4A, and 4B). The increasing miR-novel_1, miR-novel_3, and miR-novel_42 levels were confirmed by Northern blotting analysis. At the same time, the decreasing miR-novel_17 and miR- novel_76 levels also had a similar variation trend with sequencing. The miRNA cluster analysis of differentially expressed sRNAs is shown in Figure S5. It was performed based on fold-changes between DM1, NM2, DS3, and NS4. Cluster analysis of differentially expressed sRNA in different types of dormant alfalfa was displayed using different colors. One hundred and forty six miRNAs (known  miRNAs and novel miRNAs candidates) were potentially involved in the regulation of fall dormancy, as shown in Figure S5.

Discussion
For most eukaryotic cells, miRNAs are a class of regulatory sRNAs involved in gene regulation. Recent studies have shown that miRNAs play an important role in diverse abiotic and biotic stresses, including drought, cold, nutrient starvation, oxidative stress, submergence, UV-B radiation, and viruses [5-7, 10, 39, 50, 51]. Medicago truncatula, which is an annual legume species and a model plant, has decumbent shoots and the yield is very low. In recent years, knowledge of the function of miRNAs in plants has greatly advanced; high-throughput sequencing technology was used to identify Medicago truncatula miRNAs. Alfalfa (Medicago sativa L.), which is autotetraploid, allogamous and a perennial plant, has been considered as one of the most often grown perennial forage crops worldwide. However, alfalfa varieties developed from different geographic regions have distinct characteristics in FD, and few studies have been performed with the express purpose of identifying and analyzing miRNAs in response to fall dormancy of alfalfa. In alfalfa, fall dormancy is an adaptive characteristic related to biomass production and winter survival [52], but the physiological, biochemical, and molecular mechanisms causing FD and how this process occurs is still not clear at present. Dormant alfalfa varieties reduce shoot growth, even cease growth and produce short and decumbent shoots after cutting in the autumn, while nondormant varieties continue normal shoot regrowth and produce tall upright shoots in the autumn. Dormant alfalfa shows fall dormancy morphology in late summer or early autumn, and the shortening day length and falling temperatures in the seasons have been accepted as the environmental factors that induce dormancy in alfalfa [53]. Plant height is the main index indicating fall dormancy in alfalfa [54], and leaf area could also indicate these variations. From the plant height and leaf areas, we can see that Maverick (DM1) and CUF101 (NM2) grow in normal situations in May, Maverick (DS3) might enter fall dormancy in September, while CUF101 (NS4) did not show fall dormancy morphology at this time.
Medicago truncatula has decumbent shoots all year round, only for grazing pasture. There are many papers on genome-wide identification of microRNAs from different treatments of Medicago truncatula. According to the newest miRNA database (miRBase 20, released in June 2013), Medicago truncatula has the greatest amount of precursors and mature miRNAs in plants [6,55,56], even more than Arabidopsis and rice [57]. However, the functional study level of miRNAs in Medicago truncatula is not as high as Arabidopsis and rice, only those of highly conserved miRNAs, such as miR156, miR164 and miR172, have been investigated in model plants. Several conserved miRNAs were identified in the genome of Medicago truncatula, based on homology to miRNAs present in other plant species, and their expression was confirmed in several organs/tissues. Five of ten conserved miRNAs analyzed here showed differential expression in water deficits as compared to well-watered plants: miR166 was slightly up-regulated in the roots; miR169 was down-regulated in the roots and miR398a/b and miR408 were up-regulated in both shoots and roots. Conversely, miR156, miR160, miR171, miR319 and miR393 expression levels were not affected by water deficits [10].
In Arabidopsis thaliana, the predicted targets of miR156 regulated SPL transcription factors, defining an endogenous flowering pathway [58]. Likewise, miR172 regulates AP2-like genes that are involved in controlling floral organ identity in barley, maize, and rice [6,55,56]. In addition, miR156 regulates the expression of miR172 via SPL9, which directly promotes the transcription of miR172b [59]. It has been shown that decreases in N-acetylcysteine 1 (NAC1) mRNA levels due to inducible expression of miR164, result in reduction in lateral root emergence in Arabidopsis thaliana [60]. In Medicago truncatula, both highthroughput sequencing and RT-qPCR showed that miR169 was down-regulated, and the suppression of miR164 expression may contribute to increases in root/ shoot ratios under drought stress [61]. The roots of Medicago truncatula have the ability to interact with rhizobia to develop nitrogen-fixing root nodules. MiR164 may be the first miRNA found to be involved in nodule formation, which can target MtHAP2-1 encoding a transcription factor of the CCAAT-binding family with an ability to control nodule meristem function. Overexpression of miR164 leads to blockage of nodule development by down-regulating the expression of MtHAP2-1. MiR164 also has the ability to control the transcript levels, as well as the expression patterns, of their targets, suggesting that they might contribute to developmental robustness [48,62,63].
Exploring the expression profile of some key miRNAs may be very important for understanding the mechanism of fall dormancy in alfalfa. We compared three sets of different miRNA expressions; i.e., DS3 vs DM1, DS3 vs NM2 and DS3 vs NS4, and each set had 11 up-regulated miRNAs and nine down-regulated miRNAs. Among these miRNAs, conserved miRNAs miR156 and miR164 were up-regulated target genes and miR172 was a down-regulated target gene. It is conceivable that the different expressions of miR172, miR164 and miR156 miRNA families regulate floral repression and influence plant growth. Therefore, we inferred that shortening day-length and falling temperatures in autumn can cause up-regulation of miR164a, miR164d, miR156g-3p, and miR156e. This might trigger SPL transcription factors and NAC1 genes and repress growth of dormant alfalfa in autumn. In addition, the down-regulation of miR172a expression only appeared in dormant alfalfa, which strengthened floral repression and repressed plant growth (four members of miRNAs belonging to the 172miRNA family were identified). Non-dormant alfalfa continues normal growth in the autumn, and the increased expression of miR172a in three other non-dormant situations relieved floral repression and promoted plant growth.
In Medicago truncatula, the regulation of genes encoding copper proteins by miR398a/b and miR408 suggests a link between copper homeostasis and Medicago truncatula adaptation, and the up-regulation of miR398a/b and miR408 and the clear down-regulation of their respective target genes [10]. In alfalfa, among the predicted up-regulation miRNAs in three sets (DS3 vs DM1, DS3 vs NM2 and DS3 vs NS4), miR5228 expression was at the highest level by high-throughput sequencing and RT-PCR. The predicted functions of miR5228 at the transcriptome level of Medicago sativa focused on serine/threonine-protein phosphatase, Tobacco Mosaic Virus (TMV) resistance protein N-like, apoptosis inducing factor homolog A-like, and coatomer subunit beta-1-like (Table S5). Phytochrome acts as a serine kinase and can autophosphorylate itself or transphosphorylate its partner proteins. The exposure of light to plants triggers translocation of the phytochrome from the cytosol to the nucleus where it interacts with several proteins, including transcriptional factors [64]. Alterations in phytochrome A expression disrupt the circadian clock under low light intensities and inhibit the perception of flower-inducing long day conditions in Arabidopsis thaliana [65,66]. Wang et al. also demonstrated that the photoperiod effects on phytochrome B and abscisic acid (ABA) in three alfalfa varieties differing in degree of fall dormancy [67]. At the genome level of Medicago truncatula, miR5228 predicted function was sentrin-specific protease chromo-some5(chr5) in accordance with the transcriptome level of Medicago sativa. This result suggests that miR5228 might play an important role in post-transcriptional gene regulation of fall dormancy in alfalfa. Ubiquitination is an enzymatic, posttranslational modification (PTM) process in which a ubiquitin protein is attached to a substrate protein. From KEGG pathway analysis, the main factors of fall dormancy include ubiquitin-mediated proteolysis. Ubiquitination is involved in almost all organism activity, such as cell cycle control, reproductive death, gene expression, and the regulation of transcription. In this way, we inferred that increased levels of miR5228 may reduce alfalfa growth during fall dormancy. The current data show that the expression of miR5228 in dormant alfalfa (DS3) is significantly higher than in dormant alfalfa (DM1) and non-dormant alfalfa (NM2, and NS4). These results may explain why organ growth slows down in dormant alfalfa. Based on de novo transcriptome sequencing and identification of fall dormancy related genes in alfalfa, serine/threonine-protein phosphatase and ubiquitination were the main factors to target genes regulating fall dormancy functions. For the rest of the predicted up-regulation miRNAs in three sets of alfalfa, miR5279 had no hits, and the other microRNAs (miR2674, miR5299, miR2643a, miR5205b, and miR2655b) had hypothetical proteins, resistance proteins and predicted proteins. For example, miR2643a predicted functions include acting as an IAA-amino acid hydrolase ILR1-like protein, F-box/kelchrepeat protein, F-box family protein, phototropin, and heat shock protein. Later, Fan et al. and Du et al. found a negative correlation between IAA content and fall dormancy in dormant alfalfa leaves under a short photoperiod and falling temperature [68,69]. The IAA-amino acid conjugate plays a critical role in plant growth, and it is inferred that increased levels of miR2643a may reduce the content of the IAA-amino acid conjugate and further lead to stopping growth of alfalfa. F-box protein is an expanding family of eukaryotic proteins characterized by an F-box motif, which has specificity of substrate recognition in the ubiquitinmediated proteolysis. These proteins have been shown to be critical for many physiological processes, such as cell cycle transition, signal transduction, and development [70,71]. Phototropin is one kind of photoreceptor whose wavebands of the solar spectrum are 320-500 nm. It belongs to serine/threonine kinase with flavin as a chromophore and it has phototropism. Fall dormancy is affected by the photoperiod and is further influenced by phototropin. The miR2674 function annotation is resistance protein, and this was consistent with fall dormancy being a kind of biological regulation of response to stress from the analysis of GO enrichment and KEGG pathways.
In Medicago truncatula, ethylene-responsive miRNAs in roots have been identified by high-throughput sequencing. Eight miRNAs were shown to be down-regulated after exposure to ethylene, and the potential role of these miRNAs in the ethylene-induced inhibition of root elongation was discussed [72]. Among the predicted down-regulation of miRNAs in three sets of alfalfa, miR52671 and miR2629h had no hits, and only predicted a common function as a kind of hypothetical protein. MiR2676a can function as the mediator of an RNA polymerase II transcription subunit, F-box/WD-40 repeat-containing protein chr5, small ubiquitin-related modifier, and cyclin-like F-box calcium-binding mitochondrial carrier protein. MiR2593e can react on NAC-domain protein chr3, MiR5286b, miR5286a and miR2634, have similar functions with some upregulating miRNAs in the regulation of biological processes of alfalfa.
MiRNA study in Medicago truncatula is advancing quickly, however, miRNA study in Medicago truncatula is not as detailed in contrast to Arabidopsis and rice. Because most miRNAs have been identified and finding novel miRNAs from Medicago truncatula is increasingly difficult, researchers are beginning to pay more attention to the functional study of important miRNAs in diverse biochemical and physiological processes, which will become the trend of miRNA study in Medicago truncatula [57]. A more effective approach to gene silencing has been to produce sRNAs in alfalfa. Verdonk and Sullivan indicated that amiRNA silencing constructs of miR391a can be functional in alfalfa to specifically silence gene expression [33].
Some special miRNAs were found to be involved in fall dormancy of alfalfa, based on the analysis of GO classification and the KEGG pathway. Fall dormancy strongly correlates to the level of alfalfa yields. The growth of fall dormant alfalfa varieties often surpasses those of non-fall dormant alfalfa varieties in winter, but enters dormancy much earlier than non-fall dormant varieties. In this study, some novel miRNAs candidates were found to be specific to dormant alfalfa and nondormant alfalfa. These novel miRNA candidates included novel_1, novel_3, novel_42, novel_45, novel_17, novel_75, and novel_76, whose expression levels were both high and significantly different in three sets, five of which were confirmed by Northern blotting analysis. It has been suggested that these miRNAs may play a vital role in response to fall dormancy in alfalfa. However, their functions are not yet fully understood and need further study.

Conclusions
In conclusion, four sRNA libraries were generated and sequenced from alfalfa leaves in two standard varieties, at two time points, using high-throughput sequencing technology. Our study generated 7,536,934, 6,572,924, 6,548,094, and 6,627,260 clean reads, and identified 566 and 583 potential miRNAs related to fall dormancy in standard variety (Maverick and CUF101) at the genome level of Medicago truncatula and transcriptome of Medicago sativa. Based on the analysis between GO classification and the KEGG pathway, we inferred that the regulation mechanism of fall dormancy was biological regulation of responses to stress in the metabolic pathway. We also identified 28miRNAs associated with fall dormancy in standard variety (Maverick and CUF101), including 20 known miRNAs and eight novel miRNA candidates. The known fall dormancy-responsive miRNAs were found to be involved in diverse cellular processes in plants, including development, flowering, dormancy, transcription, protein degradation, cell cycle transition, signal transduction, and adaptation. These miRNA-mediated networks could play crucial roles during the dormancy of alfalfa, and our miRNA data provides valuable information regarding further functional analysis of miRNAs involved in fall dormancy in alfalfa.