Abstract
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.
Keywords: microRNAs, Hassawi-3, ILB4347, faba bean, salt stress
1. 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
[[34]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, involved in
post-transcriptional gene regulation mediated via mRNA degradation or
the repression of mRNA translation [[35]2,[36]3]. These RNAs are
largely derived from intergenic regions and are produced from
single-stranded primary miRNAs with unique hairpin structures [[37]4].
These were originally discovered in Caenorhabditis elegans [[38]5].
However, their widespread presence has been reported in animals
[[39]6], plants [[40]7], and certain viruses [[41]8]. It is generally
known that miRNAs play important regulatory roles in several biological
processes, including those during the developmental, metabolism,
pathogen defense, and stress responses in plants [[42]9,[43]10].
Recently, it was discovered that miRNAs respond to various abiotic
stresses in plants, including drought [[44]11,[45]12,[46]13,[47]14],
high salinity [[48]15,[49]16,[50]17], low temperatures [[51]16,[52]18],
oxidative stress [[53]19], hypoxic stress [[54]20,[55]21], mechanical
stress [[56]22], and UV-B radiation [[57]15], as well biotic stresses
[[58]23].
Multiple gene expression mechanisms have evolved in plants in order to
address high-salinity-induced stress. These mechanisms relate to a
broad spectrum of biochemical, cellular, and physiological processes,
including signal transduction, energy metabolism, transcription,
protein biosynthesis, membrane trafficking, and photosynthesis
[[59]24]. The transcriptional regulation of several miRNAs and genes in
response to high salt stress has been extensively studied
[[60]25,[61]26]. These studies suggested that plant responses to salt
stress could be shaped by miRNA-guided gene regulation. Microarray
analysis revealed that several miRNAs, such as miR156, miR159, miR167,
miR168, miR171, miR319, and miR396, showed differential expression
during the salt stress response of Arabidopsis sp. [[62]16] and Zea
mays [[63]27]. Wide-ranging sequencing data has been produced from next
generation sequencing (NGS) technology for the detection of
salt-responsive miRNAs in various plant species. Employing this
technology, 104 differentially expressed miRNAs were identified in the
salt-stressed functional soybean nodules [[64]28].
In this study, we employed high-throughput sequencing technology to
identify the conserved and novel miRNAs in the salt-tolerant and
salt-sensitive cultivars of V. faba growing under conditions of high
salt stress. Using advanced bioinformatics analysis, the changes in
miRNA expression following salt treatment were studied in comparison to
those of the control. To study the potential underlying mechanism of
the miRNA-mediated gene expression regulation in faba bean under salt
stress, the miRNA and target interaction networks were further analyzed
by Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and
Genomes (KEGG).
2. Materials and Methods
2.1. 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.
2.2. 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 [[65]29,[66]30].
2.3. 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 [[67]31] and other sRNA databases, except Rfam, by
using CMsearch [[68]32].
2.4. 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.
2.5. Predictions of sRNAs
The typical hairpin structure of the miRNA precursor can be used to
forecast novel miRNAs. We used RIPmiR [[69]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 [[70]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) [[71]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.
2.6. Analyzing sRNA Expression
The sRNA expression level of each gene was estimated by the Transcripts
per Kilobase Million (TPM) method [[72]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:
[MATH:
TPM=C<
/mi>×106N :MATH]
2.7. 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 [[73]37] and
TargetFinder [[74]38] were used to identify the possible targets.
2.8. Screening the Differentially Expressed sRNAs (DESs)
With reference to “the significance of digital gene expression
profiles” [[75]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:
[MATH:
P(x)=<
/mo>e−λλxx! :MATH]
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:
[MATH: 2∑i+0<
/mn>i=y
P(i|x) :MATH]
[MATH: or :MATH]
[MATH:
2 × (1−∑i+0<
/mn>i=y
P(i|x) (if
∑i+0<
/mn>i=y
P(i|x)>0.5) :MATH]
[MATH: P(y |x)=(N2N1)y(
x+y)!<
/mo>x!y!(1+N2N1
)x | y |1 :MATH]
We corrected the p-value corresponding to differential gene expression
tests using the Bonferroni method [[76]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 [[77]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.
2.9. 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
([78]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’
([79]https://www.yeastgenome.org/goTermFinder). The aforementioned
algorithm was used for this analysis, and the P value was calculated by
the following equation:
[MATH: Ρ=1−∑i=0<
/mn>m−1(M
i)
(N−M<
/mrow>n−i
)(Nn) :MATH]
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 [[80]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.
2.10. Pathway Enrichment Analysis
Pathway-based analysis is helpful in further understanding the
biological functions of the target genes of the DESs. KEGG [[81]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.
3. Results
3.1. 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
[[82]29,[83]30]. [84]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. [85]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.
Table 1.
Summary of the results of high-throughput sequencing of the sRNAs from
the leaves of two V. faba genotypes growing under salt stress
conditions and the control. The alignment statistics of the tags
following alignment to the reference genome are summarized herein.
Sample Name Sequence Type Raw Tag Count Clean Tag Count Percentage (%)
Number of Mapped Tags Percentage (%)
HC-4 SE50 51420111 47942551 93.24 36387118 75.9
HSA SE50 43734868 40682196 93.02 26272525 64.58
IC-1 SE50 48229628 43865674 90.95 21451106 48.9
IS-4 SE50 53763194 41563472 77.31 20568122 49.49
[86]Open in a new tab
Figure 1.
[87]Figure 1
[88]Open in a new tab
Sequence length distribution of the faba bean sRNAs. The size
distribution of the sRNA libraries. The y-axis represents the
sequences, while the x-axis depicts the 16–36 nt-long sequences for
each of the four sequenced libraries. Length distribution of (A) HC-4,
(B) HSA, (C) IC-1 and (D) IS-4 libraries.
After filtering, the clean tags were mapped to different sRNA
databases, such as miRBase [[89]43], Rfam [[90]44], siRNA, piRNA, and
snoRNA, among others. [91]Table 1 enlists the individual mapping rates
for each sample, while the distribution of the tags is depicted in
[92]Figure 1.
The sRNA reads were subsequently annotated by the RFam database
([93]http://rfam.sanger.ac.uk/) and miRbase v21.0
([94]http://www.mirbase.org/) for classification. [95]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.
Table 2.
Abundance of the reads of the different sRNA categories from the leaf
samples of the two faba bean genotypes growing under control and
salinity stress conditions.
Type HC-4 HSA IC-1 IS-1
Count (%) Count (%) Count (%) Count (%)
Total 47942551 100 40682196 100 43865674 100 41563472 100
Intergenic 9985569 20.83 10292988 25.3 9475423 21.6 8197558 19.72
Mature (miRNA) 3224519 6.73 1653404 4.06 6841472 15.6 5633410 13.55
Rfam other sncRNA 73325 0.15 94167 0.23 24135 0.06 35172 0.08
snRNA 1405 0 4831 0.01 2381 0.01 959 0
unmap 11040957 23.03 13725210 33.74 21005681 47.89 20270917 48.77
rRNA 1773606 3.7 1100198 2.7 244021 0.56 358242 0.86
Hairpin 38 0 9 0 80 0 147 0
snoRNA 30566 0.06 14138 0.03 9971 0.02 6058 0.01
Precursor 27019 0.06 11116 0.03 24531 0.06 40935 0.1
Repeat 21771765 45.41 13683732 33.64 6237390 14.22 7018371 16.89
tRNA 13782 0.03 102403 0.25 589 0 1703 0
[96]Open in a new tab
3.2. 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 ([97]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
([98]Table 3). Interestingly, none of the siRNAs identified from the
four samples were known siRNAs, with all of them being novel in nature.
Table 3.
Summary of the small non-coding RNAs detected from each sample.
Sample Known miRNA Count Novel miRNA Count Known piRNA Count Novel
piRNA Count Known siRNA Count Novel siRNA Count
IS-4 1301 51 0 0 0 2578
HSA 1118 38 0 0 0 1879
HC-4 1018 44 0 0 0 1919
IC-1 1465 48 0 0 0 2241
[99]Open in a new tab
3.3. 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
[100]Figure 2. Among the differentially expressed miRNAs, 284 were
upregulated, while 243 miRNAs were downregulated in HC-4 and HSA
combination ([101]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 ([102]Figure 2).
Figure 2.
[103]Figure 2
[104]Open in a new tab
Graphical representation of the number of differentially expressed
miRNAs in salt-sensitive (HC-4 vs HSA) and salt-tolerant genotypes
(IC-1 vs IS-4) growing under control and salt stress conditions.
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 ([105]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 ([106]Table S1). The
majority of miRNAs in 25 of these families were found to have been
upregulated ([107]Table 4). A total of 35 sRNAs were identified as
putative novel, faba bean miRNAs ([108]Table S1 and [109]Table 5), of
which 22 were found to be upregulated and 13 were downregulated
([110]Table 5).
Table 4.
Differentially expressed conserved miRNA families.
Conserved miRNA Family DESs in HC-4-vs.-HSA (Sensitive Genotype) DESs
in IC-1-vs.-IS-4 (Tolerant Genotype) Target Gene Family
Number of Upregulated Members Number of Downregulated Members Number of
Upregulated Members Number of Downregulated Members
miRNA156 20 6 11 7 Squamosa promoter-binding proteins
miRNA157 2 1 2 1 Squamosa promoter-binding proteins
miRNA159 27 8 8 28 MYB transcription factors/TCP transcription factors
miRNA160 4 4 11 43 Auxin Response factors
miRNA162 9 3 8 8 Dicer Like protein/E3 ubiquitin-protein ligase
RNF144A-like isoform X1
miRNA164 3 2 12 33 NAC domain protein/NAC transcription factor-like
protein
miRNA165 0 1 2 2 Homeo domain-Zip transcription
factors/homeobox-leucine zipper protein ATHB-14-like/bZIP transcription
factor
miRNA166 39 91 82 80 Homeobox-leucine zipper protein ATHB-15-like
isoform X2/bZIP transcription factor
miRNA167 38 26 27 31 Transmembrane protein, putative/translation
initiation factor eIF-2B delta subunit
miRNA168 19 13 24 9 Argonautes/protamine P1 family protein
miRNA169 3 1 2 4 HAP2/NFY transcription factors/CCAAT-binding
transcription factor
miRNA171 19 14 12 42 Scarecrow-like transcription factors/GRAS family
transcription regulator
miRNA172 1 0 0 1 AP2 domain transcription factors/AP2-like
ethylene-responsive transcription factor/myb-like transcription factor
family protein
miRNA319 6 3 6 12 MYB transcription factors/TCP transcription factors
miRNA390 3 8 4 5 TAS3-primary transcripts/LRR receptor-like kinase
family protein
miRNA391 1 0 0 1 TAS3-primary transcripts/zinc finger CCCH domain
protein
miRNA393 8 1 2 2 F-Box protein/transport inhibitor response 1
protein/Ubiquitin
miRNA394 2 0 2 2 F-Box protein/Zinc finger CCCH domain-containing
protein ZFN-like
miRNA396 36 16 31 26 Growth regulating factors/F-box protein
interaction domain protein /BZIP transcription factor bZIP80
miRNA397 2 2 4 4 Laccases
miRNA398 2 7 3 13 Cu/Zn superoxide dismutases (CSD)/BAG family
molecular chaperone regulator-like protein
miRNA399 3 3 8 4 Phosphate transporter/ubiquitin-conjugating enzyme
E2/OBP3-responsive protein
miRNA408 1 6 16 5 Plantacyanins/uclacyanin-2-like/basic blue-like
protein
miRNA2111 2 5 6 2 DNA replication factor CDT1-like
protein/calcineurin-like phosphoesterase, family protein
miRNA482 0 2 1 1
[111]Open in a new tab
DESs-Differentially Expressed sRNAs.
Table 5.
Comparative expression profiles of the novel miRNA genes in
HC-4-vs.-HSA.
miRNA ID Count (HC-4) Count (HSA) TPM (HC-4) TPM (HSA) log2 Ratio
(HSA/HC-4) Regulation Profile (Up/Down) (HSA/HC-4) p-Value FDR
novel_mir1 193 869 5.12 31.31 2.612408 Up 1.81E-153 2.31 × 10^−152
novel_mir5 24 124 0.64 4.47 2.804131 Up 3.67E-25 2.05 × 10^−24
novel_mir6 0 48 0.001 1.73 10.75656 Up 1.12E-18 5.58 × 10^−18
novel_mir7 227 1120 6.02 40.35 2.744733 Up 1.78E-208 2.63 × 10^−207
novel_mir9 0 82 0.001 2.95 11.5265 Up 2.42E-31 1.47 × 10^−30
novel_mir10 14 894 0.37 32.21 6.44384 Up 8.00E-307 1.47 × 10^−305
novel_mir13 33 115 0.88 4.14 2.234055 Up 2.47E-18 1.21 × 10^−17
novel_mir16 17 597 0.45 21.51 5.578939 Up 1.75E-194 2.51 × 10^−193
novel_mir17 491 794 13.03 28.61 1.134682 Up 2.99E-44 2.03 × 10^−43
novel_mir19 39 81 1.03 2.92 1.503324 Up 3.24E-08 1.14 × 10^−7
novel_mir22 11 32 0.29 1.15 1.987509 Up 2.35E-05 7.14 × 10^−5
novel_mir29 82 180 2.18 6.49 1.57389 Up 9.97E-18 4.82 × 10^−17
novel_mir37 492 817 13.06 29.43 1.172133 Up 6.79E-48 4.77 × 10^−47
novel_mir38 6941 74963 184.18 2700.76 3.874177 Up 0 0
novel_mir39 164 1013 4.35 36.5 3.068809 Up 2.05E-212 3.08 × 10^−211
novel_mir40 31 86 0.82 3.1 1.918572 Up 1.06E-11 4.35 × 10^−11
novel_mir41 114 2010 3.02 72.42 4.583768 Up 0 0
novel_mir42 259 1808 6.87 65.14 3.245162 Up 0 0
novel_mir43 50 130 1.33 4.68 1.815082 Up 6.45E-16 3 × 10^−15
novel_mir45 736 2609 19.53 94 2.266969 Up 0 0
novel_mir46 143 220 3.79 7.93 1.065123 Up 3.38E-12 1.41 × 10^−11
novel_mir48 919 1754 24.39 63.19 1.373407 Up 9.10E-129 1.04 × 10^−127
novel_mir8 50 0 1.33 0.001 −10.3772 Down 1.20E-12 5.13 × 10^−12
novel_mir11 926 75 24.57 2.7 −3.18587 Down 4.68E-136 5.50 × 10^−135
novel_mir15 35 0 0.93 0.001 −9.86109 Down 4.71E-09 1.73 × 10^−8
novel_mir20 35 0 0.93 0.001 −9.86109 Down 4.71E-09 1.74 × 10^−8
novel_mir23 794 252 21.07 9.08 −1.21443 Down 3.39E-35 2.17 × 10^−34
novel_mir25 59 0 1.57 0.001 −10.6165 Down 8.33E-15 3.76 × 10^−14
novel_mir27 2970 65 78.81 2.34 −5.0738 Down 0 0
novel_mir28 16 0 0.42 0.001 −8.71425 Down 0.000168 0.000477
novel_mir33 709 173 18.81 6.23 −1.5942 Down 1.26E-46 8.65 × 10^−46
novel_mir34 479 22 12.71 0.79 −4.00797 Down 1.78E-85 1.64 × 10^−84
novel_mir36 41 0 1.09 0.001 −10.0901 Down 1.72E-10 6.84 × 10^−10
novel_mir44 34 0 0.9 0.001 −9.81378 Down 8.17E-09 2.97 × 10^−8
novel_mir53 230 47 6.1 1.69 −1.85179 Down 1.71E-19 8.64 × 10^−19
[112]Open in a new tab
TPM: Transcripts per Kilobase Million.
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 ([113]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 ([114]Table S2). The majority of miRNAs in 21 of these
families were found to have been downregulated ([115]Table S2).
Additionally, 28 sRNAs were found to be novel miRNAs ([116]Table 6), of
which 11 were found to have been upregulated and 17 were downregulated
([117]Table 6).
Table 6.
Comparative expression profiles of novel miRNA genes in IC-1-vs.-IS-4.
miRNA ID Count (IC-1) Count (IS-4) TPM (IC-1) TPM (IS-4) log2 Ratio
(IS-4/IC-1) Regulation Profile (up/down) (IC-1 vs IS-4) p-Value FDR
novel_mir6 44 256 1.87 11.62 2.6354999 Up 7.82× 10^−41 4.42 × 10^−40
novel_mir8 338 6850 14.38 310.9 4.434315 Up 0 0
novel_mir9 0 207 0.001 9.4 13.198445 Up 5.11 × 10^−66 3.49 × 10^−65
novel_mir18 0 188 0.001 8.53 13.05833 Up 5.02 × 10^−60 3.27 × 10^−59
novel_mir19 200 1538 8.51 69.81 3.0362027 Up 3.94 × 10^−275 5.07 ×
10^−274
novel_mir22 0 587 0.001 26.64 14.701306 Up 7.53 × 10^−186 7.83 ×
10^−185
novel_mir30 0 200 0.001 9.08 13.148477 Up 8.25 × 10^−64 5.54 × 10^−63
novel_mir36 163 496 6.93 22.51 1.6996388 Up 5.04 × 10^−45 2.99 × 10^−44
novel_mir39 134 453 5.7 20.56 1.8508064 Up 2.78 × 10^−46 1.67 × 10^−45
novel_mir42 192 1707 8.17 77.48 3.245416 Up 0 0
novel_mir5 347 663 14.76 30.09 1.0275914 Up 2.31 × 10^−28 1.10 × 10^−27
novel_mir13 1404 298 59.72 13.53 −2.1420523 Down 2.13 × 10^−156 1.98 ×
10^−155
novel_mir14 62 10 2.64 0.45 −2.552541 Down 8.44 × 10^−10 2.72 × 10^−9
novel_mir16 116 42 4.93 1.91 −1.368015 Down 2.42 × 10^−8 7.48 × 10^−8
novel_mir20 254 60 10.8 2.72 −1.9893528 Down 4.41 × 10^−27 2.07 ×
10^−26
novel_mir27 10812 1387 459.9 62.95 −2.8690419 Down 0 0
novel_mir29 1387 564 59 25.6 −1.2045711 Down 1.14 × 10^−68 7.94 ×
10^−68
novel_mir31 245 75 10.42 3.4 −1.6157486 Down 4.83 × 10^−20 1.98 ×
10^−19
novel_mir32 4098 1680 174.31 76.25 −1.1928461 Down 3.81 × 10^−196 4.01
× 10^−195
novel_mir33 1370 157 58.27 7.13 −3.0307793 Down 2.35 × 10^−225 2.65 ×
10^−224
novel_mir4 165 11 7.02 0.5 −3.811471 Down 1.52 × 10^−34 7.93 × 10^−34
novel_mir40 451 39 19.18 1.77 −3.4377815 Down 1.63 × 10^−84 1.22 ×
10^−83
novel_mir44 240 0 10.21 0.001 −13.317695 Down 1.24 × 10^−69 8.68 ×
10^−69
novel_mir46 759 263 32.28 11.94 −1.4348377 Down 1.01 × 10^−49 6.25 ×
10^−49
novel_mir48 6461 707 274.83 32.09 −3.0983438 Down 0 0
novel_mir50 5351 1112 227.61 50.47 −2.173066 Down 0 0
novel_mir52 450 116 19.14 5.26 −1.8634561 Down 5.42 × 10^−43 3.15 ×
10^−42
novel_mir53 871 180 37.05 8.17 −2.1810656 Down 6.75 × 10^−100 5.38 ×
10^−99
[118]Open in a new tab
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 ([119]www.miRBase.com). Of these miRNAs, 5213, 5232,
5234, and 21111, were found in either of our two faba bean genotypes
([120]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 ([121]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.
3.4. 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 ([122]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 ([123]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 ([124]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 ([125]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)
were the most abundant classes ([126]Figure 3). The common cellular
component terms were cell (305), cell part (304), organelle (235), and
membrane (201).
Figure 3.
[127]Figure 3
[128]Open in a new tab
GO classifications in the HC-4 vs. HSA combination.
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 ([129]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 ([130]Figure 4).
The common cellular component terms were cell (278), cell part (278),
organelle (232), and membrane (177).
Figure 4.
[131]Figure 4
[132]Open in a new tab
GO classifications in the IC-1 vs. IS-4 combination.
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 ([133]Figure 5). The complete set of 93 pathways
for the HC vs. HSA combination after KEGG annotation has been
summarized in [134]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%).
Figure 5.
[135]Figure 5
[136]Open in a new tab
KEGG pathway classifications of the HC vs. HSA combination.
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 ([137]Figure 6). The complete set of 89
pathways identified for the IC-1 vs. IS-4 combination is summarized in
[138]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%).
Figure 6.
[139]Figure 6
[140]Open in a new tab
KEGG pathway classifications for the IC-1 vs. IS-4 combination.
4. 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 [[141]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
[[142]16], soybean [[143]11], barley [[144]46], and sugarcane
[[145]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 [[146]48,[147]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 [[148]49,[149]50].
Analysis of the four sRNA libraries of faba bean revealed that most of
the sRNAs were 21-nt and 24-nt-long ([150]Figure 1). This pattern of
length distribution has been observed in several other plant species,
including mulberry [[151]49], Populus euphratica [[152]51], potato
[[153]52], and barley [[154]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 [[155]4,[156]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 [[157]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 [[158]27], wheat [[159]56], rice
[[160]57], barley [[161]46], and sugarcane [[162]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 ([163]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
([164]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 ([165]Table 5 and [166]Table 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 [[167]58]. The GO terms of
the putative targets revealed that the targets had stimulus signaling,
binding, catalytic, transporter, and nucleic acid binding TF activities
([168]Figure 3 and [169]Figure 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 ([170]Figure 5 and [171]Figure 6, [172]Tables S5 and S6)
[[173]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 [[174]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 [[175]61,[176]62], and cold [[177]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 [[178]46,[179]64,[180]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
[[181]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 [[182]65,[183]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 [[184]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 [[185]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.
5. 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.
Supplementary Materials
The following are available online at
[186]https://www.mdpi.com/2073-4425/10/4/303/s1, Table S1:
Differentially expressed miRNAs in Vicia faba (HC-4-vs-HSA) after
salinity treatment. Table S2: Differentially expressed miRNAs in Vicia
faba (IC-1-vs-IS-4) after salinity treatment. Table S3: Target
prediction and annotation of differentailly expressed miRNAs in HC-4 vs
HSA comparison. Table S4: Target prediction and annotation of
differentailly expressed miRNAs in IC-1 vs IS-4 comparison. Table S5:
List of pathways of differentially expresed genes in HC-4-vs-HSA
comparisons. Table S6: List of Pathway of differentially expressed
genes in IC-1-vs-IS-4 comparisons.
[187]Click here for additional data file.^ (583.8KB, zip)
Author Contributions
S.M.A., A.A.A., and H.M.M. did most of the experimental work. M.A.K.,
S.M.A., and S.S.A. assisted in analyzing the results. M.A.K. and S.M.A.
drafted the manuscript. I.A.A., H.M.M., S.S.A., and A.A.A. revised the
manuscript.
Funding
We are thankful to the Deanship of Scientific Research of King Saud
University, Riyadh, KSA, for the financial assistance (Grant No.
RG-1438-036).
Conflicts of Interest
The authors declare no conflict of interest.
References