Utilization of Transcriptome, Small RNA, and Degradome Sequencing to Provide Insights Into Drought Stress and Rewatering Treatment in Medicago ruthenica

Drought is a major limiting factor in foraging grass yield and quality. Medicago ruthenica (M. ruthenica) is a high-quality forage legume with drought resistance, cold tolerance, and strong adaptability. In this study, we integrated transcriptome, small RNA, and degradome sequencing in identifying drought response genes, microRNAs (miRNAs), and key miRNA-target pairs in M. ruthenica under drought and rewatering treatment conditions. A total of 3,905 genes and 50 miRNAs (45 conserved and 5 novel miRNAs) were significantly differentially expressed in three test conditions (CK: control, DS: plants under drought stress, and RW: plants rewatering after drought stress). The degradome sequencing (AllenScore < 4) analysis revealed that 104 miRNAs (11 novel and 93 conserved miRNAs) were identified with 263 target transcripts, forming 296 miRNA-target pairs in three libraries. There were 38 differentially expressed targets from 16 miRNAs in DS vs. CK, 31 from 11 miRNAs in DS vs. RW, and 6 from 3 miRNAs in RW vs. CK; 21, 18, and 3 miRNA-target gene pairs showed reverse expression patterns in DS vs. CK, DS vs. RW, and RW vs. CK comparison groups, respectively. These findings provide valuable information for further functional characterization of genes and miRNAs in response to abiotic stress, in general, and drought stress in M. ruthenica, and potentially contribute to drought resistance breeding of forage in the future.


INTRODUCTION
Drought is the most widespread climatic extreme that has a negative impact on human and natural environments (Touma et al., 2015;Schwalm et al., 2017). It has received more attention with the increase of severe drought occurrences. The physiological acclimatization of individuals may buffer the effect of drought (Schwalm et al., 2017); thus, it is meaningful to improve adaptability in plants by increasing their drought stress tolerance.
Plants have specific adaptive responses that protect them from environmental stresses (Chinnusamy et al., 2004), which are generally controlled by many complex regulatory networks involving numerous genes. The regulation of gene expression can be carried out at the transcriptional, post-transcriptional, and epigenetic modification levels. In plants, several biological processes are regulated by small RNAs at the transcriptional and post-transcriptional levels, including plant growth and development processes, as well as biotic and abiotic stress responses (Trindade et al., 2010;Capitão et al., 2011;Liu et al., 2014;Gao et al., 2016;Shu et al., 2016;Garg et al., 2019). Small RNAs (21-26 nt) are ubiquitous, versatile repressors of gene expression in plants, animals, and many fungi, which include short interfering RNAs (siRNAs), small temporal RNAs (stRNAs), heterochromatic siRNAs, tiny non-coding RNAs, and microRNAs (miRNAs) (Finnegan and Matzke, 2003). In plants, miRNAs are endogenous non-coding small RNAs measuring 20-24 nt long (Jones-Rhoades et al., 2006). They negatively regulate gene expression at the post-transcriptional level via direct cleavage of the target mRNA or inhibition of target gene translation by recognizing and combining to their target mRNAs (Bartel, 2004;Voinnet, 2009). A single miRNA can have several target genes, and several miRNAs can regulate one gene. Therefore, the gene pool regulated by miRNAs can be very extensive.
Plant miRNAs are frequently complementary to their target mRNAs; this complementarity effectively triggers target mRNA degradation (Llave et al., 2002;Tang et al., 2003) and loss of protein-coding function (German et al., 2008;Iwakawa and Tomari, 2015). Most miRNAs regulate the expression of target genes in plants via splicing, and slicing often occurs at the tenth or eleventh nucleotide of the complementary region of the miRNA and mRNA. Thus, the identification of target genes is crucial for miRNA functional analysis. High-throughput sequencing is the most popular technique for identifying miRNAs in plants, because it allows for an easier and faster access to large numbers of miRNAs, especially those of low abundance (Allen et al., 2004). Afterwards, the functions of miRNAs may be identified through bioinformatic prediction or degradome sequencing. Degradome sequencing screens out the target genes for miRNAs by combining the advantages of high-throughput sequencing technology and bioinformatics analysis. The key point of the correlation analysis of the transcriptome and small RNA is the result of the degradome sequencing.
Medicago ruthenica L. is a cross-pollinated, diploid (2n = 16) perennial legume forage, with a remarkable ability to adapt to extreme environments. Thus, M. ruthenica is considered a valuable forage crop, which may be used as a potential source of genes to improve abiotic stress resistance in cultivated alfalfa (Campbell et al., 1997). There have been many reports on the molecular mechanisms induced by drought stress in leguminous forage Zhang et al., 2018Zhang et al., , 2019, but few on M. ruthenica. Moreover, a number of miRNAs have been identified to be associated with responses to drought stress in several species, such as Arabidopsis thaliana (Liu et al., 2008), tomato (Liu et al., 2018), Oryza sativa (Chen and Li, 2018;Nadarajah and Kumar, 2019), and Gossypium hirsutum (Lu et al., 2019). Here, we integrated transcriptome, miRNAome, and degradome results to identify drought-response genes and miRNAs in M. ruthenica, and find potential regulatory patterns of miRNA-target pairs. These results may provide novel insights into the response to abiotic stresses of M. ruthenica, and contribute to drought resistance breeding of forage in the future.

Transcriptome Sequencing in M. ruthenica Under Drought Stress
To profile the expression of genes in M. ruthenica in response to drought stress, nine libraries were constructed from three leaf samples (CK: control, DS: plants under drought stress, RW: plants rewatering after drought stress), each with three biological replicates. Then,36,207,356 to 55,818,774  All assembled unigene clusters were aligned against Gene Ontology (GO), 1 Kyoto Encyclopedia of Genes and Genomes (KEGG), 2 Pfam database, 3 SwissProt database, 4 evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) 5 databases, and NCBI non-redundant protein database (NCBI_NR) 6 using DIAMOND 23 with a threshold E-value < 0.00001 (Buchfink et al., 2015). The statistical results from six authoritative databases are listed in Table 2.

Differentially Expressed Genes in M. ruthenica Under Drought Stress
To identify the differentially regulated genes under drought stress in M. ruthenica, three comparisons of the three test conditions (CK, DS, and RW) were performed. In total, 3,905 genes were significantly differentially expressed among the three test  groups. There were 3,065, 2,508, and 205 differentially expressed genes (DEGs) in the DS vs. CK, DS vs. RW, and RW vs. CK comparisons, respectively ( Figure 1A). The overlapping DEGs among the three comparisons were shown in a Venn diagram in Figure 1B.  Figure 1C); this indicates that the 1,736 overlapping DEGs were involved in drought stress responses in M. ruthenica. Among all DEGs, 10 showed differential expression across all the treatments ( Figure 1B). The information of all the DEGs in three comparisons is shown in Supplementary Table 2.

qRT-PCR Verification of mRNA-Seq Analysis of Gene Expression
To validate the mRNA-Seq results, 10 genes were selected for quantitative real-time polymerase chain reaction (qRT-PCR) analysis. The information of these genes is shown in Supplementary Table 4. The qRT-PCR results (Figure 2) showed similar expression trends to their high-throughput sequencing analysis, which suggested that our RNA-seq data are credible.

Deep Sequencing of M. ruthenica Small RNAs
The small RNA deep sequencing of M. ruthenica was performed from three samples (CK, DS, and RW), each with three biological replicates. We counted the raw data sequencing outputs, and generated unique sequences and their corresponding copies. A total of 16,112,539, 13,520,602, and 17,459,447 reads, and 4,029,151, 2,602,204, and 3,314,079 unique sequences were generated from CK, DS, and RW, respectively (Supplementary Table 5). Meanwhile, by comparing the acquired sgRNA sequence with the mRNA, RFam (including rRNA, tRNA, snRNA, and snoRNA), and Repbase databases, we established 12,257,093, 6,163,546, and 10,663,527 valid reads, and 3,422,624, 1,675,939, and 2,468,514 valid unique sequences in the CK, DS, and RW groups. These valid data were subjected to further miRNA comparison identification and prediction analysis.

Identification of Conserved and Novel miRNAs in M. ruthenica
To identify conserved miRNAs in M. ruthenica, we compared the small RNA sequences with known plant miRNAs in the miRBase. Based on sequence homology, 532 conserved miRNAs belonging to 65 miRNA families, and 63 novel miRNAs were identified from the nine libraries (Table 5).

Drought-Responsive miRNAs in M. ruthenica
To identify miRNAs in response to drought in M. ruthenica, we analyzed and compared the read counts in nine libraries. miRNAs with a p-value < 0.05 were considered DEMs. In total, 33, 21, and 40 miRNAs were differentially expressed in the DS vs. CK, DS vs. RW, and RW vs. CK comparisons, respectively. Meanwhile, respective up-and downregulation profiles were found in each comparison: 3 and 30 miRNAs in DS vs. CK, 5 and 16 in DS vs. RW, and 21 and 19 in RW vs. CK (Figure 3). The details of these DEMs found in the three comparisons are shown in Supplementary Table 6.
Our results also revealed that 50 miRNAs showed significant differentially expressed patterns among the RW vs. DS vs. CK comparison, including 45 conserved miRNAs and 5 novel miRNAs ( Table 6). Among the 50 DEMs, the expression of 23 miRNAs, such as gma-miR171j-5p, lja-miR390a-3p, mtr-miR398a-5p, and bra-MIR408-p5, was downregulated, and then upregulated during drought stress and rehydration, respectively. On the other hand, the expression of 12 miRNAs, such as mtr-miR156b-5p, mtr-MIR156e-p3, gma-miR159a-3p, and lus-miR396b, exhibited the opposite expression profiles during the FIGURE 2 | Relative gene expression for qRT-PCR verification of mRNA-Seq results. The data were the average of three qRT-PCR replicates for each sample from three biological replicates. GAPDH of alfalfa was used as an internal reference. Error bars indicate the standard deviation of three biological replicates. same respective conditions. Moreover, 12 miRNAs, such as mtr-miR396a-5p, mtr-miR398a-3p, mtr-miR398b, mtr-miR408-3p, were always downregulated, and three miRNAs (gma-miR6300, gp1: Reads were mapped against the miRNAs/pre-miRNAs of specific species in miRBase, and the pre-miRNAs were further mapped to the genome and EST; gp2a: Reads were mapped against the miRNAs/pre-miRNAs of selected species in miRBase. The reads (as well as the miRNAs of the pre-miRNAs) were mapped against the genome, while the mapped pre-miRNAs were not. The extended genome sequences from the genome loci may form hairpins; gp2b: Reads were mapped against the miRNAs/pre-miRNAs of selected species in miRBase. The reads (as well as the miRNAs of the pre-miRNAs) were mapped to the genome, while the mapped pre-miRNAs were not. The extended genome sequences from the genome loci may not form hairpins; gp3: Reads were mapped against the miRNAs/pre-miRNAs of selected species in miRBase. The mapped pre-miRNAs and the reads were not mapped against the genome. However, the reads were mapped against the miRNAs (Matures); gp4: Reads were not mapped to pre-miRNAs of selected species in miRBase. However, the reads were mapped against the genome. The extended genome sequences may form hairpins.

Degradome Sequencing Analysis
Generally, miRNAs inhibit protein synthesis either by translational repression and/or mRNA target degradation (Eulalio et al., 2008;Filipowicz et al., 2008;Chekulaeva and Filipowicz, 2009 Table 7). For data analysis, CleaveLand4 was used to identify sliced miRNA targets detected in our study. The abundance of raw tags was plotted for each target transcript, and the cleaved target transcripts were classified into five categories (0, 1, 2, 3, and 4) according to the abundance of tags at the site position of target transcript. The numbers of cleaved target transcripts in categories 0-4 were 115,33,1,720,958,and 2,260 in CK library,97,44,2,178,1,284,and 2,716 in DS library,and 108,36,2,234,1,267,and 2,866 in RW library, respectively (Supplementary Table 8).

Identification of miRNA Targets via Degradome Sequencing in M. ruthenica
We obtained 11,390 miRNA-target pairs from the degradome sequencing of three libraries totally (Supplementary Table 8).
AllenScore reflects the penalty of mismatched bases between the miRNA and its target, and is an important measure to evaluate the matching rate between the miRNA and its target (Addo-Quaye et al., 2009;Liu et al., 2016). In our study, analysis to degradome sequencing (AllenScore < 4) showed that 104 miRNAs (11 novel and 93 conserved miRNAs) and 263 target transcripts were identified, forming 296 miRNA-target pairs in three libraries (Supplementary Table 9). It was significantly seen that the maximum targets were obtained for ppe-MIR169i-p5_2ss17GT19TG, which had 90 target transcripts in three libraries. We inferred that ppe-MIR169i-p5_2ss17GT19TG may play important roles in M. ruthenica.
In total, 200 genes, targeted by the identified miRNAs, were annotated through GO analysis in the three libraries (Supplementary Table 10). Among the biological processes, "regulation of transcription, DNA-templated" and "transcription, DNA-templated" were the most abundant. For cellular components, the most frequent categories were "nucleus" and "cytoplasm, " and "nucleus" was also the most enriched group of all GO categories. Among the molecular function categories, "DNA binding, " "DNA-binding transcription factor activity, " and "protein binding" were the most abundant ( Figure 4A).
Then, through KEGG analysis, 78 targets were annotated to 67 different pathways (Supplementary Table 10). ko04075 (plant hormone signal transduction) and ko04016 (MAPK signaling pathway -plant) pathways, which were associated with 14 and 7 unigenes, respectively, were the most abundant pathways (Figure 4B), indicating their significant roles response to drought stress in M. ruthenica. Figure 5 is the network plot for the miRNAs and their targets associated with ko04016, ko04075, and ko04626 (plant-pathogen interaction) pathways. The squamosa promoter binding protein-likes (SPLs) targeted by miR156, and auxin response factors (ARFs) targeted by miR160 were mainly involved in plant hormone signal transduction. The PRB1 (putative reverse transcriptase, RNA-dependent DNA polymerase) targeted by the novel miRNA PC-3p-39042_145 was related to ko04016, ko04075, and ko04626 pathways. The PYL9 (putative polyketide cyclase/dehydrase, START-like domain-containing protein), which was the target gene of ppe-MIR169i-p5_2ss17GT19TG, was related to ko04016 and ko04075 pathways. Probable WRKY transcription factor 75 (WRKY75), which was the target gene of ppe-MIR169i-p5_2ss17GT19TG, was involved in ko04016 and ko04626 pathways. Several target plots (T-plots) for these miRNAs and their target genes validated by degradome sequencing were showed in Figure 6. Three miRNA-target pairs with negative regulatory relationships were chosen for validation via qRT-PCR analysis. The qRT-PCR results were consistent with those of Illumina sequencing and confirmed that the three miRNA-target pairs all displayed reversed expression pattern in our study (Figure 7). This also suggests that our small RNAs and degradome sequencing data are reliable.

Correlation Analysis of miRNA Expression Profiles and Their Targets
Based on the degradome sequencing, we integrated the expression profiles of the miRNA and target genes from the different comparison groups, and obtained a table that summarizes the miRNA-target gene association analysis. After extracting information from the general table, we formulated a differential miRNA-differential target gene association analysis table, from which we obtained miRNA-target gene pairs with negative regulatory relationships. Corresponding to the DEMs, there were 38 differentially expressed targets from 16 miRNAs in DS vs. CK, 31 from 11 miRNAs in DS vs. RW, and 6 from three miRNAs in RW vs. CK (Supplementary Table 11). Then, 21, 18, and 3 miRNA-target gene pairs showed a reverse expression pattern in the DS vs. CK, DS vs. RW, and RW vs. CK comparison groups, respectively ( Table 7). We found that gma-miR171j-5p was significantly downregulated in DS vs. CK and DS vs. RW, while its target, TRINITY_DN16578_c0_g1, was significantly upregulated. This may imply that the miRNA-target gene pair were specifically triggered by water deficit.

DISCUSSION
Compared with other leguminous plants, little research has been conducted on the role of miRNAs in M. ruthenica. Here, three important high-throughput methods were applied to investigate the mechanism of drought resistance in M. ruthenica. The transcriptome data set was used as a reference sequence for small RNA and degradome sequencing analyses to identify the miRNAs and their target genes that might be associated with drought stress.
Among the 50 DEMs, three members (mtr-miR398a-3p/a-5p/b) of the miR398 family and two members (bra-MIR408-p5 and mtr-miR408-3p) of the miR408 family were all downregulated under drought, which was consistent with the results in pea (Jovanović et al., 2014). In contrast, another study indicated that miR398a/b and miR408 were increased due to water deficit in M. truncatula (Trindade et al., 2010). These differences may be caused by the difference in species, extent and duration of drought stress , and sensitivity of some miRNAs to subtle differences in plant growth conditions (Ferdous et al., 2015). This suggests that even conserved miRNAs might function in a species-specific manner (Kamthan et al., 2015). In our study, some miRNAs from the same family (such as miR398, miR408, and miR1527) showed differently expressed trends after rewatering. Similar results were observed in rice, wherein members of miR319 showed different degrees of up-or downregulation responses to drought (Zhou et al., 2010). It is possible that the miRNA gene regulators change their expression after rewatering, leading to changes in miRNA expression patterns (Ferdous et al., 2015).
Generally, a single miRNA can target several genes (Samad et al., 2017), and a single gene can be regulated by several miRNAs (Katiyar et al., 2015). We found that ppe-MIR169i-p5_2ss17GT19TG, which had the maximum number of targets among all miRNAs, targeted 90 transcripts in three libraries totally. miR169 is the largest miRNA family in A. thaliana and has 14 members, which can be divided into four groups according to their mature miRNA sequences: miRNA-169a, miRNA-169b/c, miRNA-169d/e/f/g, and miRNA-169h/i/j/k/l/m/n (Du et al., 2017). Asefpour Vakilian (2020) revealed the contribution of each miRNA to stress response in plants by implementing the feature selection algorithm on the constructed database. The algorithm showed that miRNA169 had the highest contribution to drought stress (Asefpour Vakilian, 2020).  also demonstrated that miR169a and miR169c were substantially downregulated due to drought stress, and that miR169a mainly regulates NFYA5 expression at the mRNA level , while miR169i mainly affects NFYA5 expression at the translational level (Du et al., 2017). In our study, a total of 78 ppe-MIR169i-p5_2ss17GT19TG targets were annotated by GO terms in three libraries (Supplementary Table 12). The nucleus (GO:0005634) involving 24 genes was the most enriched group of all GO categories, followed by the cytoplasm (GO:0005737), protein binding (GO:0005515), and plasma membrane (GO:0005886). KEGG analysis of ppe-MIR169i-p5_2ss17GT19TG targets showed that 30 genes annotated to 33 different pathways (Supplementary Table 12). Figure 8 was the network plot for the target genes of ppe-MIR169i-p5_2ss17GT19TG and their KEGG pathways, for instance, CRWN (protein CROWDED NUCLEI) involved in ko00270 (Cysteine and methionine metabolism), ko00280 (Valine, leucine, and isoleucine degradation), ko00290 (Valine, leucine, and isoleucine biosynthesis), and ko00770 (Pantothenate and CoA biosynthesis) pathways simultaneously. JKD (zinc finger protein JACKDAW) related to ko00052 (Galactose metabolism), ko00600 (Sphingolipid metabolism), ko00531 (Glycosaminoglycan degradation), ko00604 (Glycosphingolipid biosynthesis -ganglio series), and ko00511 (Other glycan degradation) pathways, which revealed the extensive regulatory roles of ppe-MIR169i-p5_2ss17GT19TG in M. ruthenica.
Degradome sequencing is based on high-throughput sequencing technology and bioinformatics analysis and avoids false positive results effectively, which makes it more suitable for plant miRNA target gene identification . Through the correlation analysis of the transcriptome-small RNA-degradome, we obtained differential miRNA-differential target gene association results and miRNA-target pairs with negative regulatory relationships. Among them, two miRNA-target pairs appeared in two comparison groups. gma-miR171j-5p was significantly downregulated in DS vs. CK and DS vs. RW, while its target TRINITY_DN16578_c0_g1 was significantly upregulated. Similarly, mtr-miR396a-5p was significantly downregulated in DS vs. CK and RW vs. CK, while its target TRINITY_DN23804_c0_g2 was significantly up-regulated. SCL6-II and SCL6 are the target genes of miR171 in A. thaliana (Lee et al., 2008), which play roles in plant root and leaf development, photochrome signaling, lateral organ polarity, meristem formation, vascular development, and stress response (Wang et al., 2010;Zhu et al., 2015;Asefpour Vakilian, 2020). gma-miR171o and gma-miR171q regulate GmSCL-6 and GmNSP2, respectively, which in turn influence the spatial and temporal aspects of soybean nodulation (Hossain et al., 2019). miRNA-396 is an important contributor to the plant stress response which targets four classes of stress resistance proteins: pathogen-related, nucleotide binding site resistance protein-like, dirigent-like, and ribonuclease-like proteins (Zhang et al., 2007). miRNA-396 targets cell cycle regulators, as well as those that control plant growth and differentiation. miRNA-396 is increased under stress and represses cell multiplication (Patel et al., 2019).

CONCLUSION
In conclusion, we elucidated the small RNAs and their target genes in M. ruthenica when subjected under drought stress and rehydration treatment though transcriptome, small RNA, and degradome sequencing. Although the complex miRNAmediated regulatory networks remain to be elucidated, these findings provide valuable information for further functional characterization of genes and miRNAs in response to abiotic stress, in general, and drought stress in M. ruthenica. More importantly, this study will serve as a foundation for future research on the functional roles of miRNAs and their target genes in legume forage.

Plant Material and Drought Treatments
Medicago ruthenica seeds were obtained from the Institute of Grassland Research, Chinese Academy of Agricultural Sciences at the Inner Mongolia Autonomous Region, China. This cultivar grows well in Inner Mongolia, with higher drought and cold resistance. The sanded M. ruthenica seeds were sterilized in concentrated sulfuric acid for 10 min, and then washed thoroughly with sterile water (Araújo et al., 2004). After incubating for 3 days at 4 • C in dark conditions, the seeds were placed on soaked filter paper in Petri dishes, and incubated in an artificial climate chamber (temperature, 25 ± 2 • C; photoperiod, 16 h light/8 h dark cycle; relative humidity, 50%). Two days later, the germinated seeds were transferred to pots (13 cm in diameter and 10 cm deep) with vermiculite and quartz sand. There were eight plants in each pot. Four-week-old seedlings were randomly divided into three groups: the control group (CK), the drought stress treatment group (DS), and the rewatering group (RW). The drought treatment period lasted for 15 days. During the experiment, CK plants were watered every 3 days, maintaining a relative soil moisture content of more than 80%, while DS and RW plants were subjected to water deprivation for 15 days, resulting in a relative soil moisture content of 20% on the 15th day. The RW plants were rehydrated to full soil saturation on the 15th day. To avoid any interference, no nutrients were added. Other growth conditions were maintained. Each treatment group had three biological replicates. Leaf samples of CK and DS were collected on the 15th day, and leaves of RW were collected 2 days after rewatering. Samples were immediately frozen in liquid nitrogen, and stored at −80 • C for later use.

Transcriptome Libraries Construction, Sequencing, and Analysis
Total RNA was isolated from M. ruthenica leaves using TRIzol reagent (Invitrogen, CA, United States) according to the manufacturer's protocol. Total RNA quantity and purity were analyzed using a Bioanalyzer 2100 (Agilent, CA, United States) and RNA 6000 Nano LabChip Kit (Agilent, CA, United States), respectively. Poly(A) RNA was obtained from total RNA (5 µg) using poly T oligo-attached magnetic beads after two rounds of purification. Following purification, the mRNA was fragmented into small pieces using divalent cations at elevated temperatures. Then, the cleaved RNA fragments were reverse-transcribed to create the final cDNA library by following the protocol for the mRNASeq sample preparation kit (Illumina, San Diego, CA, United States). The average insert size for the paired-end libraries was 300 bp (±50 bp). Pairedend sequencing was performed using an Illumina Hiseq4000 instrument (LC Sciences, United States) according to the manufacturer's instructions.
The adaptor contamination, low quality bases, and undetermined bases from raw data were removed using Cutadapt (Martin, 2011) and perl scripts in house. Then sequence quality was verified by FastQC, 7 including the Q20, Q30 and GC-content of the clean data. All downstream analyses were based on clean data of high quality. De novo assembly of the transcriptome was performed with Trinity 2.4.0 (Grabherr et al., 2011). Trinity groups transcripts into clusters based on shared sequence content. Such a transcript cluster is very loosely referred to as a "gene." The longest transcript in the cluster was chosen as the "gene" sequence (aka Unigene).

Differentially Expressed Genes Analysis
Salmon (Patro et al., 2017) was used to perform expression levels for genes by calculating TPM (Mortazavi et al., 2008). Statistically significant (p-value < 0.05) DEGs with a log 2 (fold change) > 1 or log 2 (fold change) < −1 we selected using the R package edgeR (Robinson et al., 2010). Next, GO and KEGG enrichment analyses were performed on the differentially expressed genes using in-house Perl scripts. The GO project is a bioinformatics resource on gene products and descriptions of functions (Gene Ontology Consortium et al., 2013). It creates annotations to describe the biological roles of individual gene products (e.g., genes, proteins, ncRNAs, complexes) by classifying them (Gene Ontology Consortium, 2015). The relationships between a gene product/or gene-product group to biological process, molecular function, and cellular component are oneto-many. GO ontologies can also use for annotation of geneexpression data (Gene Ontology Consortium et al., 2000). KEGG is a knowledge base for systematic analysis of gene functions in terms of the networks of genes and molecules (Ogata et al., 1999). KEGG is widely used for analyzing genomics, transcriptomics, proteomics, glycomics, metabolomics, and other high-throughput data (Kanehisa et al., 2017). The online Blast algorithm 8 can be used to analyze genes. If the significant similarities of Blast lead to the assignment of designated enzymes, these genes can be tagged in the corresponding KEGG pathways (Altermann and Klaenhammer, 2005).

Small RNAs Sequencing and Bioinformatics Analysis
Nine small RNA libraries were constructed from approximately 5 µg of total RNA using TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, United States). Then the libraries were sequenced by Illumina Hiseq 2500 (LC Sciences, United States) following the vendor's recommended protocol.
The adapter dimers, junk, low complexity, repeats, and common RNA families (rRNA, tRNA, snRNA, and snoRNA) were removed from raw reads using an in-house program, ACGT101-miR (LC Sciences, Houston, TX, United States). Subsequently, unique sequences (18-25 nt) were mapped to specific species precursors in miRBase 22.0 through a BLAST search; conserved miRNAs and novel 3p-and 5p-derived miRNAs were then identified. Length variation at both 3 and 5 ends and one mismatch inside of the sequence were allowed in the alignment. The unique sequences mapping to specific species mature miRNAs in hairpin arms were identified as known miRNAs. The unique sequences mapping to the other arm of known specific species precursor hairpin opposite to the annotated mature miRNA-containing arm were considered to be novel 3p-and 5p-derived miRNA candidates. The remaining sequences were mapped to other selected species precursors (with the exclusion of specific species) in miRBase 22.0 by BLAST search, and the mapped pre-miRNAs were further BLASTed against the specific species genomes to determine their genomic locations. The above two we defined as known miRNAs. 8 http://www.ncbi.nlm.nih.gov/BLAST/ The unmapped sequences were BLASTed against the specific genomes, and the hairpin RNA structures containing sequences were predicated from the flank 120 nt sequences using RNAfold software. 9

Detection of Differential Expressed miRNAs
The differential expression of miRNAs, based on normalized deep-sequencing counts, was analyzed by T-test. The criteria to classify DEMs between two-way pairing among the three experimental conditions (CK vs. DS; DS vs. RW; RW vs. CK) were as follows: the significance thresholds were set to be p < 0.05 in each test.

Construction and Analysis of Degradome Libraries
The degradome libraries construction was further optimized and simplified based on Ma et al. (2010). Total RNA from the three groups (CK, DS, and RW) was isolated and purified using TRIzol reagent (Invitrogen, CA, United States) according to the manufacturer's protocol. The RNA amount and purity of each sample were quantified using a NanoDropND-1000 (NanoDrop, Wilmington, DE, United States). RNA integrity was assessed using Agilent 2100 with RIN number > 7.0. Poly(A) RNA was purified from the plant total RNA (20 µg) twice using poly T oligo-attached magnetic beads. Because the three RNA cleavage products contain a 5 -monophosphate, the 5 adapters were ligated to their respective 5 ends using RNA ligase. Then, the first strand of cDNA for each mRNA was reverse transcribed using a 3-adapter random primer. Size selection was then performed using AMPureXP beads. Afterwards, the cDNA was amplified via PCR under the following conditions: initial denaturation at 95 • C for 3 min; 15 cycles of denaturation at 98 • C for 15 s, annealing at 60 • C for 15 s, and extension at 72 • C for 30 s; then a final extension at 72 • C for 5 min. The average insert size for the final cDNA library was 200-400 bp. Lastly, we performed 50 bp singleend sequencing using an Illumina Hiseq2500 (LC Bio, China) according to the manufacturer's recommended protocol.
Raw reads were subjected to ACGT101-DEG (LC Sciences, Houston, TX, United States) for data processing. This program mainly depends on CleaveLand4, a public software package for analyzing "degradome" data. Degradome data are a variant type of RNA-seq data, where the reads derive from the 5 -ends of uncapped RNAs (German et al., 2008;Addo-Quaye et al., 2009). These data can be used to identify miRNA and siRNA targets that are actively "sliced." CleaveLand4 handles several phases of degradome data analysis in a single command, including: Alignment of degradome data to the reference transcriptome, and parsing the output into a "degradome density file." The degradome density file reflects the counts of 5 positions across the transcriptome. Alignment of query miRNAs or siRNAs to the transcriptome to generate a list of potential target sites. This uses the program "GSTAr.pl" (Generic Small RNA Transcript Aligner), which ships with the CleaveLand4 program. GSTAr.pl uses RNA-RNA thermodynamic predictions instead of sequence similarity to identify potential target sites, making it much slower than generic aligners, but more sensitive in terms of finding all possible sites. Cross-referencing the degradome data with the alignments to identify slicing sites with evidence of slicing. This includes assessment of p-values.
All the identified target genes were annotated via GO (see text footnote 1) and classified using KEGG pathways (see text footnote 2).

Verification of mRNAs and miRNAs via Quantitative Real-Time PCR
The RNA-seq results for both mRNA and miRNA expression were verified via qRT-PCR. Total RNA was extracted from nine leaf samples using TRIzol reagent (Invitrogen, CA, United States). For the miRNA expression analysis, the U6 snRNA of M. truncatula was used as the reference gene. The alfalfa GAPDH gene was used as an internal reference for the mRNA expression analysis. All primers used for qRT-PCR are listed in Supplementary Table 13. qPCR was conducted using the ChamQ SYBR Color qPCR Master Mix (2X) on a LineGene9600plus Real Time PCR instrument. Three biological and three technical replicates were performed for each sample. Relative expression levels were calculated using the 2 − Co method (Livak and Schmittgen, 2001).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI GEO and GSE169056.

AUTHOR CONTRIBUTIONS
FM and RS conceived and designed the study. RS, WJ, XZ, and YL conducted the experiments. All authors carried out the data analysis. RS wrote the manuscript. FM, WJ, LY, ZZ, XZ, and QW revised the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was funded through an international cooperation between China and the EUC (2017YFE0111000).