Identification of Drought-Responsive MicroRNAs from Roots and Leaves of Alfalfa by High-Throughput Sequencing

Alfalfa, an important forage legume, is an ideal crop for sustainable agriculture and a potential crop for bioenergy resources. Drought, one of the most common environmental stresses, substantially affects plant growth, development, and productivity. MicroRNAs (miRNAs) are newly discovered gene expression regulators that have been linked to several plant stress responses. To elucidate the role of miRNAs in drought stress regulation of alfalfa, a high-throughput sequencing approach was used to analyze 12 small RNA libraries comprising of four samples, each with three biological replicates. From the 12 libraries, we identified 348 known miRNAs belonging to 80 miRNA families, and 281 novel miRNAs, using Mireap software. Eighteen known miRNAs in roots and 12 known miRNAs in leaves were screened as drought-responsive miRNAs. With the exception of miR319d and miR157a which were upregulated under drought stress, the expression pattern of drought-responsive miRNAs was different between roots and leaves in alfalfa. This is the first study that has identified miR3512, miR3630, miR5213, miR5294, miR5368 and miR6173 as drought-responsive miRNAs. Target transcripts of drought-responsive miRNAs were computationally predicted. All 447 target genes for the known miRNAs were predicted using an online tool. This study provides a significant insight on understanding drought-responsive mechanisms of alfalfa.


Introduction
Alfalfa (Medicago sativa L.) is an important forage species, with high nutritional quality and high yield [1]. As a legume plant, alfalfa is capable of fixing nitrogen to nitrate in nodules, by establishing a symbiotic relationship with Rhizobium in the root system, making it an ideal crop for sustainable agriculture. Thus, planting alfalfa can also improve the condition of the soil and reduce the usage of fertilizer. Alfalfa also has the potential to become a crop for bioenergy resources, with many appropriate attributes including high biomass yield potential [2,3]. However, the yield of alfalfa is often constrained by diverse abiotic stresses, including drought. Drought is a common environmental stress, affecting plants productivity [4,5]. Understanding the molecular mechanism of the alfalfa response to drought stress is a necessary step towards improving the drought tolerance of alfalfa.
Plants respond to stresses by regulating the expression of specific genes to avoid or minimize cellular damage, in order to adapt to stress conditions [6]. Gene regulation occurs at multiple levels, miRNAs in alfalfa. MiRNA qRT-PCR was also adopted for validation of the expression of selected miRNAs, which was examined by high throughput sequencing. Furthermore, characterization of target genes was performed by using bioinformatic approaches. Through high-throughput sequencing and bioinformatics analysis, known drought stress-responsive miRNAs and miRNA candidates in alfalfa were identified. This study will be very helpful for understanding post-transcriptional regulation under drought stress in alfalfa, and improving the drought tolerance of alfalfa and other legumes.

Plant Materials and Experiment Design
Medicago sativa L. cv. Aohan was used in this study. This cultivar was kindly provided by Dr. Liqiang Wan (Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing, China). Alfalfa plants were grown in pots with a diameter of three inches, containing a mixture of sand:vermiculite (1:1 v/v) in a growth chamber at 25-28 • C under a 16 h light/8 h dark photoperiod. Plants were supplied daily with Murashige-Skoog (MS) nutrient solution. Alfalfa plants at the age of eight weeks were then randomly separated into two groups, namely drought treatment and control groups. Drought stress treatment was imposed by withholding water supply for 10 days. The control plants received normal watering throughout the experiment. The volumetric water content of soil (VWCS) was detected before sampling by using a WET Sensor (Delta-T Devices Ltd., Cambridge, UK) which is a soil moisture sensor. The VWCS of control group was approximately 45%, and the VWCS of treatment group at the fifth day was 26.7%, at the tenth day this was 16.8% in the stress treatment.

Total RNA Isolation
Root and leaf samples from both drought and control plants were collected at the fifth and tenth days during the stress treatment, respectively. For each sample of leaf or root, three biological replicates were prepared, with each biological replicate collected from 10 plants. Samples were fast frozen in liquid nitrogen, stored at −80 • C.
Total RNA samples were extracted and then were prepared for sequencing, reverse transcription PCR and qRT-PCR. Equal quantities of RNA isolated from leaves and/or roots at each stress stage were pooled, using HiPure Plant RNA Mini Kit (Magen, Shanghai, China) according to the manufacturer's instructions.

Small RNA Library Construction and High-Throughput Sequencing
A total of four groups of RNA samples (WL: leaves with watering, WR: roots with watering, DL: leaves with drought stress, and DR: roots with drought stress), each with three biological replicates, were prepared. For each group, RNA at both five days and 10 days were equally pooled to make one sample. Thus, a total of 12 samples were used to construct sRNA libraries. Since we pool samples at different stress treatment time points (five days and 10 days), only robust and consistent responses could be detected.
Total RNA was isolated by 15% polyacrylamide gel electrophoresis, and RNA molecules that were less than 50 nt in length were enriched and ligated with proprietary adapters. The RNA samples ligated with adapters were reverse-transcribed and amplified by PCR to produce sequencing libraries. The 12 sRNA libraries from alfalfa leaves and roots were sequenced on an Illumina Hiseq 2500 (Santiago, CA, USA) platform at the Guangzhou RiboBio Company, China. The raw data has been deposited in the Sequence Read Archive of NCBI, with a SRA data study accession number of SRP094823.

Identification of Known and Novel MicroRNAs
The raw sequencing reads were processed to obtain unique sequences and read count/unique reads as per the procedure reported by Hackenberg et al. [36]. First, sRNA reads of 17-45 nt were annotated to Rfam databases (Rfam 11.0, rfam.janelia.org) [42], to identify and eliminate transfer RNA (tRNA), ribosomal RNA (rRNA), small nuclear RNA (snRNA) and small nucleolar RNA (snoRNA) sequences from the sRNA reads. Then we computed the rest of sequences for sequence redundancy, and mapped these sequences to miRBase (release 21, http://www.mirbase.org/) [43] without mismatches to identify known miRNAs. After removal of the known miRNAs, the remaining sequences were used to predict the novel miRNAs. The unique sequences were mapped to the M. truncatula genome version 4.0 (http://www.medicagohapmap.org/?genome) using Burrows-Wheeler Alignment (BWA) [44] to get pre-miRNA sequences for prediction of novel miRNAs.
Novel miRNAs were predicted by using Mireap [45]. Novel miRNA candidates were identified according to the criteria reported by [46]. The normalization of reads count and calculation of log Fold change were processed as [21] described.

MicroRNA Validation by qRT-PCR
The same RNA samples used for Illumina sequencing were employed in qRT-PCR analysis. Total miRNA was reverse transcribed to complementary DNA (cDNA) using the miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China). According to the manufacturers' instructions, the miRNAs were polyadenylated and reverse transcribed in one step using miRNA RT Enzyme Mix (E. coli Poly(A) Polymerase, RTase and RNasin). The Universal RT Primer was provided in the kit. Then the first-strand cDNA was prepared for qRT-PCR analysis.
qRT-PCR was performed on an Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The reaction system was constructed according using the miRcute miRNA qRT-PCR Kit (Tiangen, Beijing, China) containing SYBR ® Green detection reagents (Applied Biosystems). The cycling parameters were set according to manufacturers' recommendations. Briefly, the cycling parameters were: initial polymerase activation step for 15 min at 95 • C, 40 cycles for 20 s at 94 • C for denaturation, 34 s at 60 • C for annealing and elongation, followed by a disassociation stage. The forward primers were designed according to the miRNA sequences of interest and synthesized by Invitrogen (Carlsbad, CA, USA). The sequences of the forward primers are supplied in Table S8. The melting curves of the PCR products can be found in Figure S1. The Universal qPCR Primer was provided in the kit. The transcript abundance of each miRNA was normalized to U6 snRNA, and the 2 −∆∆Ct method was used to calculate relative expression of miRNAs [47]. In order to compare pair-wise differences in expression, a Student's t-test was performed by using Statistical Analysis Software (SAS) program.

MicroRNA Target Prediction and Function Analysis
Target genes of drought-responsive miRNAs in alfalfa were predicted using the psRNATarget online tool (http://plantgrn.noble.org/psRNATarget/). Gene annotation can also be accomplished by using this online tool. psRNATarget is a modified version of miRU. The M. truncatula spliced transcript sequences 4.0 V1 was selected as the transcript library for target search. Mature miRNA sequences responsive to drought and identified in alfalfa roots and leaves, were used as custom miRNA sequences. Default parameters for target prediction were used.

Overview of Small RNAs from Alfalfa via High-Throughput Sequencing
A total of 12 sRNA libraries comprising of four samples (WL: leaves with watering, WR: roots with watering, DL: leaves with drought stress, and DR: roots with drought stress) were generated using the Illumina HiSeq 2500 platform, each with three biological replicates. In order to obtain high quality data sets, adaptors and low-quantity reads were removed, and 12 million to 16 million clean reads at 17-45 nt in length were obtained from each of the 12 libraries. The details of raw reads and clean reads for each library are shown in Table 1. We analyzed common/specific sequences between four groups (WL, DL, WR and DR) for the total sRNA sequences. There were 40.6% and 31.1% specific sequences in group DL and DR, respectively (Figure 1), which indicates that there are changes occurring on the molecular level. These changes may be caused by drought stress, or by other phenomena such as plant development.  In order to classify the sRNA sequenced reads into different categories and identify all the miRNA sequences in the 12 libraries, we mapped the reads to specific databases. Generally, the sequenced sRNA reads mapped to the miRBase database were abundant in drought samples (DL 10.09%; DR 12.53%) compared to control samples (WL 5.93%; WR 9.04%). The percentage of the sRNA reads mapped to the tRNA database was greater in drought samples (leaves 12.71%; roots 10.98%) than in control samples (leaves 8.98%; roots 5.14%). For rRNA, the opposite was true. Both drought and control samples had large numbers of unannotated reads (WR 63.71%, WL 50.74%, DR 60.08%, DL 55.3%) ( Figure 2). The number of total sequences that matched M. truncatula genome is given in Table 1. Sequences failing to map to the genome ranged from 5.5% to 17.3%, with exception of WR (30.9%) and DR (40.1%). These unmapped reads may have been due to unavailable genome or sequencing errors. In order to classify the sRNA sequenced reads into different categories and identify all the miRNA sequences in the 12 libraries, we mapped the reads to specific databases. Generally, the sequenced sRNA reads mapped to the miRBase database were abundant in drought samples (DL 10.09%; DR 12.53%) compared to control samples (WL 5.93%; WR 9.04%). The percentage of the sRNA reads mapped to the tRNA database was greater in drought samples (leaves 12.71%; roots 10.98%) than in control samples (leaves 8.98%; roots 5.14%). For rRNA, the opposite was true. Both drought and control samples had large numbers of unannotated reads (WR 63.71%, WL 50.74%, DR 60.08%, DL 55.3%) ( Figure 2). The number of total sequences that matched M. truncatula genome is given in Table 1. Sequences failing to map to the genome ranged from 5.5% to 17.3%, with exception of WR (30.9%) and DR (40.1%). These unmapped reads may have been due to unavailable genome or sequencing errors.

Identification of Known miRNAs
Known miRNAs in alfalfa were identified by mapping sRNA sequences generated from each library to the miRNAs database (miRBase 21, released in June 2014). After a homology search and removal of miRNAs with expression levels less than 10, 287, 314, 204 and 142 miRNAs were identified from the groups WL, DL, WR and DR, respectively. There were 348 known miRNAs belonging to 80 miRNA families that were identified from the 12 libraries. The details of miRNAs of each library are listed in Tables S1 and S2. Among these miRNA families, the miR159 and miR166 families had the most reads, with exception of the replicate WR1 (miR166 and miR398 families). Of these identified miRNAs, miR166 contained the most members, including miR166a-g, miR166i and miR166u. Additionally, the most abundant miRNA was miR5213-5p followed by miR166a-3p in the DR1, DR2 and DR3 libraries. In the other libraries, miR166a-3p was the most abundant, followed by miR5213-5p or miR159 (see Tables S1 and S2).

Identification of Novel miRNAs
A total of 281 novel miRNAs were identified from the 12 sRNA libraries using Mireap software. For each library, detailed information of predicted novel miRNAs are listed in Table S3, including novel miRNAs sequences, reads length, reads number, GC contents, pre-miRNA sequences, miRNAs loci and pre-miRNA length. The read counts of these novel miRNAs ranged from 10 to 902. A total of 26 out of 281 novel miRNAs were sequenced over 100 times, while only six novel miRNAs were sequenced over 500 times. Eight novel miRNAs were found to exist in at least five libraries ( Table 2), indicating that they are miRNA candidates with a higher level of confidence. Their precursor sequences, as well as the stem-loop hairpin secondary structure are shown in Figure 3.

Identification of Known miRNAs
Known miRNAs in alfalfa were identified by mapping sRNA sequences generated from each library to the miRNAs database (miRBase 21, released in June 2014). After a homology search and removal of miRNAs with expression levels less than 10, 287, 314, 204 and 142 miRNAs were identified from the groups WL, DL, WR and DR, respectively. There were 348 known miRNAs belonging to 80 miRNA families that were identified from the 12 libraries. The details of miRNAs of each library are listed in Tables S1 and S2. Among these miRNA families, the miR159 and miR166 families had the most reads, with exception of the replicate WR1 (miR166 and miR398 families). Of these identified miRNAs, miR166 contained the most members, including miR166a-g, miR166i and miR166u. Additionally, the most abundant miRNA was miR5213-5p followed by miR166a-3p in the DR1, DR2 and DR3 libraries. In the other libraries, miR166a-3p was the most abundant, followed by miR5213-5p or miR159 (see Tables S1 and S2).

Identification of Novel miRNAs
A total of 281 novel miRNAs were identified from the 12 sRNA libraries using Mireap software. For each library, detailed information of predicted novel miRNAs are listed in Table S3, including novel miRNAs sequences, reads length, reads number, GC contents, pre-miRNA sequences, miRNAs loci and pre-miRNA length. The read counts of these novel miRNAs ranged from 10 to 902. A total of 26 out of 281 novel miRNAs were sequenced over 100 times, while only six novel miRNAs were sequenced over 500 times. Eight novel miRNAs were found to exist in at least five libraries ( Table 2), indicating that they are miRNA candidates with a higher level of confidence. Their precursor sequences, as well as the stem-loop hairpin secondary structure are shown in Figure 3.

Drought-Responsive miRNAs Identified in Alfalfa
To identify drought-responsive miRNAs, miRNAs which absent from two of three biological replicates were filtered out, then the normalized expression profiles of known miRNAs in drought-treated samples were compared to the control samples using generalized linear model analysis with the edgeR package [48]. A log Fold Change (logFC) change cut-off of 1 and a p-value ≤ 0.05 were used to obtain the differentially expressed miRNAs [49]. These differentially expressed miRNAs between the control and drought treatment group were called as "Drought-responsive miRNAs".
A total of 12 and 18 miRNAs were observed responsive to drought treatment in alfalfa leaves and roots, respectively (Figure 4). Some of drought-responsive miRNAs reported by previous studies were also detected in our analysis. For example, miRNAs miR166 and miR398 were found to

Drought-Responsive miRNAs Identified in Alfalfa
To identify drought-responsive miRNAs, miRNAs which absent from two of three biological replicates were filtered out, then the normalized expression profiles of known miRNAs in drought-treated samples were compared to the control samples using generalized linear model analysis with the edgeR package [48]. A log Fold Change (logFC) change cut-off of 1 and a p-value ≤ 0.05 were used to obtain the differentially expressed miRNAs [49]. These differentially expressed miRNAs between the control and drought treatment group were called as "Drought-responsive miRNAs".
A total of 12 and 18 miRNAs were observed responsive to drought treatment in alfalfa leaves and roots, respectively (Figure 4). Some of drought-responsive miRNAs reported by previous studies were also detected in our analysis. For example, miRNAs miR166 and miR398 were found to be down-regulated in alfalfa roots, while miR319 and miR157 were up-regulated in both roots and leaves. Detail information of the drought-responsive miRNAs such as log fold-change, counts per million, and p-value, are given in Tables S4 and S5. MiR396, miR159/319, miR160, miR482, miR157 and miR1507 in alfalfa roots were also found to be drought-responsive. MiR156, miR3512, miR5368, miR3630, miR6173, miR5213 and miR5294 were drought-responsive in alfalfa leaves. be down-regulated in alfalfa roots, while miR319 and miR157 were up-regulated in both roots and leaves. Detail information of the drought-responsive miRNAs such as log fold-change, counts per million, and p-value, are given in Tables S4 and S5. MiR396, miR159/319, miR160, miR482, miR157 and miR1507 in alfalfa roots were also found to be drought-responsive. MiR156, miR3512, miR5368, miR3630, miR6173, miR5213 and miR5294 were drought-responsive in alfalfa leaves. Our results also revealed that miR396, miR159, miR160, miR482, and miR1507 were down-regulated in alfalfa roots; miR3512, miR5368, miR3630 and miR6173 were down-regulated in leaves, whereas miR156, miR157, miR159/319, miR5213 and miR5294 were up-regulated in alfalfa leaves. The expression level of some miRNAs showed inconsistency and very high variation across replications. These miRNAs may have been be falsely identified by software, as the algorithm for discovering differential expressed miRNAs is subject to false positives. Validation of drought-responsive miRNAs by qRT-PCR experiments was therefore necessary. It is worthy to note that the expression of the homologues miRNAs belonging to the same miRNA family was consistently similar. For instance, in alfalfa roots, aqc-miR166a, hbr-miR166a and aly-miR166a-3p were all down-regulated.
We also detected other drought related miRNAs, but the expression levels of these miRNAs did not show significant changes, such as miR168, miR393, miR408 and miR2118. One of the reasons may have been due to the approach of sample pooling used in this study. The expression level of these miRNAs may have changed significantly at 5 days and/or 10 days, but we could not detect such changes in this experiment, because the significant change of expression level at one time-point could be canceled out by the expression level at another time-point after pooling the samples at five days and 10 days together.

Validation of Drought-Responsive miRNAs by qRT-PCR
To experimentally validate a number of selected miRNAs detected from the Illumina high-throughput sequencing, qRT-PCR technique was employed. Fifteen known miRNAs were Our results also revealed that miR396, miR159, miR160, miR482, and miR1507 were down-regulated in alfalfa roots; miR3512, miR5368, miR3630 and miR6173 were down-regulated in leaves, whereas miR156, miR157, miR159/319, miR5213 and miR5294 were up-regulated in alfalfa leaves. The expression level of some miRNAs showed inconsistency and very high variation across replications. These miRNAs may have been be falsely identified by software, as the algorithm for discovering differential expressed miRNAs is subject to false positives. Validation of droughtresponsive miRNAs by qRT-PCR experiments was therefore necessary. It is worthy to note that the expression of the homologues miRNAs belonging to the same miRNA family was consistently similar. For instance, in alfalfa roots, aqc-miR166a, hbr-miR166a and aly-miR166a-3p were all down-regulated.
We also detected other drought related miRNAs, but the expression levels of these miRNAs did not show significant changes, such as miR168, miR393, miR408 and miR2118. One of the reasons may have been due to the approach of sample pooling used in this study. The expression level of these miRNAs may have changed significantly at 5 days and/or 10 days, but we could not detect such changes in this experiment, because the significant change of expression level at one time-point could be canceled out by the expression level at another time-point after pooling the samples at five days and 10 days together.

Validation of Drought-Responsive miRNAs by qRT-PCR
To experimentally validate a number of selected miRNAs detected from the Illumina high-throughput sequencing, qRT-PCR technique was employed. Fifteen known miRNAs were tested by qPCR, and the results suggested that ahy-miR398, aau-miR396, mtr-miR1507-3p, aqc-miR166a, aly-miR166a-3p, ahy-miR3512, han-miR3630-3p and mtr-miR156g-3p showed similar expression patterns to those revealed by high-throughput sequencing analysis. However, we found that the expression level of aly-miR396b-5p as detected by qRT-PCR was inconsistent with that of our sequencing results, and the p-values of gma-miR319d, bdi-miR159a-3p, nta-miR482a, gma-miR5368 and mtr-miR5294a from qRT-PCR analysis showed discrepancies to the high-throughput sequencing data, these may have been caused by sequencing error or other reasons. The 15 known miRNAs are shown in Table 3. The qRT-PCR results suggest that our sequencing data are credible.

Targets of Drought-Responsive miRNAs and Their Functional Analysis
The mature sequences of the 18 miRNAs in roots and 12 miRNAs in leaves that were modulated by drought, were used to search for their targets in alfalfa (Tables S4 and S5). Target genes of some of the miRNAs identified in this study are already known, such as miR166, miR159/319, miR160, miR396, miR398, miR482, mir156 and miR157. As for the drought responsive miRNAs of which targets are unknown, we used the online tool psRNATarget to predict their targets, by matching the miRNAs to Mt 4.0, and the results are shown in Table 4.
A total of 445 target genes for the drought-responsive miRNAs, and 196 target genes for the eight novel miRNAs were predicted, and detailed information of the target genes including target ID and functional annotation are shown in Tables S6 and S7. We noted that most of these miRNAs had more than one predicted targets. Some miRNAs had no predicted target due to lack of genome information. A number of these predicted targets were involved in metabolism, growth and response to stresses. Some predicted targets play vital roles in abiotic stress responses. For instance, ahy-miR3512 was predicted to target three genes ( Table 4), one of which was spermidine synthase (Medtr8g063940.1), which is involved in growth and resistance to adverse stresses including heat, salinity and drought.

Discussion
In this study, we constructed 12 libraries from different tissues of alfalfa treated with drought stress and well watering (control). All libraries were sequenced using the Illumina Hiseq 2500 platform. Small RNA deep sequencing from leaves and roots of alfalfa and comprehensive and systematic analysis were performed. We first identified drought-responsive miRNAs and their targets, and subsequently focused on discovering novel miRNAs.
In general, an RNAseq tag density of 2-5 million reads is sufficient for miRNA expression profiling and discovery applications [49]. In the current study, approximately 14-18 million reads for each library were generated by high-throughput sequencing. Furthermore, we sequenced three biological replicates for WL, DL, WR and DR. Our ur sRNA sequencing depth is sufficient not only for profiling of the miRNAs expression, but also for the discovery of poorly-expressed novel miRNAs.
Approximately 82-94% of the reads sequenced in alfalfa leaves mapped to the M. truncatula genome, indicating that most of the miRNAs of M. truncatula and M. sativa are identical in leaves. However, only 58-73% of the reads sequenced in alfalfa roots matched to the M. truncatula. This is similar to several previous studies, on peach [50] and on M. truncatula [25]. Known miRNAs comprised 5.1-12.9% of the mapped reads, with DL and DR having higher proportions of known miRNAs (10.1-12.5%) while WL and WR had a lower frequency (5.9-9.0%). This indicates that the drought treatment activated some miRNA-related pathways. Unannotated sequences ranged from 50.7%-63.7%, which strongly implies the existence of large amounts of undiscovered sRNAs in alfalfa.
By In silico analysis, we identified 348 known miRNAs from the 12 libraries. The miR166, miR159, miR482 and miR2118 families were abundantly expressed. MiR166a-3p and miR5213-5p were the most abundant miRNAs. These highly expressed miRNAs may play important regulatory roles in gene expression. For example, miR166 is involved in plant organism morphogenesis, such as shoot apical meristem and floral development [51].
Discovering novel miRNAs is one of the advantages of next-generation sequencing compared to other technologies such as microarrays. In this study, 281 new miRNA candidates were predicted by computational methods. However, only a few miRNAs were commonly expressed in roots and/or leaves. Most of the novel miRNAs were uniquely expressed in each library, thus we cannot profile the differential expression of the novel miRNAs among libraries. For profiling the expression of the novel miRNAs, enhanced sequencing depth may be required. Moreover, the expression level of the predicted miRNA candidates was relatively lower than conserved miRNAs, and this result is in agreement with previous reports [25,50,52,53]. One possible explanation for this result is that the conserved miRNAs regulate target genes which may be involved in many important metabolic processes in Viridiplantae, while the expression of nonconserved miRNAs may be environmentally inducible or tissue-specific, thus the expression level of conserved miRNAs may be higher than non-conserved miRNAs [52,53].
Besides identifying known and novel miRNAs, high-throughput sequencing technology also provides an alternative way for evaluating the expression of miRNA genes. To obtain robust and consistently-expressed drought-responsive miRNAs, we pooled the samples from different time-points of stress together when constructing sRNA libraries according to [25,32,54,55]. To discover drought-responsive miRNAs from the sequencing data, the mature miRNA expression profile of drought-treated leaves and roots was compared with the control group, to identify miRNAs significantly modulated by drought. Finally, 18 known miRNAs in roots and 12 known miRNAs in leaves were screened as drought-responsive miRNAs.
In roots, 16 out of 18 drought-responsive miRNAs belonging to miR396, miR159, miR160, miR482, miR1507, miR166, miR156 and miR398 families were down-regulated. Since miRNAs negatively regulate their target genes, it can be predicted that targets of down-regulated miRNAs during drought stress may play positive roles for drought stress responses. Concurring with expectations, miR166, which down-regulates HD-ZIPIII transcription factors, was down-regulated. This indicates that the target gene of miR166 was up-regulated. HD-ZIPIII transcription factors are important for lateral leaf development, root development and axillary meristem initiation [56,57]. In Triticum dicoccoides [58] and in barley [59], miR166 was down-regulated under drought stress. Inversely, miR166 was up-regulated by drought in M. truncatula [60], which was different from our results. Similarly, Long et al. [21] previously reported that miRNAs in M. truncatula and M. sativa presented the opposite expression pattern under salt stress. The difference of expression pattern between M. truncatula and M. sativa may be caused by their different resistances to stresses, or the fact that the sRNA libraries were pooled from different time-points, which would have complicated the results.
MiR159/319 was also down-regulated in alfalfa root. The MYB family and TCP family are targets families of this miRNA, and both of them were identified by 5 RACE [18]. Under drought conditions, most active MYB transcription factors are involved in ABA signaling pathway. Previous studies in Arabidopsis spp. indicate that some MYBs, including MYB2, are positive regulators of ABA signaling [61,62]. With the decreased expression of miR159 resulting in the positive regulation of the ABA signaling pathway, down-stream of ABA signaling pathways are stimulated, including root development. Similarly, miR160, which targets auxin response factors (ARFs), also has positive regulatory roles in drought stress responses [63].
Known targets of miR398 are involved in respiration and oxidative stress [64]. We found that drought stress down-regulated miR398 in alfalfa roots, in accordance with the results in maize [65] and M. truncatula [27]. However, Trindade et al. [60] found the opposite in M. truncatula. The differences in the expression of miR398 showed here may be caused by differences in species responses, duration of drought stress, and the metabolic states of the individual plants in different studies [5].
MiR396, miR482 and miR1507 were down-regulated under drought stress. However, their currently identified roles only include development or disease resistance; therefore, it can be predicted that they may have additional targets that are yet to be identified.
Two of these drought-responsive miRNAs (gam-miR319d and aly-miR157a-3p) were up-regulated in alfalfa roots. In order to conserve water and protect the cell, miRNAs are expected to be up-regulated during drought stress, so that those processes involved in normal growth and metabolism can be shut down. However, we failed to obtain valuable information about stress resistance from their target gene annotation. Interestingly, as a member of miR159 family, gam-miR319d was expected to be down-regulated by drought, to cause the positive regulation of the ABA signaling pathway, but it did not show down-regulation in our study. However, the target gene of gam-miR319d, was the NB-ARC domain protein, a disease resistance protein, instead of the MYB transcription factor, as predicted using the online program psRNATarget. There are two possible reasons to explain this result: (1) gam-miR319d may have other unknown targets that play negative roles in drought adaption in alfalfa; (2) gam-miR319d plays no role in drought adaption, but can be stimulated by drought through an unknown mechanism. With the function of gam-miR319d still being unclear, further studies will be needed to explore the roles of gam-miR319d in drought stress.
In alfalfa leaves, most of the drought-responsive miRNAs were involved in development, substance synthesis and transport. MiR3512, miR5368, miR3630-3p and miR6137, whose targets remain unknown, were down-regulated in this experiment. Target genes of miR3512, miR3630-3p and miR6137 were obtained using the online program psRNATarget. These targets include zein-binding protein, polyol/monosaccharide transporter, spermidine synthase, cytochrome P450 family 71 protein and TCP family transcription factor.
MiR159/319 is down-regulated in roots and up-regulated in leaves, which may indicate that the same miRNA could play different roles in different tissues. In Arabidopsis spp., MYB33 and MYB101 transcripts are targets of miR159a [66,67]. Under drought conditions, MYB33 and MYB101 could modulate stomatal movement by regulating the ABA signal, suggesting that miR159/319 plays a positive role on drought response by decreasing stomatal conductance. Additionally, miR156, miR157, miR5213 and miR5294 were up-regulated significantly (p < 0.01), and their targets were involved in development or disease resistance, indicating these miRNAs play some roles under drought stress.
In summary, 348 known miRNAs have been identified from alfalfa leaves and roots, and a number of candidates for drought-responsive miRNAs have been identified, thus this attempt paves the path for better understanding the drought-responsive mechanisms of alfalfa. In addition, the identification of 300 novel miRNAs also provides promising resources for future research in understanding post-transcriptional regulation in alfalfa. Future studies should focus on the verification of novel miRNAs and their predicted targets by experimental approaches, and additionally, the effects of drought-responsive miRNAs on drought tolerance need to be illuminated. To our knowledge, our study is the first systematic and comprehensive identification of drought-responsive miRNAs in an alfalfa species. This study is valuable for improving drought tolerance and systems to mitigate crop losses under drought stress.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4425/8/4/119/s1. Figure S1: Melting curves of qRT-PCR products. Table S1: Members and expression abundance of known miRNAs from leaves of alfalfa treated with drought and irrigation. Table S2: Members and expression abundance of known miRNAs from roots of alfalfa treated with drought and irrigation. Table S3: Mature and precursor sequences of predicted miRNAs and the other information including reads length, reads number, GC contents and minimum free energy. Table S4: Members and expression abundance of drought-responsive miRNAs from roots of alfalfa treated with drought and irrigation. Table S5: Members and expression abundance of drought-responsive miRNAs from leaves of alfalfa treated with drought and irrigation. Table S6: Targets of drought-responsive miRNAs from leaves of alfalfa. Table S7: Targets of drought-responsive miRNAs from roots of alfalfa. Table S8: Sequences of the forward primers.