Comparative miRomics of Salt-Tolerant and Salt-Sensitive Rice.

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.

Introduction

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 [1], [2]. 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 [3]. In recent times, the disequilibrium in water cycle and poor drainage facilities in irrigated land has caused the problem of increase in soil salt [4]. The accumulation of excessive levels of salt has affected around 20 % of global irrigated land [5].

Soil salinity affects the plant’s ability to absorb water and also leads to high cellular concentrations of salts within the plant cells [6]. 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 [7]. 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 [8], [9], [10], [11], transcription factors [12], [13] and metabolic enzymes [14].

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 [15] 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 [16].

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 [17]. miRNAs are recognized as crucial effectors of diverse biological processes in plants [18], [19], [20], [21].

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 [22], [23]. 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 [24], [25], [26], [27], [28], [29]. 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 [30], [31]. 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 [32]. 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.

Architecture

Data Collection

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.

Small RNA Library Preparation

1 gram tissue was used for total RNA isolation as described earlier [33]. 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).

Library Preparation for RNAseq and Degradome Sequencing

1 gram tissue was used for total RNA isolation as described earlier [33]. 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 [34]. 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.

Workflow

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 [35]. 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 [36], [37] (http://www.mirbase.org/) using bowtie alignment tool [38] 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).

To observe the differential expression pattern of miRNAs a comparison across the normal and stressed libraries were performed and fold-change was calculated the values were plotted [39], [40].

N represents the fold-change

Identification of End Modifications and isomiRs of Known miRNAs

MiRMOD, a miRNA modification tool was used for the identification of end modified miRNAs [41]. 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.

Target Prediction

The miRNA target prediction was done using CleaveLand 4.4.3 tool [42]. 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 [43] (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) [44] (http://www.ricearray.org/).

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 [45]. Aligned transcripts were used for expression profiling in different datasets under various experimental conditions. The overall computation methodology is shown in Figure 1.

Figure 1:

A flow chart showing the stepwise analysis along with the computational methods used at each step. All tools were used with default parameters.

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

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).

Figure 2:

A flow chart showing the overview of analysis to highlight the scheme of data integration during the analysis.

Results

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.

Figure 3:

Length distribution of identified mature miRNAs in the leaf tissues of different rice varieties under normal and salt stressed conditions.

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.

Table 1:

Description of small RNA libraries and the number of miRNAs identified in each.

Abbreviation used	Tissue	Variety	Growth condition	Known miRNAs	Total no. of reads
PBNL	Leaf	Pusa Basmati	Normal	306	60,27,458
PBSL	Leaf	Pusa Basmati	Salt-stress	271	70,57,166
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
PKNL	Leaf	Pokkali	Normal	239	62,55,970
PKSL	Leaf	Pokkali	Salt-stress	299	3,23,25,958

Figure 4:

Venn diagram representing the distribution of miRNAs across the 6 leaf libraries under different experimental conditions. The miRNAs specifically expressed in the salt stressed libraries of Pusa Basmati (PBSL), Gly transgenics (GLSL) and Pokkali (PKSL) are indicated.

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).

Figure 5:

Distribution of salt-deregulated miRNAs (A) Log2 fold change as compared between the salt-stressed (SL) and normal (NL) leaf tissues across Pusa Basmati (PB), Pokkali (PK) and Glyoxalase overexpressing (GL) rice. (B) Comparative analysis of the fold changes in the deregulation patterns of the miRNAs in the salt tolerant rice as compared to salt sensitive rice. A white circle separates the two comparative sets. The lower half of the chart represents miRNA deregulations observed in PBSL/PBNL. The upper half of the chart represents deregulations of these miRNA in GLSL/GLNL (inner circle) and PKSL/PSNL (outer circle).

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.

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).

Table 2:

List of deregulated miRNA and mRNA.

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_Os10g35436.1	Up
osa-miR1425-5p	Down	LOC_Os10g35240.2	Up
osa-miR1425-5p	Down	LOC_Os10g35240.1	Up
osa-miR1425-5p	Down	LOC_Os10g35240.1	Up
osa-miR1425-5p	Down	LOC_Os08g01640.1	Down
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-miR160a,b,e,f-5p	Down	LOC_Os02g41800.2	Down
osa-miR160c,d-5p	Down	LOC_Os10g33940.1	Down
osa-miR160c,d-5p	Down	LOC_Os06g47150.1	Down
osa-miR160c,e-5p	Down	LOC_Os04g43910.1	Up
osa-miR160f-5p	Down	LOC_Os04g59430.1	Down
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-miR162b	Down	LOC_Os03g02970.1	Up
osa-miR164d	Up	LOC_Os12g41680.1	Up	No apical meristem protein, putative, expressed
osa-miR164d	Up	LOC_Os08g10080.1	Up
osa-miR166l-3p	Down	LOC_Os03g01890.2	Down	START domain containing protein, expressed
osa-miR166k-3p	Down	LOC_Os03g43930.1	Up
osa-miR166k-3p	Down	LOC_Os03g43930.2	Down
osa-miR166l-3p	Down	LOC_Os10g33960.1	Down
osa-miR166l-3p	Down	LOC_Os10g33960.2	Up
osa-miR166l-3p	Down	LOC_Os10g33960.3	Up
osa-miR166k-3p	Down	LOC_Os10g33960.4	Down
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-miR167a,c-5p	Down	LOC_Os04g57610.1	Down
osa-miR167a-5p,b	Down	LOC_Os06g46410.1	Down
osa-miR167c-5p	Down	LOC_Os06g46410.2	Down
osa-miR167a-5p	Down	LOC_Os07g33790.1	Down
osa-miR167a-5p,b,c-5p	Down	LOC_Os12g41950.1	Down
osa-miR171b	Down	LOC_Os02g44360.1	Up	Scarecrow, putative, expressed
osa-miR171f-3p	Down	LOC_Os02g44370.1	Up
osa-miR171e-3p	Down	LOC_Os04g46860.1	Down
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_Os05g03040.1	Down
osa-miR172a	Down	LOC_Os05g03040.2	Up
osa-miR172a	Down	LOC_Os07g13170.1	Down
osa-miR172d-3p	Down	LOC_Os07g13170.1	Down
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-miR172d-3p	Down	LOC_Os08g39630.1	Down
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-miR393a	Down	LOC_Os04g32460.1	NA
osa-miR396a-5p	Up	LOC_Os03g51970.1	Down	Growth regulating factor protein, putative, expressed
osa-miR396b-5p	Up	LOC_Os11g35030.1	Down
osa-miR396b-5p	Up	LOC_Os02g47280.1	Down
osa-miR396c-5p	Down	LOC_Os11g35030.2	Up
osa-miR396e-5p	Down	LOC_Os06g02560.3	Down
osa-miR396e-5p	Down	LOC_Os04g51190.3	Down
osa-miR396e-5p	Down	LOC_Os04g51190.1	NA
osa-miR396e-5p	Down	LOC_Os02g47280.2	Up
osa-miR396e-5p	Down	LOC_Os02g45570.1	Up
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-miR444c.1	Down	LOC_Os03g32499.1	Down	Expressed protein
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-miR444f	Down	LOC_Os08g06510.1	Up
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-miR812o-3p	Up	LOC_Os11g31690.2	Up	Expressed protein
osa-miR812o-3p	Up	LOC_Os09g32940.2	Up	Expressed protein
osa-miR812v	Down	LOC_Os04g47140.1	Up	Expressed protein
osa-miR812v	Down	LOC_Os02g07960.4	Up	STRUBBELIG-RECEPTOR FAMILY 3 precursor, putative, expressed

Figure 6:

Analysis of target transcripts regulated by the salt deregulated miRNAs A) Distribution pattern of the various transcription factors families targeted by the miRNAs (B–D) Gene Ontology (GO) enrichment analysis showing the target transcripts associated with the various biological process, (B) molecular function (C) and cellular component (D).

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 [46].

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.

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 [47]. 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).

Figure 7:

The distribution of isomiRs in the different libraries (A) A graphical representation of the distribution of isomiRs with 5′ end, 3′ end and mixed modifications (B) Analysis of isomiRs differing in the number of end modified nucleotides (C) Heat map of specific salt-deregulated isomiRs.

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 [48] and alter the targeted transcripts.

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 [49], [50]. 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 [49], [51]. Thus end modifications influence the putative function of the miRNAs.

Figure 8:

Frequency distribution of end modifications observed in 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.

Conclusions

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.

Click here for additional data file.