Global identification and analysis of microRNAs involved in salt stress responses in two alfalfa (Medicago sativa ‘Millennium’) lines

Abstract: Alfalfa is an important economic crop; a mutant (M) strain was identified during planting and production. M plants consistently had better relative water content and relative electrical conductivity under higher salt conditions compared with the wild type (WT) plants, suggesting that M plants have higher tolerance for salt. To understand the microRNAs (miRNAs) involved in salt stress response in alfalfa, 128 miRNAs were identified from the WT and M alfalfa plants under normal and saline conditions. Of the 128 miRNAs, 29 and 23 differentially expressed miRNAs were identified in the M vs. WT control (M-CK vs. WT-CK) and salt-stressed M vs. WT (M-salt vs. WT-salt) comparison, respectively. These miRNAs responded to salt stress and showed different expression patterns after salt treatment. Their potential target genes were predicted and further analysed by GO classification and KEGG pathway analysis, where the majority of target genes were associated with plant growth and development, and exhibited significant changes in WT and M plants. In addition, compared with the WT plants, miR172-CNGC, miR319-CAX2, miR408-NHX and miR2590-CHX14/15 showed significant upregulation in M alfalfa plants, suggesting that M plants have higher ion transport levels. The differential expression profiles of miRNAs and putative target genes were further validated by quantitative real-time polymerase chain reaction. It is speculated that these miRNAs are involved in the increased salt tolerance of the M alfalfa plants.


Introduction
High soil salinity is a widespread and serious problem facing irrigated lands throughout the world; soil salinity can lead to lower yields and limited capacity to grow crops. Salinity is one of the most significant abiotic factors that plants are exposed to. Salt tolerance of plants is a very complex phenomenon involving many biological processes and pathways which operate at the wholeplant level (Deinlein et al. 2014). When plants are subjected to salt stress, an early response involving osmotic stress and the accumulation of Na + and Cl − occurs (Zhu 2001), which inhibits plant growth by affecting metabolic processes and photosynthetic efficiency. Plants experiencing salt exhibit phenotypes of drought, wilting, and even death (Hasegawa 2013). Molecular mechanisms of salt tolerance in plants have been studied in many plant species including model plants, crops, and halophytes (Bartels and Sunkar 2005). Various genes are involved in the response to salt stress that regulates the ion balance to facilitate the intra-and inter-cellular Na + homeostasis of plants. The plasma membrane Na + /H + antiporter SOS1 is responsible for mediating Na + efflux and controls long-distance Na + transport in plants, and SOS2 and SOS3 regulate SOS1 transport activity in Arabidopsis (Qiu et al. 2002). In addition, one particularly important class of exchangers is the NHX-type, which regulates vacuolar Na + /H + exchange, thereby maintaining ion balance (Jiang et al. 2010).
Small RNAs (sRNAs) are important regulators of growth and development, and of biotic and abiotic stress responses of plants (Khraiwesh et al. 2012;Dong et al. 2016). Several classes of sRNAs, including microRNAs (miRNAs) and small interfering RNAs (siRNAs), regulate stress response by altering the expression of target genes at the post-transcriptional level (Leung and Sharp 2010). miRNAs are 20-24 nt short noncoding RNAs that are produced from a precursor hairpin structure, and they downregulate genes by targeting mRNAs for cleavage or translational repression (Bartel 2004). Many miRNAs are associated with the salt response. In Arabidopsis, 12 miRNAs (miR156, miR158, miR159, miR165, miR167, miR168, miR169, miR171, miR319, miR393, miR394 and miR396) are upregulated in response to high salt stress, while no significantly downregulated miRNAs are known to date (Liu et al. 2008). In Populus euphratica Olivier, a total of 211 known miRNAs and 162 new miRNAs were identified under salt treatment (Li et al. 2013). In Vigna unguiculata (L.) Walp. (cowpea), 18 conserved miRNAs belonging to 16 miRNA families are differentially expressed under salt stress using a comparative genomic approach (Paul et al. 2011). Furthermore, Jia et al. (2009) identified miR398 in response to abscisic acid (ABA) and salt stress, and showed its differential and dynamic regulation in Populus tremula L. and Arabidopsis thaliana (L.) Heynh. (Jia et al. 2009). miR169 members (miR169g, miR169n and miR169o) are induced by high-salinity treatment, and transiently inhibit the expression of NF-YA transcription factor (Zhao et al. 2009). Taken together, miRNAs are widely involved in the salt stress response in plants.
Alfalfa is an important economic crop that is used as animal forage. However, salinity significantly limits alfalfa production and quality. Although miRNAs regulate tolerance to salt stress, this has not been widely studied in alfalfa. Here, a mutant of alfalfa with high salt tolerance was identified during planting and production. The aim of this study was to identify miRNAs involved in salt stress response in alfalfa, and their downstream targets, using high-throughput sequencing technology.

Plant materials and growth conditions
Wild type (WT) and mutant (M) alfalfa (Medicago sativa L. 'Millennium') were kept in plastic pots and cultured in a greenhouse under 16 h light/8 h dark (25°C/20°C) cycles with a relative humidity of 70%-75%. All materials were obtained from Zhejiang Agriculture and Forest University (Zhejiang, People's Republic of China). After 3 mo of culturing, at least six uniformly sized individual plants were treated with a 250 mmol L −1 NaCl solution, and at least six individuals were treated with water as control. Leaves were harvested at 0 and 72 h of the treatment. All the samples were collected from three randomly selected plants, immediately frozen in liquid nitrogen, and stored at −80°C.

Determination of physiological parameters
A total of 1 g of fresh leaves (fresh weight, FW) were collected from the different salt treatment stages (0,24,36,48,72, and 96 h after treatment). Relative water content (RWC) was measured according to Yamasaki and Dillenburg (1999). All samples were floated in distilled water for about 24 h to measure the turgid weight (TW). Then the leaves were dried for 30 min at 105°C, and the dried leaves were used to measure dry weight (DW). The RWC was calculated as RWC (%) = [(FW − DW)/(TW − DW)] × 100. For the relative electrical conductivity (REC), fresh leaves (0.5 g) were obtained from different salt treatment stages (0,24,36,48,72, and 96 h after treatment). The REC was measured using a conductivity meter (DDS-11A) following methods described by Li et al. (2000). Three biological replicates were measured for each parameter.

RNA extraction and sRNA sequencing
Total RNA was isolated from 0.5 g of leaves using the RNAiso Plus reagent (TaKaRa, Dalian, People's Republic of China) according to the manufacturer's protocol. RNase-free DNase I (TaKaRa) was used to remove residual DNA. Equal amounts of leaf RNA from the three independent biological replicates were pooled for sRNA library construction. The 0 and 72 h samples from the WT and M plants were sequenced using the Illumina HiSeq™ 2000 platform. The datasets are available in the NCBI repository (http://trace.ncbi.nlm.nih.gov/Traces/sra_sub/ sub.cgi?) under the accession number SRP150034. Adapter sequences were removed from the raw sequencing reads, and 10% of the nitrogen sequences and low-quality (Q < 5) reads were removed to obtain highquality sequences. In addition, the length distributions of clean reads, total reads, and unique sequences were summarized.

miRNA annotation and target prediction
The set of sRNAs (18-30 nt) from the libraries were filtered to remove ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nucleolar RNAs (snoRNAs) and small nuclear RNAs (snRNAs). The remaining sequences were compared with equivalent plant noncoding sRNA sequences deposited in either GenBank (http://www. ncbi.nlm.nih.gov/genbank/) or the Rfam (http://www. sanger.ac.uk/science/tools/rfam) database. The filtered sRNA sequences were aligned in the miRNA database (http://www.mirbase.org) to identify known miRNAs using BLASTN with a maximum of two mismatches. The secondary structures of all known miRNA precursors were obtained using the Mfold version 3.5 software. In addition, the novel miRNAs were predicted using the MIREAP software (http://sourceforge.net/projects/ mireap/). miRNAs that had a |log2(fold change)| ≥ +1 and a P value < 0.05 were considered significantly differentially expressed. In addition, a reads per minute (RPM) value of zero was assigned a nominal RPM of 0.01. Targets for differentially expressed miRNAs were predicted using the alfalfa genome (ftp://ftp.ncbi.nlm.nih.gov/ genomes/all/GCF/000/219/495/GCF_000219495.3_MedtrA17_ 4.0/GCF_000219495.3_MedtrA17_4.0_genomic.fna.gz). Target prediction was conducted following methods from Allen et al. (2005). The predicted target genes were functionally classified using the Gene Ontology (GO) database (http://www.geneontology.org) and assigned to a pathway through the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/kegg1.html).

Real-time quantitative polymerase chain reaction (RT-qPCR) assay
Total RNA samples were reverse transcribed using a PrimeScript miRNA qPCR starter kit version 2.0 (Takara) following the manufacturer's protocol. Each RT-qPCR reaction was 20 μL RT-qPCR and contained 10 μL SYBR R Premix Ex Taq TM II (Takara), 5 μL cDNA (100 ng μL −1 ), 1.0 μL 10 μM forward primer, and 1.0 μL UnimiR qPCR Primer. The primers of target genes were designed using Primer version 5.0, and the sequences are shown in Supplementary Tables S1 and S2 1 . The PCR programme was as follows: 95°C for 2 min; 40 cycles of 95°C for 15 s and 55°C for 15 s; and 72°C for 20 s. The gene Elongation Factor 1-alpha (EF-1α) (GenBank accession number: JZ818475) was used as the reference gene (Castonguay et al. 2015). Three biological replicates were tested for each sample, and relative expression levels were calculated using the 2 −ΔΔCT method (Livak and Schmittgen 2001).

M plants have higher salt resistance than WT plants
After salt treatment WT plants began to wilt and turn yellow, and some plants died (Fig. 1). In contrast, M plants exhibited higher salt tolerance (Fig. 1). Next, the physiological parameters of the plants, including RWC and REC, were measured after salt treatment. RWC rapidly decreased in WT plants relative to M plants after salt treatment ( Fig. 2A). REC values were lower in M plants when compared with WT plants at any length of salt treatment time (Fig. 2B). These results suggested that M plants have higher salt tolerance than WT plants.
The length distribution of the sRNAs is shown in Fig. 3. Reads with a length of 21 or 24 nt were abundant in the four libraries, followed by those with a length of 22 nt or 23 nt (Fig. 3). When these libraries were compared with the NCBI GenBank and miRbase databases, many of the sRNAs were matched to different types of sRNAs, including miRNA, siRNA, rRNA, tRNA, snRNA, snoRNA, and repeats. Noncoding sRNAs were removed except for miRNAs for subsequent analysis ( Table 2).

Analysis of differentially expressed miRNAs from WT and M plants
A total of 128 miRNAs were obtained from four libraries (Supplementary Table S3 1 ). To explore the expression patterns of miRNAs, differentially expressed miRNAs were identified using the following criteria: |log2(fold change)| ≥ 1 and P value < 0.05. The number of differentially expressed miRNAs is shown in Fig. 4. Twenty-nine miRNAs were differentially expressed between the WT-salt vs. WT-CK comparison, and 16 downregulated and 13 upregulated miRNAs were identified under salt stress (Supplementary Table S4 1 ).
Similarly, 28 differentially expressed miRNAs were found between the M-salt vs. M-CK. Of these miRNAs, 13 were downregulated, and 15 were upregulated in the M-salt library (Supplementary Table S5 1 ). Furthermore, a total of 29 differentially expressed miRNAs were obtained from the M-CK vs. WT-CK comparison, including 17 downregulated and 12 upregulated miRNAs (Table 3). Meanwhile, 23 differentially expressed miRNAs were identified between the M-salt vs. WT-salt comparison, including 10 miRNAs that were upregulated    Note: log2FC, log 2-fold change; FDR, false discovery rate; WT-CK, control wild type; M-CK, control mutant. and 13 miRNAs that were downregulated in the M plants ( Table 4).

Identification of target genes associated with differentially expressed miRNAs
To further analyse the response to salt stress in the M and WT plants, target genes of the miRNAs that were identified in the previous analysis were characterized using the four libraries through NR (non-redundant), Swiss-Prot, GO, KEGG, Clusters of Orthologous Groups (COG), EuKaryotic Orthologous Groups (KOG) and Pfam databases. A total of 1833, 874, 1698 and 491 potential target genes were predicted for the WT-salt vs. Table S6 1 ). The GO functional categories of the target genes were analysed, including cellular components, molecular functions, and biological processes for WT and M alfalfa under salt treatment (Fig. 5). In the cellular components, GO terms associated with cell, cell part, and organelle were the most abundant in the four comparisons. For molecular functions, the most abundant GO terms were binding, catalytic activity, and transport activity. Metabolic process and single-organism process were the two most dominant GO categories in the biological processes. In the KEGG database, these potential genes were assigned to 84 pathways ( Supplementary  Fig. S1 1 ). The results of COG and KOG were further analysed and shown in Supplementary Figs. S2 and S3 1 .

WT-CK, M-salt vs. M-CK, M-CK vs. WT-CK and M-salt vs. WT-salt comparisons, respectively (Supplementary
Salt stress impacts plant growth, development, and ion transport. Therefore, the targets of differentially expressed miRNAs involved in growth, development, and ion transport were identified in the WT and M alfalfa (Table 5).

RT-qPCR validation of miRNAs and target genes
Validation of the differential expression of miRNAs and their target genes that were identified using RNA-Seq analysis was performed using RT-qPCR. Ten differentially expressed miRNAs were identified in the WT and M plants (Fig. 6). These miRNAs included miR156, miR166, miR168, miR169, miR319, miR390, miR395, miR396, miR398 and miR408. The other differentially expressed miRNAs are shown in Supplementary  Fig. S4 1 . Their target genes that were associated with ion transport, growth, and development were analysed including SPL2 (MTR_8g080680) Discussion miRNAs are an important class of endogenous sRNAs that are widely distributed in plants, where they regulate growth, development, and environmental stress responses by silencing the mRNAs of target genes. Here, the leaf materials from two alfalfa lines with or without salt stress were used to identify the miRNAs involved in plant response to salt stress.

Bioinformatic analysis of miRNAs under salt stress in WT and M alfalfa
Based on the length distribution of sRNAs, 18-30 nt sRNAs were dominant, and 21 and 24 nt sRNAs were abundant in both WT and M libraries, followed by 22 and 23 nt (Fig. 3). This is consistent with previous studies in other species, such as Triticum compactum Host (Liu and Sun 2017), rice (Sunkar et al. 2008), and potato (Zhang et al. 2013). In addition, four libraries of alfalfa exposed to salt stress were sequenced, and a total of 128 miRNAs were identified in the WT and M libraries (Supplementary  Table S3 1 ). Among these libraries, most miRNAs were differentially expressed in WT and M plants after salt stress (Fig. 4). These results indicated that salt stress might regulate miRNA expression, and that these miRNAs directly or indirectly affect plant growth and development to reduce salt damage (Frazier et al. 2011). Additionally, the differentially expressed miRNAs in M plants are likely involved in the M plant's ability to adapt to saline conditions.

Analysis of miRNA responses to salt stress in M alfalfa
The M alfalfa plants were more resistant to salinity than the WT plants ( Figs. 1 and 2). The expressions of miRNAs associated with salt stress were compared between the M-CK vs. WT-CK and M-salt vs. WT-salt libraries. In total,   (Tables 3 and 4). Compared with WT-CK, miR156, miR159, miR160, miR166, miR167, miR319, miR396, miR398 and miR408 were downregulated in the M plants. A similar result was observed in the WT-salt vs. M-salt comparison, where miR156, miR166, miR169, miR319, miR396, miR398, miR408 and miR894 were also downregulated. These homologous miRNAs are involved in environmental stress responses (Kantar et al. 2011), and they contribute to salinity stress response in many plant species (Sunkar et al. 2012;Xie et al. 2014). They are necessary and sufficient for development and response to environmental stress in Arabidopsis (Stief et al. 2014), and are significantly downregulated during the early stages of salt stress (Jatan et al. 2018). The expression of miR1507, miR1509, miR1510 and miR2590 becomes strongly downregulated after stress treatment in plants (Jiang et al. 2014;Khandal et al. 2017). miR2661 regulates the plant-pathogen interaction and may regulate resistance against Phytophthora infestans (Mont.) de Bary in tomatoes (Luan et al. 2015). In addition, miR172 and miR482 were more upregulated in M-CK than in WT-CK plants (Table 3), but their expression was similar in M-salt vs. WT-salt. After salt treatment, miR160, miR168, miR390, miR395, and miR397 were upregulated in the M plants (Table 4). Overall, these miRNAs were differentially expressed in WT and M plants under normal and salt conditions, and they are likely involved in salt stress response, as well as growth and development of alfalfa, through translational repression or mRNA decay.

Targets genes associated with salt response in M alfalfa
Usually miRNAs downregulate or silence their target genes, thereby responding to environmental stresses and modulating plant growth and development (Chen et al. 2010). A series of target genes, such as SPL, GRAS, ARF, GRF, TCP, WRKY, MYB and bHLH (Table 5), were associated with the regulation of salt stress response and the growth and development in the alfalfa plants that were included in this study. Under salt stress these genes were upregulated in M plants compared with WT plants, except for SPL genes. miR156 post-transcriptionally regulate SPL genes, and 35S:MIM156 can increase the sensitivity of stress tolerance in Arabidopsis (Cui et al. 2014). The expression of the SPL2 gene (MTR_8g080680) was significantly decreased in WT and M plants after salt stress (Fig. 7), implying that salt stress negatively affected the growth and development of M and WT plants. The putative target of miR166, GRAS family transcription factor GAI (MTR_3g065980), was upregulated in M-salt compared with WT-salt (Fig. 7). This transcription factor is involved in salt and drought stress response (Golldack et al. 2014). In addition, other potential targets such as miR168-AGO (MTR_4g094858), miR396-GRF5 (MTR_4g125490) and miR397-PYL8 (MTR_3g090980), also showed differential expression patterns in the WT and M groups after salt stress (Fig. 7). AGO protein allows RNA silencing to carry out its biological functions (Vaucheret 2008). GRF protein binds to the promoter of dehydration-responsive element binding protein2a (DREB2A) and represses its expression to increase osmotic and heat tolerance (Kim et al. 2012). The 35S:PYL8/RCAR3 transgenic plants show ABA-resistant drought response and strongly inhibit early root growth (Saavedra et al. 2010). In alfalfa, PYL8 (MTR_3g090980) expression was strongly induced after salt stress, demonstrating that PYL8 may upregulate ABA signalling during salt stress in WT and M plants. Furthermore, we found that the target genes miR172-CNGC (MTR_7g012260), miR319-CAX2 (MTR_2g105640), miR408-NHX (MTR_7g099800) and miR2590-CHX14/CHX15 (MTR_5g074360, MTR_5g074380) showed significant differential expression in the M plants after salt treatment, and these target genes were involved in ion transport. Cyclic nucleotide-gated ion channel (CNGC) is an important component of Ca 2+ signal transduction which regulates various biological processes, such as development and response to salt stress (Saand et al. 2015). CAX2 is an exchanger that maintains intracellular Ca 2+ ion MTR_6g016250 TOM20 plant-specific import receptor subunit TOM20 MTR_2g062920 -ATP-citrate synthase alpha chain protein 2-like miR6214 MTR_8g064470 -Ubiquitin-like-specific protease 1D miR8736 MTR_2g084020 -Peroxidase homeostasis (Kumar et al. 2013). NHX and CHX14/15 are involved in Na + transport, and NHX is a salt tolerance determinant that affects cellular ion homeostasis (Deinlein et al. 2014). CHX14/15 were strongly induced by salt stress, and they are involved in K + accumulation and homeostasis (Cellier et al. 2004). Taken together, these differently expressed miRNAs play important roles in regulating salt stress in M plants by downregulating target genes, and they likely contribute towards the higher salt tolerance observed in M plants.

Conclusion
The M alfalfa plants are more resistant to high salinity than WT plants. Using sRNA-seq results, changes in miRNA expression patterns in M and WT plants under normal and high salt conditions were studied. Differentially expressed miRNAs and their target genes were identified, and these genes not only responded to salt stress but also contributed to plant growth and development. Furthermore, CNGC (miR172), CAX2 (miR319), NHX (miR408) and CHX14/CHX15 (miR2590) were identified as genes involved in ion exchange and transport that were differentially expressed in M compared with WT. Therefore, miRNAs may regulate these genes to provide salt resistance displayed by M alfalfa plants.