Increase in soil salt causes osmotic and ionic stress to plants, which inhibits their growth and productivity. Rice production is also hampered by salinity and the effect of salt is most severe at the seedling and reproductive stages. Salainity tolerance is a quantitative property controlled by multiple genes coding for signaling molecules, ion transporters, metabolic enzymes and transcription regulators. MicroRNAs are key modulators of gene-expression that act at the post-transcriptional level by translation repression or transcript cleavage. They also play an important role in regulating plant’s response to salt-stress. In this work we adopted the approach of comparative and integrated data-mining to understand the miRNA-mediated regulation of salt-stress in rice. We profiled and compared the miRNA regulations using natural varieties and transgenic lines with contrasting behaviors in response to salt-stress. The information obtained from sRNAseq, RNAseq and degradome datasets was integrated to identify the salt-deregulated miRNAs, their targets and the associated metabolic pathways. The analysis revealed the modulation of many biological pathways, which are involved in salt-tolerance and play an important role in plant phenotype and physiology. The end modifications of the miRNAs were also studied in our analysis and isomiRs having a dynamic role in salt-tolerance mechanism were identified.
Plants have the capability to perceive and respond to multiple environmental cues by the coordinated regulation of gene-expression. The genetic networks are regulated in a spatio-temporal manner at transcriptional, post-transcriptional and translation levels, for maintaining the cellular homeostasis , . Extreme fluctuations in the environmental conditions such as drought (water deficit), salinity, nutrient deficiency, frost and heat are perceived as abiotic stresses, which limit the plant productivity . In recent times, the disequilibrium in water cycle and poor drainage facilities in irrigated land has caused the problem of increase in soil salt . The accumulation of excessive levels of salt has affected around 20 % of global irrigated land .
Soil salinity affects the plant’s ability to absorb water and also leads to high cellular concentrations of salts within the plant cells . Increased uptake of Na+ and Cl− ions disrupts the osmotic and ionic balance leading to inhibition of many biochemical, molecular and physiological processes. Plant’s response to salinity stress involves activation of complex physiological traits and metabolic pathways controlled by multiple genes . During the last couple of decades, number of genes involved in the salt-stress tolerance pathway have been identified. These include signal transducers, ion transporters , , , , transcription factors ,  and metabolic enzymes .
The Glyoxalase (Gly) enzymes have been identified as important players in salt-stress tolerance as they can detoxify the stress induced accumulation of the cytotoxic aldehyde methylglyoxal (MG). In a coordinated action, GlyI catalyzes the initial step of the pathway  by converting MG to SD-lactoylglutathione (SLG) in presence of GSH while GlyII catalyzes the next step, in which the cytotoxic SLG is hydrolyzed to D-lactate and GSH is released. This pathway has been reported to confer salt-stress tolerance in plants by removing the cellular toxicity and regulating the GSH homeostasis .
The molecular reprogramming to combat salt-stress also involves the recruitment of several miRNAs. They are 21–24 nucleotide molecules processed endogenously in the nucleus from single stranded RNA by the action of DICER-like (DCL) proteins. Then mature miRNAs are transported to the cytoplasm where they associate with the ARGONAUTE (AGO) containing RNA Induced Silencing Complex (RISC) to regulate gene expression. The miRNA loaded RISC recognizes the target transcripts in a sequence-dependent manner to bring about transcript cleavage or translation repression . miRNAs are recognized as crucial effectors of diverse biological processes in plants , , , .
Deep sequencing is a powerful tool for prediction and identification of miRNAs and their target mRNAs. This approach relies on complex computational algorithms that also assist in understanding their biological functions. We have employed algorithms based on the established criteria for identifying the miRNAs and determining their digital expression patterns in leaf tissues of different rice varieties grown under normal and salt-stressed conditions. Rice is among the most salt-sensitive crops and its sensitivity varies at different developmental stages. The emergence and early seedling stages are highly sensitive to salt-stress , . However, rice plants are adapted to grow in a wide environmental range. In particular, Pokkali is relatively tolerant to salinity during early germination stages, whereas Pusa Basmati is a salt-sensitive rice variety.
In this study the miRNA profiles of these two contrasting varieties were obtained and compared with the profiles observed for the salt-tolerant Gly over-expressing transgenic rice lines , , , , , . Targets of the miRNAs were predicted using degradome data and their expression pattern was examined using the RNAseq data. Gene enrichment and KEGG pathway analysis was performed to identify the associated cellular and metabolic pathways , . It revealed the modulation of many biological pathways, which are involved in salt-tolerance and play an important role in plant phenotype as well as physiological traits. The 5′ and 3′ end modifications of the miRNAs were also studied and isomiRs having a dynamic role in salt-tolerance mechanism were identified . The changes in the expression patterns of the miRNAs and their alteration in response to salt-stress were followed across rice varieties and transgenic lines varying in their sensitivity to salt-stress. The analysis revealed that adaptation to salt-stress involves the stringent regulation of biochemical and cellular pathways.
2.1 Data Collection
2.1.1 Plant Materials and Growth Conditions
Mature seeds of rice Pusa Basmati (wild type), Gly overexpressing Pusa Basmati (transgenic lines) and Pokkali were used in this study. The seeds were surface sterilized with 10 % commercial bleach for 5 min, washed thoroughly with sterile water and placed on germinating sheets. The seeds were grown under controlled conditions, temperature (28.2 °C), relative air humidity (70 %), and 16/8 h light/dark cycle. For further analysis, leaf tissues were harvested from normal or 200 mM NaCl stressed 15-days-old seedlings.
2.1.2 Small RNA Library Preparation
1 gram tissue was used for total RNA isolation as described earlier . The small RNA (sRNA) was enriched by LiCl precipitation and used for library construction. For each tissue three biological replicates were used for the analysis. The libraries thus made were used for deep sequencing on GAII sequencer (Illumina).
2.1.3 Library Preparation for RNAseq and Degradome Sequencing
1 gram tissue was used for total RNA isolation as described earlier . Paired-end sequencing of transcriptome were performed on Illumina HiSeq 2000 platform as per sequencing guidelines.
Degradome sequencing was done to identify the miRNA mediated cleaved product of mRNA. The construction of libraries were as per established protocols . Reverse transcription and PCR libraries were sequenced on GAII sequencer (Illumina) using the 5′ adapter only, resulting in the sequencing of the first 36 nucleotides of the inserts that represented the 5′ ends of the mRNAs.
2.2.1 Analysis of NGS Data for MiRNA Identification
Each sRNA library yielded 36 nt long sequences with base quality scores and read counts. NGS tool kit was used for quality filtering from raw sRNA reads . The reads qualifying the quality control parameter were processed for removing the adaptor from 3′ end. The adaptor trimmed reads from sRNA library were size filtered to select for 17–27 nt reads for further analysis. These were aligned to miRbase release 21.0 ,  (http://www.mirbase.org/) using bowtie alignment tool  with zero mismatch to identify the miRNA. Only perfectly aligned mature miRNAs and their precursors were chosen for further analysis. Since the cDNA library was used for sequencing, the resultant sequences were having thymine (T) rather than uracil (U) in the individual reads.
Digital expression status of the miRNA was generated by normalizing the expression value as TPM (transcript per million).
N represents the fold-change
2.2.2 Identification of End Modifications and isomiRs of Known miRNAs
MiRMOD, a miRNA modification tool was used for the identification of end modified miRNAs . The 17–27 nt reads were aligned to the precursor sequence with perfect match in reference to mature miRNA. For identification of isomiRs the reads mapping outside to the identified miRNA region and having template independent 5′ and 3′ modification were shortlisted. The abundant sequence with their corresponding miRNA was considered as putative isomiRs.
2.2.3 Target Prediction
The miRNA target prediction was done using CleaveLand 4.4.3 tool . The degradome and transcriptome data was also used for target prediction. All predicted targets with p-value ≤0.05 were considered for further analysis. Target prediction of isomiRs was done by psRNATarget tool with default parameter  (http://plantgrn.noble.org/psRNATarget/). To study the full gene annotation of target mRNA, and their associated pathways GO enrichment was performed using online available database AgriGO (http://bioinfo.cau.edu.cn/agriGO/). The KEGG pathway analysis was performed using direct search at Rice Oligonucleotide Array Database (ROAD)  (http://www.ricearray.org/).
2.2.4 Analysis of NGS Data for mRNA Identification
The expression pattern of the targets of the mRNAs was determined using the RNAseq data. The transcriptome data contained maximum 101 nt long reads. Quality filtered transcriptome reads were mapped against TIGR genome cDNA (OSA1R5) using burrow wheeler alignment software to identify the mRNA transcripts . Aligned transcripts were used for expression profiling in different datasets under various experimental conditions. The overall computation methodology is shown in Figure 1.
Digital expression status of the transcripts was generated by normalizing the expression value as RPKM (Reads Per Kilobase Million).
where, C = Number of reads mapped to a gene
N = Total mapped reads in the experiment
L = Exon length in base-pairs for a gene
2.3 Data Integration
Three types of data sets were used for this analysis including the sRNAseq, RNAseq and degradome data. The sRNAseq was used for miRNA identification, isomiR prediction and identification of 5′/3′ end modifications. Whereas the degradome and RNAseq data sets were used for target prediction and obtaining their digital expression profiles, respectively. This information was useful in identifying the miRNA regulated targets and comparing the expression patterns of the miRNAs to those of the target sequences. The analysis was extended to identify the miRNA controlled pathways that are deregulated under salt stress (Figure 2).
3.1 Analysis of Small RNA Sequences and MiRNA
The deep sequencing data from six small RNA (sRNA) libraries representing leaf tissues isolated from Pusa Basmati wild type (PB), Gly over-expressing transgenics in Pusa Basmati background (GL) and Pokkali (PK) seedlings grown under normal (NL) and salt-stressed (SL) conditions were analyzed. A total of 443 mature miRNAs, represented in at least one library, were identified on aligning with miRBase release 21. The preliminary analysis revealed that the identified mature miRNAs range from 20 to 24 nt in length with the 21 nt species being dominantly abundant in each sample followed by 24 nt species (Figure 3). It was seen that salt-stress resulted in an increase in the total number of miRNAs in PK whereas it caused a decline in the number of miRNAs in PB and GL. However, the percentage decrease was higher in the wild type plants as compared to the Gly-transgenics.
The expression of identified miRNAs varied across the rice varieties under different conditions (Table 1). It was observed that 186 miRNAs were found to express diversely in all libraries under control (normal) conditions while 184 miRNAs were abundant in >2 libraries. The expression of specific miRNAs was restricted to particular library only (Figure 4) with 9, 12 and 17 miRNAs expressing in PBSL, GLSL and PKSL libraries, respectively. The details of the miRNA expression reads are available in Supplementary File-1.
|Abbreviation used||Tissue||Variety||Growth condition||Known miRNAs||Total no. of reads|
|GLNL||Leaf||Gly overexpressing Pusa Basmati||Normal||346||2,63,71,105|
|GLSL||Leaf||Gly overexpressing Pusa Basmati||Salt-stress||324||2,67,80,311|
To understand if the variation in miRNA expression patterns is related to the salt-tolerant/susceptible physiology further comparative studies were performed. The expression levels of miRNAs were compared between the NL and SL tissues of respective rice varieties and their fold-change was plotted (Supplementary File-1). It was observed that a large pool of miRNA was down-regulated in the salt-sensitive PBSL (Figure 5A). The level of down-regulation ranged up to 6-fold or more. In the salt-tolerant PKSL tissues a relatively similar level of down-regulation pattern was observed though for a smaller pool of miRNA. In contrast, small number of miRNAs was up-regulated in the salt-sensitive PBSL (Figure 5A). The level of up-regulation ranged up to 4-fold. In the salt-tolerant PKSL tissues a higher volume of up-regulation was observed with the variation levels ranging beyond 6-fold. Whereas the salt-tolerant GL exhibited a less variant pool of deregulated miRNAs. The order of miRNA up-regulation was higher though the level of fluctuation was mainly within the 2-fold range (Figure 5A).
Considering the stark variation in the miRNA profiles between the salt-susceptible and salt-tolerant varieties, an in depth analysis was performed to understand the deviations in the expression patterns (Figure 5B). The changes in expression patterns of the miRNAs between the two varieties were studied using PBSL/PBNL as the reference set. A complex pattern emerged indicating that majority of the salt-induced PB miRNAs were down-regulated or absent in the PK datasets (Figure 5B outer circle). Within the 35 PB miRNAs showing 0–2 fold up-regulation, 9 were down-regulated and 6 were absent in PK. Amongst them osa-miR390-5p was >6-fold down-regulated while osa-miR5817 and osa-miR156j-3p were >4-fold down-regulated in PK. Similarly, 8 of the 26 miRNAs showing 2–4 fold up-regulation in PB were down-regulated in PK while 5 were absent. Likewise 27 of the 70 miRNAs, including osa-miR535-5p and osa-miR168a-5p, which were 0–2 fold down-regulated in PB, were up-regulated in PK, whereas expression of 4 miRNAs (osa-miR2871a-5p, osa-miR5814, osa-miR5825, osa-miR5161) were further down-regulated to more than 6-fold. Similarly, of 50 miRNAs which were down-regulated by 2-4 fold, 10 miRNAs (osa-miR166k-3p, osa-miR166l-3p, osa-miR167d-5p, osa-miR167e-5p, osa-miR167f, osa-miR167g, osa-miR167h-5p, osa-miR167i-5p, osa-miR167j, osa-miR1425-5p) were up-regulated by >6-fold while 3 members of osa-miR444 (osa-miR444a-3p.2, osa-miR444d.2, osa-miR444e) were up-regulated by 4–6 fold.
However, comparisons made against the GLSL/GLNL datasets revealed a much more constant pattern indicating a typical 0–2 fold up-regulation in the expression patterns of most of the salt deregulated miRNAs (Figure 5B, inner circle). It was observed that among the 61 PB miRNAs exhibiting 0–4 fold up-regulation, 18 miRNAs were down-regulated or not expressed in the GL. This set included osa-miR5158 which showed a >6-fold down-regulation. Also some miRNAs like osa-miR2121a,b, osa-miR444a-5p, osa-miR166i-3p, osa-miR396e-3p, which exhibit >4-fold up-regulation in the wild type plants were down-regulated in the GL transgenics. Similarly, out of 177 miRNAs showing down-regulation in PBSL/PBNL, 122 were an up-regulated in GLSL/GLNL. Among them osa-miR1320-3p was up-regulated by more than 4-fold.
3.2 mRNA Target Expression Analysis
The degradome data sets were used to predict the targets of the miRNAs identified in the six libraries (Supplementary File-2). The analysis identified 269 targets for the 171 miRNAs differentially regulated across libraries. The targeted sequences were searched in the RNAseq data to check for their digital expression status. This validated the expression profiles of 144 transcripts targeted by 81 miRNA across different degradome datasets (Table 2). However, 28 predicted target transcripts were not found in the RNAseq datasets. Fifteen percentage targets were dominantly active in transcription regulation by encoding important transcription factors belonging to 14 diverse families such as MYB and NAC. These transcription factors are known to play an active role in response to abiotic stresses (Figure 6A).
|miRNA ID||Regulation up/down||Target gene||Regulation up/down||Target description|
|osa-miR1320-5p||Down||LOC_Os05g47550.1||Up||ANTH/ENTH domain containing protein, putative, expressed|
|osa-miR1425-5p||Down||LOC_Os10g35640.1||Down||Rf1, mitochondrial precursor, putative, expressed|
|osa-miR1425-5p||Down||LOC_Os10g08580.1||Down||FAD binding domain of DNA photolyase domain containing protein, expressed|
|osa-miR160a,b,c,e-5p||Down||LOC_Os06g47150.3||Up||Auxin response factor 18, putative, expressed|
|osa-miR162a||Down||LOC_Os03g02970.1||Up||Dicer, putative, expressed|
|osa-miR162a||Down||LOC_Os02g27030.1||Up||Cysteine proteinase 1 precursor, putative, expressed|
|osa-miR162b||Down||LOC_Os05g51650.1||Up||LSM domain containing protein, expressed|
|osa-miR162b||Down||LOC_Os03g07820.1||Down||Exostosin family protein, putative, expressed|
|osa-miR164d||Up||LOC_Os12g41680.1||Up||No apical meristem protein, putative, expressed|
|osa-miR166l-3p||Down||LOC_Os03g01890.2||Down||START domain containing protein, expressed|
|osa-miR166k-3p||Down||LOC_Os06g47230.1||Up||Coiled-coil domain-containing protein 72, putative, expressed|
|osa-miR167a-5p,b||Down||LOC_Os02g06910.1||Down||Auxin response factor, putative, expressed|
|osa-miR171b||Down||LOC_Os02g44360.1||Up||Scarecrow, putative, expressed|
|osa-miR172d-5p||Down||LOC_Os03g27310.1||Down||Histone H3, putative, expressed|
|osa-miR172a||Down||LOC_Os04g55560.3||Down||AP2 domain containing protein, expressed|
|osa-miR172d-3p||Down||LOC_Os06g04020.1||Down||Histone H1, putative, expressed|
|osa-miR172a||Down||LOC_Os08g39630.1||Down||Helix-loop-helix DNA-binding domain containing protein, expressed|
|osa-miR1856||Down||LOC_Os07g27810.1||Down||Fiber protein Fb34, putative, expressed|
|osa-miR1878||Up||LOC_Os03g15890.3||Up||RNA recognition motif containing protein, expressed|
|osa-miR2863c||Up||LOC_Os05g11780.1||Up||Mitochondrial carrier protein, putative, expressed|
|osa-miR2871b||Down||LOC_Os10g13810.1||Down||Glycosyltransferase family 43 protein, putative, expressed|
|osa-miR2873a||Down||LOC_Os10g10990.3||Down||Transcription initiation factor IIF, alpha subunit domain containing protein, expressed|
|osa-miR2873b||Down||LOC_Os12g19381.1||Down||Ribulose bisphosphate carboxylase small chain, chloroplast precursor, putative, expressed|
|osa-miR393a||Down||LOC_Os05g05800.1||Up||OsFBL21 – F-box domain and LRR containing protein, expressed|
|osa-miR393a||Down||LOC_Os04g32460.2||Down||OsFBL16 – F-box domain and LRR containing protein, expressed|
|osa-miR396a-5p||Up||LOC_Os03g51970.1||Down||Growth regulating factor protein, putative, expressed|
|osa-miR398b||Down||LOC_Os07g46990.1||Up||Copper/zinc superoxide dismutase, putative, expressed|
|osa-miR444b.2c.1,c.2||Down||LOC_Os02g36924.1||Down||OsMADS27 – MADS-box family gene with MIKCc type-box, expressed|
|osa-miR444b.1,c.1||Down||LOC_Os02g49090.1||Up||WD domain, G-beta repeat domain containing protein, expressed|
|osa-miR444b.2,c.1||Down||LOC_Os02g49840.2||Up||OsMADS57 – MADS-box family gene with MIKCc type-box, expressed|
|osa-miR444d.2||Down||LOC_Os02g49840.2||Up||WD domain, G-beta repeat domain containing protein, expressed|
|osa-miR444a-3p.2,b.1||Down||LOC_Os02g49840.3||Down||OsMADS57 – MADS-box family gene with MIKCc type-box, expressed|
|osa-miR444f||Down||LOC_Os02g49840.4||Up||OsMADS57 – MADS-box family gene with MIKCc type-box, expressed|
|osa-miR444b.2||Down||LOC_Os04g38780.1||Down||Transcription factor, putative, expressed|
|osa-miR444d.1||Up||LOC_Os07g48640.2||Up||Short-chain dehydrogenase/reductase, putative, expressed|
|osa-miR444a-3p.2||Down||LOC_Os08g06510.1||Up||Zinc finger, C3HC4 type domain containing protein, expressed|
|osa-miR444a-5p||Up||LOC_Os08g25624.2||Up||Phosphate/phosphate translocator, putative, expressed|
|osa-miR444c.2||Down||LOC_Os08g33488.1||Up||OsMADS23 – MADS-box family gene with MIKCc type-box, expressed|
|osa-miR444a-3p.1||Up||LOC_Os11g47809.1||Up||metallothionein, putative, expressed|
|osa-miR535-5p||Down||LOC_Os06g49010.1||Down||OsSPL12 – SBP-box gene family member, expressed|
|osa-miR5788||Down||LOC_Os12g02330.2||Up||LTPL13 – Protease inhibitor/seed storage/LTP family protein precursor, expressed|
|osa-miR5788||Down||LOC_Os11g02400.1||Up||LTPL8 – Protease inhibitor/seed storage/LTP family protein precursor, expressed|
|osa-miR5813||Down||LOC_Os03g10340.1||Up||40S ribosomal protein S3a, putative, expressed|
|osa-miR812v||Down||LOC_Os02g07960.4||Up||STRUBBELIG-RECEPTOR FAMILY 3 precursor, putative, expressed|
miRNA analysis identified specific molecules whose expression was restricted to a single library and they targeted essential metabolic enzymes and could be responsible for regulating the plants’ physiology. The PBNL specific miRNAs down-regulated stress-responsive dehydrin (LOC_Os11g26780.1 targeted by osa-miR5793), beta-amylase (LOC_Os08g09250.2 targeted by osa-miR1861c), plant-specific domain TIGR01589 containing protein (LOC_Os05g38680.1 targeted by osa-miR5519) and peptidase T1 family (LOC_Os02g08520.1 targeted by osa-miR5519). The PKNL specific osa-miR5532 targets myb/SANT domain protein (LOC_Os01g08680.1), which was found to be associated with QTL AQFW179 for spikelet fertility. This transcript was found to up-regulated in PKSL RNAseq datasets. Similarly, 9 PBSL preferential miRNAs were predicted to target 8 transcripts including a calcium/calmodulin dependent protein kinase (LOC_Os03g22050.4), ribosome inactivating protein (LOC_Os01g06740.1) and plasma membrane ATPase, (LOC_Os12g44150.1) involved in oxidative phosphorylation pathway. However, these transcripts were up-regulated in PKSL. The negative regulation of these transcripts can be clearly associated with the salt-sensitive phenotype of PBNL.
Gene ontology (GO) analysis of the predicted targets revealed their functional significance. Fifty three percentage targets were involved in regulation of biological processes followed by cellular process (Figure 6B–D).
Considering the negative regulation of the targets by the miRNA, the anti-correlation between the expressions of both molecules was analyzed. It was observed that 51 miRNAs were exhibiting a relatively similar deregulation pattern in the PKSL/PKNL and GLSL/GLNL datasets. For instance Osa-miR398b is up-regulated in PBSL while it is down-regulated in PKSL and GLSL, respectively. It regulates LOC_Os01g43540.1, a putative copper/zinc superoxide dismutase, which is involved in oxidative stress response. The target transcript was down-regulated in PBSL and up-regulated in PKSL and GLSL tissues. This explains the salt-tolerant physiologies in relation to earlier studies that a reduced expression of osa-miR398b increases the expression of CSD1 and CSD2 mRNA resulting in enhanced oxidative stress tolerance .
Some variation in the miRNA regulation was also observed as exemplified by osa-miR1870-3p. It is involved in regulating Glyoxylate and dicarboxylate metabolism pathway specifically associated with glutamine synthesis. The miRNA shows salt-induced 2-fold up-regulation in GL but a 0.9-fold down-regulation in PK while there is no significant deregulation in PB. Similarly, osa-miR166k-3p and osa-miR166g-3p target LOC_Os03g43930, a START domain containing protein, which is a putative HD-ZIP (homeobox-leucine zipper protein) transcription factor involved in response to abiotic stress and abscisic acid. Both miRNAs are down-regulated in PKSL and PBSL, whereas their corresponding target is up-regulated. However, the fold up-regulation of the target transcript is more in PKSL as compared to PBSL.
3.3 Identification isomiRs and Expression in Different Experimental Condition
The alignment based analysis of sequencing reads to the precursor miRNA templates identified specific sRNA sequences that align outside the canonical miRNA origins. It has been reported that during biogenesis the pre-miRNA processing region can shift resulting in variation at either ends . The new molecules containing pre-miRNA sequence dependent modifications are termed as “isomiRs”. It was observed that most of the isomiRs contained processing modifications at 3′ with respect to the canonical miRNA sequence (Figure 7A). This largely included deletion of nucleotides and in few cases addition of nucleotides (Figure 7B). In our analysis 103 isomiRs were identified in the different libraries (Supplementary File-3).
Some isomiRs resulted from modifications at both 5′ and 3′ end. A comparative analysis indicated that the number of molecules showing modifications at both ends (mixed modifications) was less in PBNL but their volume considerably increases in PBSL (Figure 7A). The modification profile in GLNL and PKNL was similar to PBSL. On perception of salt-stress the mixed modifications increased in GLSL while not much change was seen in PKSL (Figure 7A). This indicates that PK seems to be pre-adapted to regulate gene expression under salt-stress. The GL express only one pathway to regulate salt-stress and this alone may not be enough to control the increase in the isomiRs.
Several dominant isomiR families were identified that showed differential expression under specific conditions (Figure 7C). Some iso-miRs were present in a specific library (Figure 7C) like iso-miR2873a (addition of A at 5′ and AC at 3′) was abundant in PKSL alone. Likewise iso-miR396a,b-5p-1 (TG deletion at 3′ end) accumulated in PBNL whereas other variants of iso-miR396 family like iso-miR396a,b-5p-2 (CTG deletion at 3′ end) iso-miR396a,b-5p-3 (ACTG deletion at 3′ end) and iso-miR396a,b-5p-4 (G deletion at 3′ end) were expressed in specific libraries. It was interesting to note that there was a deviation in the transcripts targeted by miR396 and its isomiRs, though GRF transcript emerged as a common target. The biological significance of isomiRs is not clear, but it is obvious that addition or deletion at the 5′ end and/or 3′ end of the miRNA can affect their AGO sorting  and alter the targeted transcripts.
3.4 Identification of 3′ and 5′ End Modifications
The post-processing editing mechanisms also modulate the mature miRNA sequence by deleting, adding and substituting the nucleotides in a pre-miRNA template-independent manner , . Such modifications have functional importance and some were shown to regulate the stability of miRNAs. The end modifications also affect the miRNA: target interactions, so they may have potentially important role in stress regulation. About 1.5 % of the salt deregulated miRNAs were found to contain end modifications and 110 modification events were identified in at least one library (Supplementary File-4). Among the end modified reads 83 % were trimmed at 3′ end and 6 % were trimmed at 5′ end. Whereas around 7 % and 4 % reads contained additional bases at 3′ and 5′ end, respectively. The analysis showed that sequences containing A at their 3′ end are modified in all ways (Figure 8). The sequences containing AAA or ATC or GCT or TCC at their 3′ end are predominantly trimmed from their 3′ end while sequences containing TCCA at their 3′ end are predominantly trimmed from their 5′ end (Figure 8). The addition of bases at 3′ end normally includes TTT or TGGA while AAG is particularly added at 5′ end. It is hypothesized that the first nucleotide at the 5′ end is responsible for AGO sorting and the terminal nucleotides play an important role in strand selection during RISC loading , . Thus end modifications influence the putative function of the miRNAs.
Among the 51 salt-deregulated miRNAs 21 miRNAs were modified at both ends. It was observed that osa-miR167c-5p which targets auxin response transcription factor (LOC_Os06g46410.2) was the most modified miRNA, as 4080 modified reads aligned to it. It was observed that the 5′ end modifications were comparatively less than 3′ events. This suggested that the modification patterns were not uniform and various members of a miRNA family members contained different type and volume of modifications.
Adaptation to salt-stress involves the interplay of several genes associated with various biochemical and cellular pathways and the miRNAs are key modulators of gene expression. This analysis identified 443 mature rice miRNAs from six sRNA libraries. It was seen that salt-stress induced the expression of a large number of miRNAs in the salt-tolerant PK while repressed the expression of several miRNAs in salt-sensitive PB. The up-regulation of the miRNAs seemed to play an important role in enhancing the tolerance of the plants to salt-stress. This was supported by the observations made with the salt-tolerant GL transgenics, which exhibited a large pool of up-regulated miRNAs. The targets of the deregulated miRNAs were searched across the degradome datasets and their expression status was checked in the RNAseq data. The identified miRNAs and their predicted targets exhibited varied expression profiles across the rice varieties. Anti-correlation between the expressions of miRNAs and their corresponding targets was observed in several cases.
The analysis of the datasets also identified several isomiRs and post-processing end modifications of the canonical miRNAs. The addition and/or deletion of nucleotides possibly alter their AGO recruitment and the targeted transcripts. Expression values of such molecules also varied across datasets, indicating that these were not random events but were possibility genetically controlled. The significance of these variations in relation to plant physiology needs to be thoroughly understood. Further in depth analysis of the data is required to identify the regulatory nodes operative under salt-stress.
Thus comparative and integrated analysis of different types of deep sequencing data sets adopted in this study produces a global picture of the regulatory networks operative in rice plants under salt stress. The global profiles of miRNAs and their targets clearly indicate a difference in the genetic regulation of salt-susceptible and salt-tolerant rice varieties. By comparing the data within the genetically similar backgrounds using the Gly-transgenics in which salt-tolerance was artificially engineered, the changes in the regulatory networks are apparent. It also indicates that manipulating one pathway alone may not be sufficient to completely alter the physiology of the plants.
Funding source: Department of Biotechnology, Ministry of Science and Technology
Award Identifier / Grant number: BT/PR628/AGR/36/674/2011
Funding statement: The research was supported by financial grants received from the Department of Biotechnology, Ministry of Science and Technology, Government of India (No. BT/PR628/AGR/36/674/2011).
The authors thank Prof. S.K. Sopory and Dr. S.L. Singla-Pareek for providing the Gly transgenics. The authors thank Dr. Deepti Mittal, Dr. Osmani Chacon and Dr. Neha Sharma for library preparation.
Conflict of interest statement: Authors state no conflict of interest. All authors have read the journal’s Publication ethics and publication malpractice statement available at the journal’s website and hereby confirm that they comply with all its parts applicable to the present scientific work.
 Barrera-Figueroa BE, Wu Z, Liu R. Abiotic stress-associated microRNAs in plants: discovery, expression analysis, and evolution. Frontiers Biol. 2013;8:189–97.10.1007/s11515-012-1210-6Search in Google Scholar
 Seki M, Narusaka M, Ishida J, Nanjo T, Fujita M, Oono Y, et al. Monitoring the expression profiles of 7000 Arabidopsis genes under drought, cold and high salinity stresses using a full length cDNA microarray. Plant J. 2002;31:279–92.10.1046/j.1365-313X.2002.01359.xSearch in Google Scholar PubMed
 Mae T. Physiological nitrogen efficiency in rice: nitrogen utilization, photosynthesis, and yield potential. In: Ando T, editors. Plant nutrition for sustainable food production and environment. Dordrecht: Springer Netherlands, 1997:51–60.10.1007/978-94-009-0047-9_5Search in Google Scholar
 Serrano R, Mulet JM, Rios G, Marquez JA, De Larrinoa IF, Leube MP, et al. A glimpse of the mechanisms of ion homeostasis during salt stress. J Exp Bot. 1999;50:1023–36.10.1093/jxb/50.Special_Issue.1023Search in Google Scholar
 Hasanuzzaman M, Nahar K, Alam MM, Roychowdhury R, Fujita M. Physiological, biochemical, and molecular mechanisms of heat stress tolerance in plants. Int J Mol Sci. 2013;14:9643–84.10.3390/ijms14059643Search in Google Scholar PubMed PubMed Central
 Kumar R, Mustafiz A, Sahoo KK, Sharma V, Samanta S, Sopory SK, et al. Functional screening of cDNA library from a salt tolerant rice genotype Pokkali identifies mannose-1-phosphate guanyl transferase gene (OsMPG1) as a key member of salinity stress response. Plant Mol Biol. 2012;79:555–68.10.1007/s11103-012-9928-8Search in Google Scholar PubMed
 Singh AK, Ansari MW, Pareek A, Singla-Pareek SL. Raising salinity tolerant rice: recent progress and future perspectives. Physiol Mol Biol Plants. 2008;14:137–54.10.1007/s12298-008-0013-3Search in Google Scholar PubMed PubMed Central
 Uddin MI, Qi Y, Yamada S, Shibuya I, Deng XP, Kwak SS, et al. Overexpression of a new rice vacuolar antiporter regulating protein OsARP improves salt tolerance in tobacco. Plant Cell Physiol. 2008;49:880–90.10.1093/pcp/pcn062Search in Google Scholar PubMed
 Verma D, Singla-Pareek SL, Rajagopal D, Reddy M, Sopory S. Functional validation of a novel isoform of Na+/H+ antiporter from Pennisetum glaucum for enhancing salinity tolerance in rice. J Biosci. 2007;32:621–8.10.1007/s12038-007-0061-9Search in Google Scholar PubMed
 Chinnusamy V, Zhu J, Zhu JK. Salt Stress Signaling and Mechanisms of Plant Salt Tolerance.. In: Setlow Jane K., editors. In Genetic Engineering: Principles and Methods. Boston, MA: Springer US, 2006:141–177.978-0-387-25855-3.10.1007/0-387-25856-6_9Search in Google Scholar PubMed
 Kumari S, nee Sabharwal VP, Kushwaha HR, Sopory SK, Singla-Pareek SL, Pareek A. Transcriptome map for seedling stage specific salinity stress response indicates a specific set of genes as candidate for saline tolerance in Oryza sativa L. Funct Integr Genomics. 2009;9:109.10.1007/s10142-008-0088-5Search in Google Scholar PubMed
 Singla-Pareek S, Reddy M, Sopory S. Genetic engineering of the glyoxalase pathway in tobacco leads to enhanced salinity tolerance. Proc Natl Acad Sci U S A. 2003;100:14672–7.10.1073/pnas.2034667100Search in Google Scholar PubMed PubMed Central
 Thornalley PJ. The glyoxalase system: new developments towards functional characterization of a metabolic pathway fundamental to biological life. Biochem J. 1990;269:1.10.1042/bj2690001Search in Google Scholar PubMed PubMed Central
 Saxena M, Roy SD, Singla-Pareek SL, Sopory SK, Bhalla-Sarin N. Overexpression of the glyoxalase II gene leads to enhanced salinity tolerance in Brassica juncea. Open Plant Sci J. 2011;5:23–8.10.2174/1874294701105010023Search in Google Scholar
 Aukerman MJ, Sakai H. Regulation of flowering time and floral organ identity by a microRNA and its APETALA2-like target genes. Plant Cell. 2003;15:2730–41.10.1105/tpc.016238Search in Google Scholar PubMed PubMed Central
 Mallory AC, Reinhart BJ, Jones‐Rhoades MW, Tang G, Zamore PD, Barton MK, et al. MicroRNA control of PHABULOSA in leaf development: importance of pairing to the microRNA 5′ region. EMBO J. 2004;23:3356–64.10.1038/sj.emboj.7600340Search in Google Scholar PubMed PubMed Central
 Lutts S, Kinet J, Bouharmont J. NaCl-induced senescence in leaves of rice (Oryza sativa L.) cultivars differing in salinity resistance. Ann Bot. 1996;78:389–98.10.1006/anbo.1996.0134Search in Google Scholar
 Ding D, Zhang L, Wang H, Liu Z, Zhang Z, Zheng Y. Differential expression of miRNAs in response to salt stress in maize roots. Ann Bot. 2009;103:29–38.10.1093/aob/mcn205Search in Google Scholar PubMed PubMed Central
 Alvarez-Gerding X, Espinoza C, Inostroza-Blancheteau C, Arce-Johnson P. Molecular and physiological changes in response to salt stress in Citrus macrophylla W plants overexpressing Arabidopsis CBF3/DREB1A. Plant Physiol Bioch. 2015;92:71–80.10.1016/j.plaphy.2015.04.005Search in Google Scholar PubMed
 Mostofa MG, Fujita M, Tran LSP. Nitric oxide mediates hydrogen peroxide-and salicylic acid-induced salt tolerance in rice (Oryza sativa L.) seedlings. Plant Growth Regul. 2015;77:265–77.10.1007/s10725-015-0061-ySearch in Google Scholar
 Mostofa MG, Hossain MA, Fujita M. Trehalose pretreatment induces salt tolerance in rice (Oryza sativa L.) seedlings: oxidative damage and co-induction of antioxidant defense and glyoxalase systems. Protoplasma. 2015;252:461–75.10.1007/s00709-014-0691-3Search in Google Scholar PubMed
 Alvarez-Gerding X, Cortés-Bullemore R, Medina C, Romero-Romero JL, Inostroza-Blancheteau C, Aquea F, et al. Improved salinity tolerance in Carrizo Citrange Rootstock through overexpression of glyoxalase system genes. BioMed Res Int. 2015;2015:827951.10.1155/2015/827951Search in Google Scholar PubMed PubMed Central
 Gaidatzis D, van Nimwegen E, Hausser J, Zavolan M. Inference of miRNA targets using evolutionary conservation and pathway analysis. BMC Bioinform. 2007;8:69.10.1186/1471-2105-8-69Search in Google Scholar PubMed PubMed Central
 Bokvaj P, Hafidh S, Honys D. Transcriptome profiling of male gametophyte development in Nicotiana tabacum. Genome Data. 2015;3:106–11.10.1016/j.gdata.2014.12.002Search in Google Scholar PubMed PubMed Central
 Tan GC, Chan E, Molnar A, Sarkar R, Alexieva D, Isa IM, et al. 5′ isomiR variation is of functional and evolutionary importance. Nucleic Acids Res. 2014;42:9424–35.10.1093/nar/gku656Search in Google Scholar PubMed PubMed Central
 Mittal D, Mukherjee SK, Vasudevan M, Mishra NS. Identification of tissue-preferential expression patterns of rice miRNAs. J Cell Biochem. 2013;114:2071–81.10.1002/jcb.24552Search in Google Scholar PubMed
 Ma Z, Coruh C, Axtell MJ. Arabidopsis lyrata small RNAs: transient MIRNA and small interfering RNA loci within the Arabidopsis genus. Plant Cell. 2010;22:1090–103.10.1105/tpc.110.073882Search in Google Scholar PubMed PubMed Central
 Patel RK, Jain M. NGS QC toolkit: a platform for quality control of next-generation sequencing data. In: Nelson Karen E, editors. Encyclopedia of Metagenomics. Boston, MA: Springer US, 2015:544–8.Search in Google Scholar
 Griffiths-Jones S, Grocock RJ, Van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(Suppl 1):D140–D4.10.1093/nar/gkj112Search in Google Scholar PubMed PubMed Central
 Kozomara A, Griffiths-Jones S. miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011;39:D152–7.10.1093/nar/gkq1027Search in Google Scholar PubMed PubMed Central
 Peng T, Lv Q, Zhang J, Li J, Du Y, Zhao Q. Differential expression of the microRNAs in superior and inferior spikelets in rice (Oryza sativa). J Exp Bot. 2011;62:4943–54.10.1093/jxb/err205Search in Google Scholar PubMed
 Chen Z, Li F, Yang S, Dong Y, Yuan Q, Wang F, et al. Identification and functional analysis of flowering related microRNAs in common wild rice (Oryza rufipogon Griff.). PLoS One. 2013;8:e82844.10.1371/journal.pone.0082844Search in Google Scholar PubMed PubMed Central
 Kaushik A, Saraf S, Mukherjee SK, Gupta D. miRMOD: a tool for identification and analysis of 5′ and 3′ miRNA modifications in next generation sequencing small RNA data. Peer J. 2015;3:e1332.10.7717/peerj.1332Search in Google Scholar PubMed PubMed Central
 Addo-Quaye C, Miller W, Axtell MJ. CleaveLand: a pipeline for using degradome data to find cleaved small RNA targets. Bioinformatics. 2009;25:130–1.10.1093/bioinformatics/btn604Search in Google Scholar PubMed PubMed Central
 Cao P, Jung K-H, Choi D, Hwang D, Zhu J, Ronald PC. The rice oligonucleotide array database: an atlas of rice gene expression. Rice. 2012;5:1.10.1186/1939-8433-5-17Search in Google Scholar PubMed PubMed Central
 Llorens F, Bañez-Coronel M, Pantano L, del Río JA, Ferrer I, Estivill X, et al. A highly expressed miR-101 isomiR is a functional silencing small RNA. BMC Genomics. 2013;14:104.10.1186/1471-2164-14-104Search in Google Scholar PubMed PubMed Central
 Cloonan N, Wani S, Xu Q, Gu J, Lea K, Heater S, et al. MicroRNAs and their isomiRs function cooperatively to target common biological pathways. Genome Biol. 2011;12:1.10.1186/gb-2011-12-12-r126Search in Google Scholar PubMed PubMed Central
 Ebhardt HA, Tsang HH, Dai DC, Liu Y, Bostan B, Fahlman RP. Meta-analysis of small RNA-sequencing errors reveals ubiquitous post-transcriptional RNA modifications. Nucleic Acids Res. 2009;37:2461–70.10.1093/nar/gkp093Search in Google Scholar PubMed PubMed Central
 Wyman SK, Knouf EC, Parkin RK, Fritz BR, Lin DW, Dennis LM, et al. Post-transcriptional generation of miRNA variants by multiple nucleotidyl transferases contributes to miRNA transcriptome complexity. Genome Res. 2011;21:1450–61.10.1101/gr.118059.110Search in Google Scholar PubMed PubMed Central
 Saraf S, Sanan-Mishra N, Gursanscky NR, Carroll BJ, Gupta D, Mukherjee SK. 3′ and 5′ microRNA-end post-biogenesis modifications in plant transcriptomes: Evidences from small RNA next generation sequencing data analysis. Biochem Biophys Res Commun. 2015;467:892–9.10.1016/j.bbrc.2015.10.048Search in Google Scholar PubMed
The online version of this article offers supplementary material (DOI: https://doi.org/10.1515/jib-2017-0002).
©2017, Neeti Sanan-Mishra, published by De Gruyter, Berlin/Boston
This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.