Identification and Characterization of Salt-Responsive MicroRNAs in Vicia faba by High-Throughput Sequencing

Salt stress has detrimental effects on plant growth and development. MicroRNAs (miRNAs) are a class of noncoding RNAs that are involved in post-transcriptional gene expression regulation. In this study, small RNA sequencing was employed to identify the salt stress-responsive miRNAs of the salt-sensitive Hassawi-3 and the salt-tolerant ILB4347 genotypes of faba bean, growing under salt stress. A total of 527 miRNAs in Hassawi-3 plants, and 693 miRNAs in ILB4347 plants, were found to be differentially expressed. Additionally, 284 upregulated and 243 downregulated miRNAs in Hassawi-3, and 298 upregulated and 395 downregulated miRNAs in ILB4347 plants growing in control and stress conditions were recorded. Target prediction and annotation revealed that these miRNAs regulate specific salt-responsive genes, which primarily included genes encoding transcription factors and laccases, superoxide dismutase, plantacyanin, and F-box proteins. The salt-responsive miRNAs and their targets were functionally enriched by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, which showed that the miRNAs were involved in salt stress-related biological pathways, including the ABC transporter pathway, MAPK signaling pathway, plant hormone signal transduction, and the phosphatidylinositol signaling system, among others, suggesting that the miRNAs play an important role in the salt stress tolerance of the ILB4347 genotype. These results offer a novel understanding of the regulatory role of miRNAs in the salt response of the salt-tolerant ILB4347 and the salt-sensitive Hassawi-3 faba bean genotypes.


Introduction
Salt stress is one of the major abiotic stresses that negatively threaten plant growth and development. Worldwide, about 20% of agricultural lands and 50% of croplands are exposed to salt stress [1]. High salinity also results in secondary stresses, such as nutritional imbalance and oxidative stress, that lead to cellular damage, inhibition of growth, and a reduction in crop yield. Vicia faba is an economically important pulse. Owing to its high nutritional value, it is commonly used as food material and is widely distributed throughout the world. Understanding the role of V. faba microRNAs (miRNAs) in response to salt stress may help screen gene function and regulation in leguminous plants, and could contribute to more effective plant breeding.
The miRNAs comprise a class of endogenous non-coding RNAs that are approximately 21-24 nucleotides (nts) long, and are recognized as an important class of regulatory molecules,

Plant Materials, Growth Conditions, and Salt Stress Treatments
Two genotypes of V. faba, namely the salt-sensitive Hassawi3 and the salt-tolerant ILB4347, growing under control and stress conditions were used in this study. The plants from both the genotypes were washed in sterile deionized water and left to germinate in sand at 23 • C for 5 days. The seedlings that had germinated homogeneously were selected and planted in a mixture of sterilized sand and peat moss at a ratio of 3:1. The experiment was conducted in a greenhouse with a 16 h light/8 h dark cycle at a 26 • C/20 • C day/night temperature and 50-80% relative humidity at the College of Science, King Saud University, Riyadh, Saudi Arabia. The seedlings were irrigated with tap water and allowed to grow for 14 days. Hoagland nutrient solution (1/10 strength) was applied once a week. NaCl stress was initiated when the seedlings were nearly 21 days old and had attained two to three true leaves. Treatments constituted the control group and the groups that were administered NaCl at 150 mM concentrations. For both genotypes, the samples for RNA extraction were collected after two weeks of treatment, when the stressed plant started to wilt in comparison to the plants in the control setup and the plants that received 150 mM NaCl and had been exposed to maximum stress. The collected samples were quickly frozen in liquid nitrogen and stored at −80 • C until RNA extraction. A pool of three replicates was used for small RNA (sRNA) sequencing. HC-4 and HSA represented the Hassawi3 genotype growing under control and stress conditions, respectively. Similarly, IC-1 and IS-4 represented the ILB4347 genotype growing under control and stress conditions, respectively.

sRNA Library Construction and Sequencing
The total RNA was extracted from the foliar tissues of the salt-sensitive Hassawi-3 and the salt-tolerant ILB4347 genotypes of faba bean using the total RNA extraction kit (SV Total RNA isolation system, Promega, Madison, WI, USA), according to the manufacturer's instructions that involved the on-column digestion of any residual DNA in the samples of the control and stress-induced plants with RNase-free DNase I. The quality and quantity of the extracted RNA were measured using a NanoDrop ND-2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Then, the RNAs were sent to BGI (Shenzen, China) for sRNA library construction and high-throughput sequencing using BGISEQ-500 technology [29,30].

Data Filtering and Mapping Reads
The raw reads were first cleaned; this involved the removal of low-quality tags, tags without a 3 primer, tags with 5 primer contaminants, tags without the insert, tags with poly A, and tags that were shorter than 18 nts. The length distribution of the cleaned reads was next categorized for analyzing the composition of the sRNA data, and the 16-28-nt long sRNAs were selected for further analyses. The tags that remained after filtering comprised the "clean tags", which were stored in FASTQ format.
Then, the clean, high-quality sRNA tags were mapped to the reference genome using AASRA [31] and other sRNA databases, except Rfam, by using CMsearch [32].

Classification of sRNAs
Some sRNA tags can map to more than one category after alignment and annotation. In order to ensure that each unique sRNA maps to only one category after annotation, we followed the following priority rule: miRNA > piRNA > snoRNA > Rfam > other sRNAs.

Predictions of sRNAs
The typical hairpin structure of the miRNA precursor can be used to forecast novel miRNAs. We used RIPmiR [33] to predict novel miRNAs by analyzing their secondary structures, the Dicer cleavage site, and the minimum free energy of the unannotated sRNA tags, which could be mapped to the genome.
The piRNAs were predicted by Piano [34], which is based on an algorithm that uses piRNA-transposon interaction information, and the support vector machine (SVM) on these features. Small interfering RNAs (siRNAs) [35] are 22-24-nt-long double-stranded RNAs, in which each strand is 2-nt longer than the other at the 3 end. We aligned the tags on the basis of this structural feature to identify the siRNAs that met this criteria, as these tags could be potential siRNA candidates.

Analyzing sRNA Expression
The sRNA expression level of each gene was estimated by the Transcripts per Kilobase Million (TPM) method [36]. The TPM method can eliminate the influence of sequencing discrepancy while calculating sRNA expression. The gene expression thus calculated can therefore be directly used to compare the differences in gene expression among different samples. The following formula is used to calculate TPM:

Target Prediction
The computational prediction of miRNA targets is a critical step in the initial identification of miRNA:mRNA target interactions for experimental validation. In this study, psRobot [37] and TargetFinder [38] were used to identify the possible targets.

Screening the Differentially Expressed sRNAs (DESs)
With reference to "the significance of digital gene expression profiles" [39], a stringent algorithm was developed for identifying the differentially expressed genes between the two samples. If the number of tags from an sRNA is denoted as "x", and given the fact that the expression product of each RNA occupies only a small portion of the library, then the Poisson distribution can be computed from "x" by the following equation: If the total number of clean tags from sample 1 and sample 2 is N1 and N2, respectively, and if an sRNA "A" contains "x" tags in sample 1 and "y" tags in sample 2, then the probability of sRNA A being equally expressed between the two samples can be calculated as follows: We corrected the p-value corresponding to differential gene expression tests using the Bonferroni method [40]. As DESs analysis generates large multiplicity problems in which thousands of hypotheses are simultaneously tested, such as whether a gene "x" is differentially expressed between two groups, false positive (type I) and false negative (type II) errors are corrected by the false discovery rate (FDR) method [41]. Assuming that among "R" number of differentially expressed genes, "S" number of genes really show differential expression, and the remaining "V" genes are false positive, and supposing that the error ratio Q = V/R must remain below a certain cutoff, say 5%, then the value of FDR should be preset to a number no larger than 0.05. In this study, we preset the value of FDR to ≤0.001, and the absolute value of Log2Ratio was ≥1, which were set as default thresholds for judging the significance of the differences in gene expression. More stringent criteria with smaller values of FDR and higher fold-change values can be used to identify DESs.

GO Enrichment Analysis
Gene ontology (GO) enrichment analysis provides all the GO terms that are significantly enriched in a list of genes, and filters the genes with specific biological functions. In this method, all the genes are first mapped to the GO terms in the database (http://www.geneontology.org/), and then the gene numbers for each term are calculated. Lastly, a hypergeometric test is performed to identify the significantly enriched GO terms in the input gene list, based on GO: TermFinder' (https://www.yeastgenome.org/goTermFinder). The aforementioned algorithm was used for this analysis, and the P value was calculated by the following equation: where, N is the total number of genes with GO annotation, n is the number of DESs target genes in N, M is the total number of genes that are annotated by certain GO terms, and m is the number of DESs target genes in M. The calculated p-value was corrected by the Bonferroni Correction method [40], considering the corrected p-value ≤0.05 as a threshold. The GO terms that fulfill this condition are commonly defined as the significantly enriched GO terms. This analysis is performed to identify the main biological functions of the targets of the DESs.

Pathway Enrichment Analysis
Pathway-based analysis is helpful in further understanding the biological functions of the target genes of the DESs. KEGG [42] is an important public pathway-related database used for performing pathway enrichment analysis. Using this analysis, the significantly enriched metabolic pathways and signal transduction pathways of the DESs target genes are identified by comparing them with the whole genome as the background. The formulae used in pathway-based analysis are the same as those used in GO analysis, where "N" represents the total number of genes with KEGG annotation, "n" represents the total number of DESs target genes in N, "M" represents the total number of genes annotated by specific pathways, and 'm' depicts the number of DESs target genes in M.

sRNA Sequencing
In order to study the salt stress-responsive miRNAs in the salt-sensitive Hassawi-3 and salt-tolerant ILB4347 genotypes of faba bean, four sRNA libraries were constructed from the leaves. These libraries were then sequenced using BGISEQ-500 technology [29,30]. Table 1 briefly summarizes the sequencing information data derived from each sample. Upon sequencing the four libraries, a total of 197147801 raw reads were obtained; with the number of reads from the HC-4, HSA, IC-1, and IS-4 libraries being 51420111, 43734868, 48229628, and 53763194, respectively. The resulting raw sequence reads were computationally processed to remove the 3 adapter, which yielded a total of 174053893 clean reads from the four libraries; with the number of reads from the HC-4, HSA, IC-1, and IS-4 libraries being 47942551, 40682196, 43865674, and 41563472, respectively. The clean tags were used for analyzing the distribution of the 16-30 nt-long sRNAs. Generally, the length of sRNAs varies between 16-30 nt. Figure 1 shows the results of the length distribution analysis for the samples of sRNA. The length of sRNAs for HC-4 ranged between 17-34 nt, with the 19-30 nt-long RNAs being the most abundant. Similarly, for HSA, IC-1, and IS-4, the lengths of the sRNAs were 17-33 nt, 17-29 nt, and 17-31 nt, respectively. The sRNA length distribution (16-30 nt) of each library indicated that the species with a 21-24 nt length were the most abundant and diverse. The 24-nt-long RNAs were the most diverse among all the four sequencing libraries, followed by the 21-nt-long RNAs that were the second most abundant group of RNAs among the four libraries, with the exception of the HC-4 library, where the 21-nt-long RNAs were the most abundant.
After filtering, the clean tags were mapped to different sRNA databases, such as miRBase [43], Rfam [44], siRNA, piRNA, and snoRNA, among others. Table 1 enlists the individual mapping rates for each sample, while the distribution of the tags is depicted in Figure 1.  The sRNA reads were subsequently annotated by the RFam database (http://rfam.sanger.ac.uk/) and miRbase v21.0 (http://www.mirbase.org/) for classification. Table 2 shows the proportions of different types of sRNAs. The sRNAs were classified into different The sRNA reads were subsequently annotated by the RFam database (http://rfam.sanger.ac.uk/) and miRbase v21.0 (http://www.mirbase.org/) for classification. Table 2 shows the proportions of different types of sRNAs. The sRNAs were classified into different annotation categories of non-coding RNAs, including miRNA, tRNA, rRNA, snRNA, and snoRNA, among others. In order to ensure that each unique sRNA mapped to only one annotation category, the following priority rule was used: miRNA > piRNA > snoRNA > Rfam > other sRNA.

Annotation of sRNAs and miRNA Identification
A total of 1062 and 1156 miRNAs were identified from HC-4 and HSA, respectively; from which, 44 and 38 novel miRNAs were identified from HC-1 and HSA, respectively (Table 3). A total of 1513 and 1352 miRNAs were identified from IC-1 and IS-4, respectively; of which, 48 and 51 were novel miRNAs, identified from IC-1 and IS-4, respectively (Table 3). Interestingly, none of the siRNAs identified from the four samples were known siRNAs, with all of them being novel in nature.

Identification of Differentially Expressed miRNAs
DESs screening is primarily used to identify the sRNAs that are differentially expressed between different samples, followed by subsequent analysis. The DESs in each pairwise alignment are shown in Figure 2. Among the differentially expressed miRNAs, 284 were upregulated, while 243 miRNAs were downregulated in HC-4 and HSA combination ( Figure 2). Similarly, 298 of the differentially expressed miRNAs were upregulated in IC-1 and IS-4 combination, while 395 were found to be downregulated ( Figure 2).
The miRBase database v21.0 was used to identify the conserved miRNAs in our study from the perfect or near-perfect matches, with the latter having up to two mismatches. Based on sequence similarity, our analysis identified 492 known differentially expressed miRNAs in the salt-sensitive Hassawi-3 genotype (HC-4-vs.-HSA), belonging to 39 miRNA families (Table S1). Of these 39 miRNA families, seven families, namely those of miR166, miR167, miRNA396, miRNA159, miRNA171, miRNA168, and miR156, contained 130, 64, 52, 35, 33, 32, and 26, miRNAs, respectively; while 16 families contained two to 12 miRNAs, and 16 other families were represented by a single miRNA (Table S1). The majority of miRNAs in 25 of these families were found to have been upregulated (Table 4). A total of 35 sRNAs were identified as putative novel, faba bean miRNAs (Table S1 and  Table 5), of which 22 were found to be upregulated and 13 were downregulated (Table 5). The miRBase database v21.0 was used to identify the conserved miRNAs in our study from the perfect or near-perfect matches, with the latter having up to two mismatches. Based on sequence similarity, our analysis identified 492 known differentially expressed miRNAs in the salt-sensitive Hassawi-3 genotype (HC-4-vs.-HSA), belonging to 39 miRNA families (Table S1). Of these 39 miRNA families, seven families, namely those of miR166, miR167, miRNA396, miRNA159, miRNA171, miRNA168, and miR156, contained 130, 64, 52, 35, 33, 32, and 26, miRNAs, respectively; while 16 families contained two to 12 miRNAs, and 16 other families were represented by a single miRNA (Table S1). The majority of miRNAs in 25 of these families were found to have been upregulated (Table 4). A total of 35 sRNAs were identified as putative novel, faba bean miRNAs (Table S1 and Table 5), of which 22 were found to be upregulated and 13 were downregulated (Table 5).
The distribution miRNAs in these two genotypes showed that some of the miRNAs present in other legumes are also present in vicia faba or vice versa. The miRNAs such as 1507, 1508, 1509, 1512, 1514, 1521, 2086, 2109, 2199, miR2111, R2118, R5213 and R5232, 4414, 5213,5232, and 5234,5770 were generally found in legumes such as Cicer arietinum, Vigna unguiculata, Glycine max, Medicago truncatula, Glycine soja, Lotus japonicus, Phaseolus vulgaris, Lotus japonicus, and Acacia auriculiformis (www.miRBase.com). Of these miRNAs, 5213, 5232, 5234, and 21111, were found in either of our two faba bean genotypes (Table S2 and Table S2). However, the miR2111 was not only reported in legumes, but was also found in other plants such as in Populus trichocarpa, Vitis vinifera, Arabidopsis thaliana, and Malus domestica (www.miRBase.com). The reason for the presence of a limited number of legume-specific miRNAs in our study is that our small RNA library is not comprehensive owing to one tissue source under abiotic stress. Therefore, miRNAs in other tissues under biotic and abiotic stresses could speculate whether other miRNAs frequently present in legumes are also expressed in faba bean. Similarly, in the salt-tolerant genotype ILB4347 (IC-1-vs.-IS-4), 665 known miRNAs that belonged to 31 miRNA families were found to have been differentially expressed (Table S2). Of these 31 miRNA families, nine families, of miR166, miR167, miRNA396, miRNA160, miRNA171, miRNA164, miRNA159, miRNA168, and miR408, contained 162, 58, 57, 54, 54, 45, 36, 33, and 21 miRNAs, respectively; while 12 families contained two to twelve miRNAs, and 10 families were represented by a single miRNA (Table S2). The majority of miRNAs in 21 of these families were found to have been downregulated (Table S2). Additionally, 28 sRNAs were found to be novel miRNAs (Table 6), of which 11 were found to have been upregulated and 17 were downregulated ( Table 6).
The distribution miRNAs in these two genotypes showed that some of the miRNAs present in other legumes are also present in vicia faba or vice versa. The miRNAs such as 1507, 1508, 1509, 1512, 1514, 1521, 2086, 2109, 2199, miR2111, R2118, R5213 and R5232, 4414, 5213,5232, and 5234,5770 were generally found in legumes such as Cicer arietinum, Vigna unguiculata, Glycine max, Medicago truncatula, Glycine soja, Lotus japonicus, Phaseolus vulgaris, Lotus japonicus, and Acacia auriculiformis (www.miRBase.com). Of these miRNAs, 5213, 5232, 5234, and 21111, were found in either of our two faba bean genotypes (Table S2 and Table S2). However, the miR2111 was not only reported in legumes, but was also found in other plants such as in Populus trichocarpa, Vitis vinifera, Arabidopsis thaliana, and Malus domestica (www.miRBase.com). The reason for the presence of a limited number of legume-specific miRNAs in our study is that our small RNA library is not comprehensive owing to one tissue source under abiotic stress. Therefore, miRNAs in other tissues under biotic and abiotic stresses could speculate whether other miRNAs frequently present in legumes are also expressed in faba bean.

Target Prediction and Functional Analysis of the miRNAs
In the HC-4-vs.-HSA combination, a total of 4996 putative targets were predicted for 527 miRNAs; 4972 targets were of known miRNAs and 72 were of novel miRNAs (Table S3). For the IC-1 vs. IS-4 combination, a total of 5785 putative targets were predicted for 693 miRNAs; 4910 targets were of known miRNAs and 62 were of novel miRNAs (Table S4). Most of the predicted target genes encoded some stress-related transcription factors (TFs), including those of the auxin response factor (ARF) family, the AP2-like ethylene-responsive TF, squamosa promoter-binding proteins, myb domain proteins (MYBs), NAC domain-containing proteins, nuclear transcription factor Y, and bZIP proteins (Tables S3 and S4), as well as other proteins such as Argonaute (AGO2), glutamate decarboxylase, laccases, superoxide dismutase, plantacyanin, and F-box proteins.
For understanding the biological functions of miRNAs, all the putative target genes were subjected to GO functional classification by the Blast2GO software. GO-based enrichment analysis was carried out to further probe the possible role of miRNAs in response to salt stress. In the HC-4 vs. HSA combination, a total of 3708 potential miRNA targets were mapped to 1651 biological processes, 661 molecular functions, and 1396 cellular components (Figure 3). Some of the significant biological processes with the highest target numbers were the cellular process (348), metabolic process (318), biological regulation (148), regulation of the biological process (134), response to stimulus (102), and signaling (43). Among the categories for molecular function, binding (297), catalytic activity (270), transporter activity (23), nucleic acid binding TF activity (18), molecular transducer activity (8), and protein binding TF activity (6)  In the IC-1-vs.-IS-4 combination, a total of 3354 potential miRNA targets were mapped to 1488 biological processes, 574 molecular functions, and 1292 cellular components (Figure 4). Some of the significant biological processes with the highest target numbers were the cellular process (326), metabolic process (291), biological regulation (121), regulation of the biological process (102), response to stimulus (93), and signaling (34). Among the categories for molecular function, binding In the IC-1-vs.-IS-4 combination, a total of 3354 potential miRNA targets were mapped to 1488 biological processes, 574 molecular functions, and 1292 cellular components (Figure 4). Some of the significant biological processes with the highest target numbers were the cellular process (326), metabolic process (291), biological regulation (121), regulation of the biological process (102), response to stimulus (93), and signaling (34). Among the categories for molecular function, binding (253), catalytic activity (246), transporter activity (21), nucleic acid binding transcription factor activity (17), molecular transducer activity (4), and structural molecule activity (13) were the most abundant classes (Figure 4). The common cellular component terms were cell (278), cell part (278), organelle (232), and membrane (177).
Genes 2019, 10 13  For the HC vs. HSA combination, pathway analysis was performed by KEGG annotation, which returned 93 KEGG pathways that were classified into five main categories ( Figure 5). The complete set of 93 pathways for the HC vs. HSA combination after KEGG annotation has been summarized in supplementary Table S5. The annotation results showed that most of the genes were involved in "metabolic function" (327 unigenes), followed by those involved in "genetic information processing" (153), "environmental information processing" (67), "cellular processes" (52), and "organismal systems" (12). The only four pathways that were categorized under the "environmental information processing" category were from plants (24, 3.69%). The majority of representative pathways for the unigenes under the "metabolism" category were involved in the biosynthesis of secondary metabolites (ko01110; 45, 6.92%). The two pathways categorized under the "organismal systems" (environmental adaptation) category were involved in plant-pathogen interactions (ko04626; 10, 1.54%) and plant circadian rhythms (ko04712; 2, 0.31%). For the HC vs. HSA combination, pathway analysis was performed by KEGG annotation, which returned 93 KEGG pathways that were classified into five main categories ( Figure 5). The complete set of 93 pathways for the HC vs. HSA combination after KEGG annotation has been summarized in Supplementary Table S5. The annotation results showed that most of the genes were involved in "metabolic function" (327 unigenes), followed by those involved in "genetic information processing" (153), "environmental information processing" (67), "cellular processes" (52), and "organismal systems" (12). The only four pathways that were categorized under the "environmental information processing" category were from plants (24, 3.69%). The majority of representative pathways for the unigenes under the "metabolism" category were involved in the biosynthesis of secondary metabolites (ko01110; 45, 6.92%). The two pathways categorized under the "organismal systems" (environmental adaptation) category were involved in plant-pathogen interactions (ko04626; 10, 1.54%) and plant circadian rhythms (ko04712; 2, 0.31%). Genes 2019, 10 14  For pathway analysis, KEGG annotation was performed for the IC-1 vs. IS-4 combination, which returned 89 KEGG pathways that were classified into five main categories ( Figure 6). The complete set of 89 pathways identified for the IC-1 vs. IS-4 combination is summarized in supplementary Table S6. The majority of genes were annotated under the "metabolic function" (243) category, followed by the "genetic information processing" (123) category, "environmental information processing" (54), "cellular processes" (33), and "organismal systems" (9) categories. The only four pathways under the "environmental information processing" category were categorized under the plant hormone signal transduction (ko04075; 34, 6.53%), plant MAPK signaling pathway (ko04016; 6, 1.15%), ABC transporter (ko02010; 3, 0.58%), and phosphatidylinositol signaling system (ko04070; 13, 2.5%) categories. The majority of representative pathways for the unigenes under the "metabolism category" were involved in metabolic pathways (ko01100; 76, 14.59%), and the biosynthesis of secondary metabolites (ko01110; 30, 5.76%). The two pathways categorized under the "organismal systems" (environmental adaptation) category were involved in plant-pathogen interactions (ko04626; 9, 1.73%). For pathway analysis, KEGG annotation was performed for the IC-1 vs. IS-4 combination, which returned 89 KEGG pathways that were classified into five main categories ( Figure 6). The complete set of 89 pathways identified for the IC-1 vs. IS-4 combination is summarized in Supplementary Table S6. The majority of genes were annotated under the "metabolic function" (243) category, followed by the "genetic information processing" (123) category, "environmental information processing" (54), "cellular processes" (33), and "organismal systems" (9) categories. The only four pathways under the "environmental information processing" category were categorized under the plant hormone signal transduction (ko04075; 34, 6.53%), plant MAPK signaling pathway (ko04016; 6, 1.15%), ABC transporter (ko02010; 3, 0.58%), and phosphatidylinositol signaling system (ko04070; 13, 2.5%) categories. The majority of representative pathways for the unigenes under the "metabolism category" were involved in metabolic pathways (ko01100; 76, 14.59%), and the biosynthesis of secondary metabolites (ko01110; 30, 5.76%). The two pathways categorized under the "organismal systems" (environmental adaptation) category were involved in plant-pathogen interactions (ko04626; 9, 1.73%). Genes 2019, 10 15

Discussion
Salinity adversely affects plant growth and development. Therefore, in order to cope with salt stress, plants adapt themselves by modulating various stress-responsive genes at the transcriptional and post-transcriptional levels. In recent years, the role of miRNAs in the regulation of gene expression has been increasingly understood. The miRNAs are said to have pivotal roles in plant responses to abiotic and biotic stresses [45]. Numerous studies have demonstrated that miRNA-mediated gene regulation plays an important role in the salt stress response of several plant species, including Arabidopsis [16], soybean [11], barley [46], and sugarcane [47].
In this study, we investigated millions of sRNA sequences from faba bean for understanding the effects of salt stress on plant growth and development under miRNA regulation. Our results revealed that salt stress modulated the expression of miRNAs and their predicted targets in faba bean.
More than 40 million sRNA reads were produced by high-throughput sequencing from each library of faba bean genotype growing under control and salt stress conditions. The sRNAs of faba bean comprised different classes of RNAs, including snRNAs, tRNAs, snoRNAs, rRNAs, and miRNAs. The majority of these were unannotated in the existing sRNA libraries owing to the scarcity of information on this species. The annotation performed herein was in accordance with numerous previous studies [48,49]. The results indicated that to date, numerous sRNAs remain to be identified. Among the different types of sRNAs, miRNAs are considered as being crucial in the regulation of gene expression at the post-transcriptional level. In plants, the genes under miRNA regulation are involved in vegetative and reproductive growth, as well as in the plant stress response [49,50]. Analysis of the four sRNA libraries of faba bean revealed that most of the sRNAs were 21-nt and 24-nt-long (Figure 1). This pattern of length distribution has been observed in several other plant

Discussion
Salinity adversely affects plant growth and development. Therefore, in order to cope with salt stress, plants adapt themselves by modulating various stress-responsive genes at the transcriptional and post-transcriptional levels. In recent years, the role of miRNAs in the regulation of gene expression has been increasingly understood. The miRNAs are said to have pivotal roles in plant responses to abiotic and biotic stresses [45]. Numerous studies have demonstrated that miRNA-mediated gene regulation plays an important role in the salt stress response of several plant species, including Arabidopsis [16], soybean [11], barley [46], and sugarcane [47].
In this study, we investigated millions of sRNA sequences from faba bean for understanding the effects of salt stress on plant growth and development under miRNA regulation. Our results revealed that salt stress modulated the expression of miRNAs and their predicted targets in faba bean.
More than 40 million sRNA reads were produced by high-throughput sequencing from each library of faba bean genotype growing under control and salt stress conditions. The sRNAs of faba bean comprised different classes of RNAs, including snRNAs, tRNAs, snoRNAs, rRNAs, and miRNAs. The majority of these were unannotated in the existing sRNA libraries owing to the scarcity of information on this species. The annotation performed herein was in accordance with numerous previous studies [48,49]. The results indicated that to date, numerous sRNAs remain to be identified. Among the different types of sRNAs, miRNAs are considered as being crucial in the regulation of gene expression at the post-transcriptional level. In plants, the genes under miRNA regulation are involved in vegetative and reproductive growth, as well as in the plant stress response [49,50]. Analysis of the four sRNA libraries of faba bean revealed that most of the sRNAs were 21-nt and 24-nt-long ( Figure 1). This pattern of length distribution has been observed in several other plant species, including mulberry [49], Populus euphratica [51], potato [52], and barley [53]. The 24-nt-long miRNA class was the most diverse among all the four libraries, except for HC-4, where the 21-nt-long miRNAs were the most abundant. Other studies have indicated that the length of miRNAs primarily varies between 21 nt and 22 nt, while siRNAs are mostly 24-nt-long [4,54]. Our study therefore revealed that the more enriched faba bean sRNAs could in reality be miRNAs that had been cleaved by DCL1. The importance of the miRNA length distribution lies in the fact that it allows easy alignment with the RNA-induced silencing complex (RISC) for regulating gene expression by inhibiting its translation or by degrading the target mRNAs, depending on the miRNA-mRNA similarities.

Salt Stress-Responsive miRNAs and Their Targets in Faba Bean
Recently, miRNAs have evolved as a valuable genetic tool for interpreting stress tolerance at the molecular level and for finally regulating the stress response in crops [55]. The identification of salt stress-responsive miRNAs and their subsequent functional analysis will aid the understanding of the stress-responsive mechanism in plants. Earlier studies in maize [27], wheat [56], rice [57], barley [46], and sugarcane [47] have demonstrated that miRNAs such as miR156, miR169, miR160, miR159, miR168, miR171, miR172, miR393, and miR396 are the major salt stress-responsive miRNAs in plants. In our study, we identified 527 differentially expressed miRNAs in the HC-4 vs. HSA combination; 284 of which were upregulated and 243 were downregulated (Table S1). A total of 693 differentially expressed miRNAs were identified in the IC-1 vs. IS-4 combination; 298 of which were upregulated and 395 were downregulated (Table S2). The aforementioned salt stress-responsive miRNAs were also identified in our study, which suggests the presence of common salt stress-related miRNAs. Additionally, other conserved miRNAs and a few novel miRNAs were identified as salt stress-responsive sRNAs in both genotypes (Tables 5 and 6). Therefore, these results serve as novel information on the salt stress-responsive miRNAs of plants. It is well-established that miRNA-mediated complex regulatory networks are involved in the regulation of gene expression at the post-transcriptional level in plants [58]. The GO terms of the putative targets revealed that the targets had stimulus signaling, binding, catalytic, transporter, and nucleic acid binding TF activities (Figures 3 and 4). KEGG analysis revealed that some targets of the salt stress-responsive miRNAs in the two faba bean genotypes considered in this study were mapped to salt stress-related pathways, such as plant hormone signal transduction, flavonoid biosynthesis, ABC transporter activity, ubiquitin-mediated proteolysis, and DNA repair ( Figures 5  and 6, Tables S5 and S6) [59].
The large number of upregulated and downregulated miRNAs induce a more pronounced change in the expression of multiple downstream genes. The response to salt stress in crops is accompanied by a wide range of intracellular processes, including signal perception, signal transduction, transcription, and protein biosynthesis, among others. Although different species of plants respond to stress by employing numerous miRNA-mediated regulatory strategies [60], some studies have reported that hub miRNAs, such as miR169, miR171, miR393, miR396, and miR398, among others, are linked to several abiotic stresses, such as drought, salinity [61,62], and cold [63]. These studies report that the miRNA targets are involved in sensing stress, signal transduction, and other stimuli. Previous studies have demonstrated that miR171 and miR393 are upregulated under conditions of salt stress in A. sp., wheat, and barley [46,64,65]. In this study, we observed that the expression of miRNAs in both the miR393 and miR171 families was upregulated and downregulated under conditions of salt stress in all the genotypes of faba bean, except for the miR393 family, in which the miRNA expression was only upregulated for the HC-4 vs. HSA combination. These results suggest the existence of common regulatory mechanisms for salt stress response in plants, and these miRNAs might control the same targets in different plants [66]. The miR393 family targets genes encoding the F-box protein family that includes TIR1 and AFB2 in A. sp. and rice, and inhibits the growth of lateral roots under abscisic acid (ABA) treatment or osmotic stress [65,67].
In this study, miR166, miR171, miR398, miR396, and miR1432 were identified in both genotypes of faba bean. The miR171 family is known to target the myeloblastosis (MYB) family of TFs that might have a role in the regulation of osmotic balance under drought and salt stress conditions. Alptekin and coworkers (2017) reported that both salt stress and drought affect the osmotic balance of plant cells [62]. In this study, miR396 was found to be differentially expressed in both the genotypes of faba bean. Using expressed sequence tags (ESTs), Kantar and coworkers (2011) reported that the miR396 family targets the growth factor-like (GRL) TF and its putative heat-shock protein, as the changes in their expression were consistent with the changes in the expression of miR396. The study demonstrated that the expression of GRL TF and its putative heat-shock protein was upregulated following the downregulation of miR396 under stress [68]. The function of heat-shock proteins is to protect other proteins from degradation under stress. The downregulation of miR396 and the subsequent regulation of its targets would therefore improve the tolerance of faba bean to salt stress. A better understanding of miRNAs and their functions under stress conditions would improve the efficacy and reliability of the application of miRNA-mediated gene regulation for increasing plant stress tolerance.

Conclusions
To the best of our knowledge, this is the first study to identify the salt stress-responsive miRNAs in Vicia faba using high-throughput sequencing technology. Four sRNA libraries were generated and sequenced from two faba bean genotypes under control and salt stress conditions. A total of 1062 and 1156 miRNAs were identified in the HC-4 and HSA samples, respectively; while a total of 1513 and 1352 miRNAs were identified in IC-1 and IS-4, respectively. Differential expression analysis revealed that 527 miRNAs were differentially expressed in the HC-4 vs HSA combination; whereas 693 miRNAs were differentially expressed in the IC-1 vs. IS-4 combination. GO enrichment analysis revealed that the targets of these differentially expressed salt stress-responsive miRNAs were highly enriched in the corresponding GO terms, and possessed stimulus signaling, binding, catalytic, transporter, and nucleic acid binding TF activities, among others. KEGG analysis suggested that several of the miRNAs in faba bean participate in stress-related pathways, including plant hormone signal transduction, the MAPK signaling pathway, flavonoid biosynthesis, ubiquitin-mediated proteolysis, apoptosis, and ABC transporter activity. This study enhanced the existing genetic information and resources of salt-responsive miRNAs in faba bean. This will not only improve our understanding of the roles of miRNAs in post-transcriptional gene regulation during salt stress response, but will also aid studies on the miRNA-based genetic improvement of salt tolerance in cultivated faba bean genotypes.