Mapping the microRNA Expression Profiles in Glyoxalase Over-expressing Salinity Tolerant Rice.

In the recent years, glyoxalase pathway has been an active area of research in both human and plants. This pathway is reported to confer stress tolerance in plants, by modulating the glutathione homeostasis to achieve detoxification of a potent cytotoxic and mutagenic compound, methylglyoxal. The microRNAs (miRNAs) are also reported to play significant role in stress tolerance for plants. However, the cross-talk of miRNAs with the metabolism regulated by glyoxalase in the salinity-tolerance is unexplored. We therefore investigated whether expression profiles of miRNAs are altered in response to glyoxalase overexpression, and if any of these are also responsible for modulating the stress responses of plants. In this study, the Next Generation Sequencing (NGS) was employed to profile miRNA expression levels from glyoxalase overexpressing transgenic lines. The associated targets of differentially expressed miRNAs were predicted and their functional annotation was carried out using Gene Ontology (GO) and KEGG Orthology (KO), which showed their involvement in several crucial biological pathways. The analysis of NGS datasets also identified other isoforms or isomiRs of selected miRNAs, which may have an active role in developing tolerance against salt stress. Different aspects of miRNA modifications were also studied in glyoxalase overexpressing lines.

INTRODUCTION

Improved rice varieties and hybrids developed by researchers all over the world have helped to ameliorate the quality as well as quantity of rice production. However, poor irrigation practices have increased the threat of soil salinity which severely affects the plant developmental processes like the germination of seed, seedling vigor, flowering time and seed set [1]. High soil salt content disrupts the osmotic and ionic balance at both cellular as well as whole-plant levels. All major crops are sensitive to high concentrations of sodium (Na+). The search for genes imparting salt tolerance has identified molecules that are either involved in the stress signaling pathways or are crucial for the structural and functional protection [2, 3] or those that help to restore the osmoticum or ion homeostasis [4-6].

In a classical study, it was shown that Salt Overly Sensitive (SOS) signaling pathway that comprises of three main proteins SOS1, SOS2 and SOS3, acts on ion homeostasis and enables plants to achieve some level of tolerance against high salt [7, 8]. The SOS1 encodes a plasma membrane Na+/H+ antiporter and has a role in regulation of Na+ efflux at cellular level. It also mediates the long distance transport of Na+ [9, 10]. The SOS2 gene encodes a serine/threonine kinase that is turned on by salt stress drawn Ca2+ signals. The SOS3 protein contains myristoylation site at its N-terminus [11]. The binding interaction of SOS3 with SOS2 protein activates the kinase [12], which in turn phosphorylates SOS1 and increases Na+ efflux, thereby decreasing Na+ toxicity [13, 14]. The salt tolerance can also be attained by over-expression of SOS1 in plants [15]. The role of several MAP kinases has also been reported in response to salt stress. Recent reports have shown MPK6 mediated regulation of SOS1 phosphorylation by NaCl and Phosphatidic Acid [16].

Another strategy to keep a check on ion concentration in cytoplasm, involves maintaining the cellular levels of Na+/K+. Both the ions share common mechanism of transport therefore K+ competes with Na+ for entry into the cytoplasm. Under salinity stress an increased amount of Na+ ion results in decreased uptake of K+ [1, 17]. This induces the expression of transporters such as HKT (histidine kinase transporter) and NHX [18-20], which favour a continuous uptake of K+ [18-20]. These transporters belong to large gene families and the role of specific isoforms (e.g. LeNHX3 and LeNHX4) in ion homeostasis is now established in different plants [21].

Among the various strategies the glyoxalase (Gly) system is involved in chemical detoxification by the coordinated action of two enzymes- Glyoxalase I (GlyI, lactoylglutathione lyase; EC4.4.1.5) and Glyoxalase II (GlyII, hydroxyacylglutathione hydrolase; EC 3.1.2.6), which act together to convert methylglyoxal (MG) and other 2-oxoaldehydes in their 2-hydroxyacids with the help of a cofactor Glutathione (GSH) [22]. Under natural physiological conditions, Gly I catalyzes the initial step of the pathway [22] by converting MG to SD-lactoylglutathione (SLG) in presence of GSH. In the next step, catalyzed by Gly II, the cytotoxic SLG is hydrolysed to D-lactate and GSH is released [22, 23]. This pathway has been reported to confer salt stress tolerance in plants by removing the cellular toxicity and regulating the GSH homeostasis [24]. It was shown that with the exposure of plants to various abiotic stresses, concentration of MG increases [25]. In general, most important source of MG is glycolysis at triose-phosphate step and catabolism of amino acids and acetone [26]. High levels of MG is noxious for the plant as it reacts with the major macro-molecules such as DNA, RNA and proteins [27]. The levels of Gly I and Gly II transcripts were shown to increase with the duration of salt treatment in cotyledonary leaves as well as in roots of Brassica juncea and rice [24, 28]. Studies involving over-expression of the Gly genes in plants have confirmed its role in enhancing tolerance towards salt, drought and heavy metal stress [29-31] in many plant species. The mechanism behind the Gly mediated salt tolerance is being studied.

Increasing evidence has established that miRNAs play crucial roles under various abiotic stresses by regulating gene expression. miRNAs are small non-coding RNAs, typically derived from a segment of the genome (pri-miRNA) which lies in the intergenic, or intronic, regions [32, 33]. In plant system the pri-miRNA is processed in the nuclei, through a pathway mediated by DICER-like (DCL) proteins, into hairpin like structure called precursor miRNA (pre-miRNA) which is subsequently cleaved into 20-24 nt long duplexes. This mature miRNA binds to ARGONAUTE (AGO) protein complex termed as RNA Induced Silencing Complex [RISC] [34-38] and targets the complementary transcripts to negatively regulate gene expression by translation repression or degradation of targeted mRNAs [39, 40].

Recently, computational approaches have gained importance for rapid, rigorous and cost effective identification of miRNAs from NGS dataset [41]. The algorithms have been successful in predicting the targets of the miRNAs [42-44]. Identification of miRNAs and their target genes therefore is an important step towards understanding their biological functions. We have standardized computational algorithms for identifying the expression patterns of rice miRNAs in different tissues and these computational results were experimentally validated [45]. In this study, the same algorithms were used to identify the miRNA expression profiles and predict novel miRNAs and their target genes from deep sequenced libraries of Gly over expressing transgenic rice lines. The associated targets of the differentially expressed miRNAs were predicted and their functional annotation was carried out using Gene Ontology (GO) and KEGG Orthology (KO) analysis. The analysis of datasets also identified the isomiRs and modified forms (isoforms) of selected miRNAs. The role of miRNAs and their variants in Gly mediated tolerance against salt stress is discussed.

MATERIALS AND METHODS

Plant Materials and Growth Conditions

Seeds of salinity susceptible, wild type Pusa Basmati lines and salinity tolerant, Pusa Basmati lines overexpressing Gly I and Gly II [46] were surface sterilized with a solution of 10% commercial bleach (0.525% sodium hypochlorite) 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 tissue samples were harvested from normal or 200 mM salt (NaCl) stressed 15-days-old seedlings.

Small RNA Library Preparation

Leaf tissues of 15-days-old wild type (PB) and Gly seedlings (Gly) grown under normal (NL) and salt stressed (SL) conditions were used for RNA isolation. One gram of leaf tissue was rinsed briefly with sterile H2O containing 0.1% (v/v) Diethylpyrocarbonate (DEPC, Sigma) and used for total RNA isolation as described earlier [45]. The small RNA was enriched by LiCl precipitation and used for library construction as described earlier [45]. The libraries thus made were used for deep sequencing on GAII sequencer (Illumina).

Analysis of Sequencing Data

Pre-processing

The data was delivered by sequencer as sequences of 36 bases in length with the base quality scores and read counts. Quality Filtering Analysis (performed by NGSQC tool kit) showed that above 97% reads were of high quality. The raw sequence reads were processed to remove the 3' adapter by searching for 5 bases of the adapter sequence using the criteria of exact match. The reads containing these sequences were further checked and the sequences aligning to the complete adapter sequence were selectively removed from the read. After adapter trimming only those reads were retained which were having read length of 18-24 nt. These sequences were clustered based on sequence identity such that all identical reads in each small RNA library were grouped into single sequence tag each containing information on the unique sequence and its abundance or expression value in the sequenced library.

Identification of Known miRNA

The small RNA sequences of normal and salt stressed Gly overexpressing libraries were matched against the rice sequences in miRBase 21.0, using Bowtie, to identify known miRNAs [47]. The sequences showing 100% match to known mature rice miRNAs were selected. Their digital expression status was generated by normalizing the expression of miRNAs in each sample and expressing the values as transcript per million (TPM) [41, 48].

To observe the differential expression pattern of miRNAs a comparison of normalized expression (NE) across the normal and stressed libraries were performed and fold-change was calculated. The values were plotted as log2-ratio [49, 50].

The miRNAs, which were ≥2 fold deregulated, were considered for subsequent analysis.

Identification of isomiRs and End Modifications of Known miRNAs

miRBase database [51] was used to retrieve the precursor sequences of all the known miRNAs identified in the sequenced libraries. PERL script was used to align the miRNA reads of length range from 17-27 nt to the precursor sequences using the criteria of perfect match. The expression of aligned reads was checked. The reads, which showed higher abundance than that of the corresponding known mature miRNA sequence, were assigned as candidate isomiRs.

For identification of 3’ and 5’ end modifications the sRNA reads were aligned to known miRNA sequences individually using bowtie. These aligned reads were identified and subjected to the miRMOD tool for identification of precursor independent end modifications.

Identification of Novel microRNA

Novel miRNAs with several modifications were identified by following the steps [41]. All the unique tags were aligned to the genome (retrieved from BGI - Indica genome) with zero mismatch, using Bowtie 0.12.7. The reads showing alignment with genome were further aligned on known miRNA (mirBase Rel. 21) [51] and other small non-coding RNA species (retrieved from BGI - Indica genome data) like rRNA, snRNA, sRNA and tRNAs in the Rfam10.1 [52]. The aligned reads were excluded from further analysis. For the remaining reads genomic alignment were fetched and flanking region (200 base up and down stream) were retrieved from genomic sequence. The flanking region were considered as candidate pre-miRNAs and further screened for meeting the structure and energy requirements of plant miRNA [53]. The most abundant read aligning to the stem region of the putative precursor is considered as putative miRNA. To further validate the putative miRNA its corresponding miRNA* sequence was predicted and traced back in sequenced libraries. The reads for which corresponding star sequence was present in any of the libraries were predicted as miRNAs.

QTL Mapping for Novel miRNA Families

The QTL coordinates of rice were retrieved from Gramene QTL database [54]. These coordinates were used to extract the QTL sequences from MSU 6 Rice Genome. Next the precursor sequences of the novel miRNAs and their isoforms were searched in the retrieved QTL sequences for exact match. The smallest and largest associated QTLs of each aligned miRNA were identified and the trait corresponding to the smallest QTL was assigned to the miRNA.

Target Prediction

Potential gene targets for known and novel rice miRNAs were searched using online web tool psRNA Target, Oryza sativa (rice), transcript, TIGR genome cDNA OSA1 Release 5.0 (OSA1R5), using default parameters. Targets mRNA sequences with a score of ≤ 3 were considered as potential targets. The predicted target sequences were searched against PlantTFDB [55] to identify sequences with transcription regulatory activities.

To completely examine the function of the target sequences, GO term and KEGG pathway annotation was performed. The candidate targets were analyzed by Rice Oligonucleotide Array Database [56] to calculate fold-enrichment in each GO category. If the calibrated hypergeometric p-value obtained was below or equal to 0.05, the biological function of the targets could be confirmed. The regulator pathway annotation was performed based on scoring and visualization of the pathways collected in the KEGG database by the help of KEGG Orthology (KO) enrichment analysis from Rice Oligonucleotide Array Database (ROAD) (http://ricearray.org/). Further the number of target transcripts for each miRNA and number of miRNAs targeting a single transcript were calculated and a network view was created to visualize the miRNA and target gene interactions and enrichments using Cytoscape 2.8.3 [57].

RESULTS AND DISCUSSION

Analysis of Small RNA Sequences and miRNA Expression Profiles

Deep sequencing of the small RNA libraries of leaf tissue from Gly seedlings growing under normal (GNL) or salt stressed (GSL) conditions generated total 26371105 and 26780310, reads, consisting of 2898648 and 3385190 unique sequences with their length between 18-nt to 24-nt, respectively (Table ; Supplementary Fig. ). The size distribution analysis of the unique reads indicated that 24-nt sequences represented the largest and the most diverse class with most sequences containing Adenosine or Guanosine residues at the first nucleotide position (Fig. ). The diversity seen in the 24-nt small RNAs is consistent with their origin mainly as short interfering RNAs (siRNAs), which may reflect their role in chromatin maintenance [45, 58]. The overall pattern was similar for both GNL and GSL libraries and also consistent with previous reports from other species such as Arabidopsis thaliana, Medicago truncatula, Solanum lycopersicum, Cucumis sativa, Zea mays and Populus trichocarpa [59-65].

On comparing the total abundance based on the number of reads, the 21-nt sequences emerged as the largest size class (Fig. ) followed by the 24-nt and 20-nt long sequences, respectively. Among the 21-nt and 20-nt long sequences there was abundance of Thymidine at the first nucleotide while among the 24-nt size class the majority of the sequences contained Adenosine residues at the first nucleotide position. The 21-nt size is characteristic for canonical miRNAs and high ratio of abundances/unique reads, reflects the high expression and homogeneity among the sequences indicating an important regulatory role for these small RNAs in plant biology.

The deep sequencing data sets obtained from the Gly over-expressing lines comprised 357 out of 713 mature rice miRNAs listed under miRBase (Release 21). We could identify

319 miRNAs from GNL and 334 miRNAs in the GSL libraries. Among these the 21-nt size class, that is known to regulate post-transcriptional gene expression, emerged as the abundant molecules as observed in (Fig. ). This size distribution pattern is similar to earlier observations on the rice miRNAs [45, 66, 67]. The NGS data sets also provided the digital expression status of each miRNA based on the normalized expression values (Supplementary Table ).

Salt Stress-regulated miRNAs

The detailed coverage of mature miRNAs and their expression levels provided by each sequencing data set enabled a comparison between the normal (NL) and salt stressed (SL) tissues of the wild-type plants and Gly overexpressing (G) lines. A four-way comparison of all libraries depicted that 231 (64.7%) of the 357 known miRNAs identified in Gly lines were common to all libraries while expression of 72 (~20%) known miRNAs was induced in the GNL library. Among the common miRNAs 59 miRNAs were up regulated and 13 miRNAs were down regulated on salt stressed. It was hypothesized that the targets for these miRNAs may play a significant role in imparting salt tolerance to the Gly lines.

To understand the significance of salt stress induced deregulations of the miRNA molecules in the wild-type and Gly overexpressing lines, the fold deregulations were plotted and compared (Fig. ). It was observed that greater number of miRNAs were down regulated by salt stress in the WT plants (SLPB/NLPB) as compared to the Gly lines (GSL/GNL), while the number of up regulations was more in the Gly transgenics (Fig. ). This indicates that over-expression of Gly braces the miRNA machinery to regulate the damage induced by salt stress. This was evident on comparing the deregulations between the normal tissues of WT and Gly lines. A uniform distribution of down regulated and up regulated miRNAs was observed (GNL/NLPB), however, imposition of stress caused an upward shift in the profile (GSL/SLPB) (Fig. ).

It was observed that 11 miRNAs were specifically up regulated and 18 miRNAs were down regulated in the GNL tissues as compared to wild-type (NLPB) tissues (Fig. ). Among these osa-miR5828 and osa-miR440 were 4-fold down regulated, leading to the up regulation of a F-box protein and potassium channel AKT1, respectively. This indicates the activation of pathways favouring controlled degradation of cellular proteins [68] and an import of K+ ions [69] both of which can help the plants to override the toxic effects of Na+ accumulation. The levels of osa-miR1846 were 4-fold up regulated and it is predicted to target an anion-transporting ATPase [70]. In the GSL tissues as compared to salt stressed wild-type (SLPB) tissues 4 miRNAs were specifically up regulated and 17 miRNAs were down regulated (Fig. ).

Several previous studies have highlighted many salt stress associated miRNAs, 39 of which were identified, in the Gly tissues in this study [41, 71, 72]. Expression levels of some miRNAs like osa-miR398 [73], osa-miR812 [74], osa-miR528 [75], osa-miR397 [76], osa-miR167 [77], osa-miR160 [77], osa-miR408 [78] etc. were highly up regulated in the salt stressed Gly plants (GSL). This supports the observation that the Gly over-expressing plants are better prepared to tolerate stress.

Prediction of Novel miRNAs and their Targets from Gly Plants

Our analysis predicted 87 novel miRNAs from the Gly libraries. A search for their putative miRNA* sequences in the libraries [79] validated 31 sequences as putative miRNAs. The sequence and read values of all the novel miRNAs is provided in (Table ) while their structures are represented in (Supplementary Fig. ). The complete information including miRNA, miRNA* and precursor sequences with their chromosome location for all the isoforms of these novel sequences is provided in (Supplementary Table ). Many of the novel miRNA sequences had more than one hit on genome. The read length distribution analysis revealed that the novel miRNAs were also enriched in the 21-nt size class [58].

Among the predicted novel miRNAs, 8 sequences were found to exhibit homology to the pre-miRNA sequences of known miRNAs reported in miRBase. These included osa_13, osa_24 and osa_25, which were identified as novel members of osa-miR2120, osa-miR2863 and osa-miR5788 families, respectively. The osa_26 - osa_29 were found to arise from the 3p arm of osa-miR2866, osa-miR169i,j,m and osa-miR1846e whereas, osa_30 and osa_31 mapped to the 5p arm of osa-miR159f and osa-miR171h.

The annotation of a miRNA is incomplete without the identification of the targets which they regulate. The plant miRNAs are nearly perfect complementary to their targets transcripts, so the computational approach of identifying the targets has been very successful and effective. The targets for the novel miRNAs incuded a large number of unannotated “hypothetical proteins” and “expressed proteins”. On considering only the annotated transcripts 712 miRNA: mRNA interactions could be identified (Supplementary Table ) and these included 9 transcription factor families like GRAS, G2-like, WRKY, MIKC and FAR1 associated with salt stress.

The novel miRNAs also exhibited stress-induced deregulations with 6 sequences down regulated and 18 sequences up regulated under salt stress. The novel miRNA osa-miR23 which targets pentatricopeptide repeat protein, PPR1106-17, was found to be 4-fold down regulated in GSL. This indicates an up regulation of the corresponding PPR transcripts which is supported by the fact that the PPR proteins play an important role in mitigating the stress responses in plants by regulating the ABA signaling pathway [80].

Enrichment Analysis

Gene Ontologies (GO) provide a controlled vocabulary to assess the biological process (P), molecular function (F) and

cellular component (C) a gene product is putatively contributing [81, 82]. KEGG Orthology (KO) consists of manually assigned ortholog groups which belongs to KEGG pathway and BRITE hierarchy nodes, which acts as a basis for the representation of pathway maps in KEGG and BRITE functional hierarchies [83]. Hence, to identify biological patterns in the list of miRNA target, analyses of both GO and KO was performed using Rice Oligonucleotide Array Database (ROAD).

It was observed that within the miRNA target list 99 GO terms (having p-value >= 0.05) were enriched including 51 categories under Molecular Function, 39 terms under biological process and 9 in cellular components (Fig. ). KO enrichment analysis showed the involvement of the target

transcripts under various biological pathways some of which were associated with the cellular signaling pathways such as PPAR signaling pathway, MAPK signaling pathway etc.

One such example is provided by osa-miR1879 which targets an antioxidant gene Catalase (LOC_Os06g51150) and is >5 fold down regulated in the Gly libraries. This target transcript is associated with “MAPK signaling pathway” and “Glyoxylate and dicarboxylate metabolism”. Thus the miRNA down regulation will lead to accumulation of catalase and this may work as an adaptive mechanism to quench the H2O2 under salt stress and offer defense against oxidative damage [84]. The catalase enzyme is also a down-stream component of the Mitogen-activated protein kinase (MAPK) cascade. This is an evolutionarily conserved signal transduction component concerned with conveying the extracellular signals to the nucleus for suitable cellular alterations [85]. The cascade constitutes three important components, a MAPK kinase (MAPKKK), a MAPK kinase (MAPKK) and a MAPK connected to each other by sequential phosphorylation reactions. It is active in response to various stresses as exemplified by the association of MAPKK4, MAPKK6, MAPKK1 MEKK1and OsMAPK5 in salt stress [78, 86, 87]. The signaling pathway is known to regulate phosphorylation of SOS1 and induce activity of catalase [88, 89].

QTL Identification

To understand the biological role of the novel miRNAs QTL (Quantitative Trait Locus) analysis was performed. This search resulted in total 1506 QTL associations for 14 novel precursors (osa_1, osa_2_iso1,4, osa_8_iso2,3, osa_10, osa_16, osa_17_iso2, osa_19_iso1,3,7, osa_20_ iso2, osa_24 which includes three isoform of osa_19, two of osa_2 and two of osa_8). The smallest QTL was considered as the basis of inferring the phenotypic traits associated with the novel miRNAs.

It was found that osa_24 was linked to the abiotic stress QTL trait named “phosphorus uptake” and osa_19_iso1 linked to the QTL trait named “Salt Sensitivity”. osa_19_iso1 showed more than 7 fold down regulation in GSL datasets, indicating a role in the plant’s stress tolerance mechanism. Among the other traits osa_1 associated with “leaf height”, osa_2_iso1 and osa_16 were associated with “Carbohydrate content” of plant, osa_2_iso4 was linked to “vegetative and reproductive growth time” and osa_20_iso2 was associated with “lodging incidence”. It was observed that while osa_8_iso2 was associated with “cold tolerance” its isoform associated with “days to heading” (Fig. ). The details of all the identified QTLs for these novel miRNAs is provide in (Supplementary Table ). These findings suggest that any change in expression levels of these miRNAs may alter the target transcripts regulating the related responses in plants.

Identification of isomiRs and their Expression Profiles

The introduction of deep sequencing techniques and sequence-based analysis of small RNA libraries has now made it possible to identify and analyze various end modifications on the miRNAs [79, 90]. It is observed that miRNAs do not exist as a single typical canonical sequence but they have a range of sequence variants, which are termed as “isomiRs”. The isomiRs may be produced as a result of processing variations, at unconventional sites, during the biogenesis of miRNAs as the sequences can be well aligned to the pre-miRNA templates [91]. This shows that the miRNA as published in miRBase is not always the most common sequence; other cutting variants are also encountered under different physiological conditions. The biological significance of the end modifications is not clear though it is evident that additions or deletions at the 5’ end of the miRNA can affect its AGO sorting [85, 86] and can lead to alteration in the target transcripts as well [92].

We observed that the isomiRs could be identified for 81 known miRNAs using the pre-miRNA as the template sequence. Normally more than one sequence variant was present for each miRNA so we selected the most abundant form as the isomiR (Supplementary Table ). The alignments of selected miRNAs and isomiRs (osa-miR396, osa-miR156, osa-miR820 and osa-miR167) are represented in (Fig. ). Subsequent analysis showed that 54 sequences had higher expression in the GNL library as compared to their corresponding canonical mature sequences (Fig. ). Among these expression levels of 14 were more than 4 fold high (Fig. ). Similarly 60 isomiRs had higher expression in GSL library with 15 sequences having greater than 4 fold expression.

It was found that isomiRs of osa-miR160f-3p, osa-miR815a-c, osa-miR5151 and osa-miR397b were present only in Gly libraries (GNL and GSL) and their expression was further induced by salt stress. The expression of isomiR of osa-miR397b showed more than 12 fold higher expression in the GSL library. However osa-miR397b has been reported as a stress-down regulated miRNA [93, 94] and its expression decreased further in the salt stressed Gly plants. miR397 is predicted to target β-fructofuranosidase, an enzyme involved in starch and sucrose metabolism. This suggests that sucrose metabolism, would be up regulated under stress and this may help in maintaining the cellular osmotic balance by accumulation of sugar osmolytes [95]. The miR397 is also predicted to target the laccase gene family, which was reported to reduce root growth under dehydration [96]. The isomiR of osa-miR397 carry deletions at the 5’ end resulting in the change of target to ascorbate oxidase, which regulates the amount of cellular ascorbic acid (a cellular antioxidant) and brings about decarboxylation of auxin to impact development. This clearly indicates that the isomiR enables the plant to overcome the stressful conditions.

The isomiRs of osa-miR1423-3p and osa-miR530-5p miRNA were specifically absent in the salt stressed wild-type (SLPB) library. Similarly the isomiRs of osa-miR390-3p, osa-miR1436, osa-miR399j, osa-miR5161, osa-miR 7694-3p osa-miR6253 and osa-miR167a-3p were present only in GSL data sets. This indicates that they may be playing an important role in Gly induced stress tolerance. However the biological significance of the processing and function of the isomiRs requires detailed investigation

Identification of 3’ and 5’ End Modifications

It is well established that different post-processing editing mechanisms also modify the miRNA sequences by addition, substitution or deletion of nucleotides in a precursor-independent manner [97-99]. Though the functional importance of such modifications is yet to be explored but there are several reports which show that these modifications are playing a crucial role in stability of miRNAs and may also affect the miRNA:target interactions. Thus they may have a potential role in regulating plant stress responses. We have

used the miRMOD tool to investigate such modifications in all 357 known miRNAs identified in the Gly libraries [100].

79 miRNAs with 381end modifications and 83 miRNAs with 387 end modification were identified in GNL and GSL libraries, respectively miRNA. osa-miR156j-5p had maximum number of detected modifications (Supplementary Table ). It was observed that the 5’ end modifications are less common as compared to 3’modifications. The functional significance of the 5’ modifications are easier to comprehend as they directly affect the miRNA functions by changing the AGO sorting and target recognition (Fig. ). It was observed that the 5’ modifications majorly involved addition of “A” or “C” residues, while addition of A’s was the common modification at the 3′ end of miRNAs.

CONCLUSION

Glyoxalase (Gly) regulates the glutathione homeostasis and reduces the stress-induced cytoxicity of the cells by detoxifying MG. The implementation of sRNA sequencing and its analysis equips with an unrivalled opportunity to gain a detailed understanding of the global miRNA profiles. In this study, we observed that a greater number of miRNAs were up regulated by salt-stress in the Gly lines (GSL/GNL) and these included some of the known salt stress induced miRNAs reported in literature. This indicates that the deregulation in the miRNA expression profiles correlated well with the salt tolerant phenotypes of the Gly overexpressing plants. The targets for these miRNAs were predicted and their GO and KO analysis has provided insights in the role of miRNA regulatory networks in controlling the adaptation to salinity. The analysis also identified an enrichment of the isomiRs and modified forms (isoforms) of specific miRNAs. The target predictions suggested that miRNA modifications favour a change in the target transcripts such that the cellular machinery is geared to withstand stress. It needs to be deciphered how the MG and GSH levels initiate the molecular signaling cascades in plant adaptation to stress. This study necessitates further experimental validations to understand the mechanism operating in Gly plants which induce modulations in the miRNA biogenesis and function.

Table 1

Summary of reads obtained at each step of sequence analysis. Normal Glyoxalase (GNL) and salt stressed glyoxalase leaf (GSL) library.

Library	GNL	GSL
Read Length	36	36
Total No. of reads	26371105	26780310
Total number of High Quality (HQ) reads	25839193 (97.98%)	26129532 (97.57%)
Number of Primer/Adaptor contaminated HQ reads	671726	706017
%age of reads with non-ATGC bases	0.38	1.3
Total no. of unique Reads (18-24 length) after adapter trimming	2898648	3385190
Perfect match to indica genome	1804026 (62.24%)	2093148 (61.83%)

Table 2

List of predicted miRNAs with their calculated read values as transcripts per million.

___	miRNA	GNL	GSL	Description
Osa_1	AGAAGCTGTGGACTGTTTAGGGCA	16.078	15.87	Novel
Osa_2	TACCTCGTTAGATTCGTCTCA	7.356	8.812	Novel
Osa_3	AAATTCTGAGAAGCAGCTGGTTGG	4.815	7.057	Novel
Osa_4	AGAAGCTATGGACTGTTTGGGGCA	3.868	4.742	Novel
Osa_5	ACGTGACTTTGACTGCAGAGGATC	2.275	3.435	Novel
Osa_6	TGCCACGTGGACGATGACATG	0.569	1.195	Novel
Osa_7	ATCATGAGTAAATTGGTTCCCGGT	0.4550	1.008	Novel
Osa_8	TTGGATGGATGCCCGGGCAT	0.0758	0.037	Novel
Osa_9	AATGTATGACGCTGTTGACTTT	0.0758	0.224	Novel
Osa_10	AGGGATCCTCTGACTTAACTGCTG	1.0238	0.672	Novel
Osa_11	CGTGACATTCTCGGTGACATTGAG	0.1138	0.261	Novel
Osa_12	ACATTGTGAAACGGAGGAAGTAC	0.1138	0.112	Novel
Osa_13	TTACCAACCGGGACTAAAGATC	0.0758	0.112	miR2120 family
Osa_14	CAACTGACGGTAGTGGCATGAGA	0.0758	0	Novel
Osa_15	AAGGGATTTTTGGTGGATGTGACA	0.0379	0.187	Novel
Osa_16	ACGTAGATTACCGACAAAACCGA	0.0379	0.261	Novel
Osa_17	TCAAGGACCGTAGAATTGCTC	0.0379	0	Novel
Osa_18	TTGTACTAGGATGTGTCACAT	0	0.037	Novel
Osa_19	TACCTCGTTAGATTCGTCTCG	25.217	25.73	Novel
Osa_20	TGACATGTCGGAGTTGATGAGCAT	0.1138	0.224	Novel
Osa_21	TTTCGCTGGCGGATACTTAAGAGG	0.1138	0.075	Novel
Osa_22	ATCAGCAGTACTGCAGATAATCCC	0.0379	0	Novel
Osa_23	ACGATACCAAATAATTTTCGATCG	0	0.037	Novel
Osa_24	GTTCATTTAGTTGAATTAGAGTGG	0.1896	0.299	miR2863 family
Osa_25	ATTAGAATATGTTACATCCAACGG	0	0.037	miR5788 family
Osa_26	TGCTGAGCACAAACTAGAGAA	1.3651	2.651	3p arm of miR2866
Osa_27	GGCAGTCTCCTTGGCTAGC	0.1896	0.224	3-p arm of miR169g,h,j,l
Osa_28	AGGCAGTCTCCTTGGCTAGC	0.113761	0.149	3-p arm of miR169i,j,m
Osa_29	GGTGACTCCGGCCTCCTCGCC	0	0.112	3p of miR1846e
Osa_30	AGCTCCCTTCGATCCAATCC	0.3034	0.336	5-p arm of miR159f
Osa_31	TTGGTATTGTTTCGGCTCATG	0.3034	0.037	5p arm of miR171h