Genome-wide changes in microRNA expression during short and prolonged heat stress and recovery in contrasting rice cultivars

Highlight Heat susceptible and tolerant rice genotypes have different landscapes of miRNAs and the tolerant genotypes show efficient recovery mediated by miRNAs. Roots regulate miRNAs more sensitively than shoots during heat stress.


Introduction
Plants exposed to stress use multiple gene regulatory mechanisms, including post-transcriptional regulation of gene expression, to restore and re-establish cellular homeostasis. MicroRNAs (miRNAs) have emerged as ubiquitous critical regulator molecules in nearly all eukaryotes (Bartel, 2004). They are a group of small endogenous non-coding RNAs (19-25 nucleotides) that negatively modulate gene expression at the post-transcriptional levels by directing the cleavage of target genes or by inhibiting translation depending on the extent of the complementarity between the miRNA and its target (Reinhart et al., 2002;Chen, 2004; Jones-Rhoades and Bartel, 2004).
Plant miRNAs have diverse functions in growth and development processes, including root development (Guo et al., 2005), leaf morphogenesis (Kim et al., 2005), transition of growth from vegetative to reproductive stage (Achard et al., 2004;Lauter et al., 2005), and floral development and differentiation (Chen, 2004). The association of miRNAs with abiotic stress responses is now well established (Sunkar and Zhu, 2004;Jeong et al., 2011;May et al., 2013;Agarwal et al., 2015), but despite miRNAs being the key components of gene regulation during these stresses, the study of high temperature-responsive miRNAs is very limited in rice and other crops (Jeong and Green, 2013). Understanding miRNA expression and regulation during elevated temperature will help in better understanding the molecular pathways associated with high temperature response in the model monocot crop rice. The use of miRNAs for trait improvement has been successfully demonstrated wherein overexpression of OsmiR397 enhanced rice yield by increasing grain size and promoting panicle branching .
In this study, 16 small RNA libraries were prepared from root and shoot tissue of heat susceptible and tolerant rice genotypes. Two contrasting rice genotypes, Nagina 22 (N22) as tolerant and Vandana as susceptible, were selected based on earlier reports (Jagadish et al., 2008;Sailaja et al., 2014Sailaja et al., , 2015. The libraries represented small RNAs expressed during optimum temperature (control), high temperature (short and prolonged stress) and recovery. These libraries were sequenced to identify the miRNAs expressed in shoot and root of rice genotypes in control and stressed environments. The target genes of these miRNAs were predicted. Expression of miRNAs and target genes was verified using quantitative PCR.

Materials and methods
Plant materials and heat stress treatments Rice (Oryza sativa) seeds of cultivars N22 and Vandana were used for small RNA sequencing and miRNA analysis. The details of germination, growth conditions and heat stress treatments are the same as described previously (Sailaja et al., 2014). Germinated seedlings were transferred to Yoshida medium (Yoshida et al., 1976) and maintained at 13 h of light and 11 h dark at optimum day/night temperature (30 °C/24 °C). Eight-day-old seedlings were subjected to high temperature (42 °C/36 °C day/night) treatments for 24 h (short) and 5 days (long). To get the recovery samples, high temperaturetreated plants (for 4 days) were shifted back to optimum temperature for 24 h. In total, 16 samples were collected from root and shoot of N22 and Vandana representing control (30 °C/24 °C), short duration heat stress (SDS) treatment (42 °C/36 °C day/night for 24 h), long duration heat stress (LDS) treatment (42 °C/36 °C day/night for 5 days), and 24 h of recovery (REC) at optimum temperature (30 °C/24 °C) after a heat stress treatment (42 °C/36 °C) of 4 days. On the 13th day after germination, the 16 samples were harvested from four batches (control, SDS, LDS, and REC) including root and shoot tissues of N22 and Vandana.
Sequencing of small RNA Total RNAs were extracted from 16 samples using TRIzol Reagent (Life Technologies). RNA was isolated separately from root and shoot tissues subjected to control and stress treatments. Small RNA enriched samples were sequenced using Illumina sequencing. Samples of both the cultivars were harvested from two biological replicates. The small RNA quality statistics were generated using FASTX Toolkit (v 0.0.13) (http://hannonlab.cshl.edu/fastx_ toolkit/). The quality box plot graph was generated to visualize the quality score of each base with respect to its base position. Quality filtering was performed using the fastq quality filter program, and quality cutoff (q) of 20 and minimum percentage (p) value of 100 were assigned for preprocessing. Raw data were subjected to quality filtering (q20/p100) and the resultant reads were retained. During the quality filtering process, the low quality reads were discarded from the total reads. Adapter clipping was performed for those high quality (HQ) score reads using the FASTX toolkit. The sequence length distribution of HQ score reads was analysed individually using the FastQC tool (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/). The sequence data have been submitted to NCBI. The BioProject ID is PRJNA322758.

Analyses of miRNAs
The outline of the bioinformatics analysis work flow is given in Fig. 1. The O. sativa Nipponbare reference genome was downloaded from the MSU Rice Genome Annotation Project Database and Resource (Kawahara et al., 2013). The HQ score reads of 16 samples of O. sativa were mapped to the reference genome of O. sativa Nipponbare using bowtie (Langmead et al., 2009). The HQ score reads were mapped to miRNAs of O. sativa using bowtie (Langmead et al., 2009). To detect the miRNAs, the miRanalyzer standalone version (Hackenberg et al., 2011) was used. The number of mismatches allowed for alignment to known miRNAs, transcribed libraries and genome was set to zero. For several pairwise combinations, the differentially transcribed miRNAs with expression exhibiting two-fold or higher were chosen for downstream analysis. Gitools was used for drawing heat maps of differentially expressed miRNAs (Perez-Llamas and Lopez-Bigas, 2011).

Target prediction of miRNAs and function analysis
The target prediction of miRNAs was performed using the plant miRNA target prediction program psRNATarget (Dai and Zhao, 2011). To study Gene Ontology (GO), the best target genes were picked. Also the GO Slim term from the MSU Rice Genome Annotation Project Database and Resource (Kawahara et al., 2013) was downloaded for downstream analysis. The target genes were mapped to the GO Slim term using in-house perl script. Quantitative real-time PCR For quantification of miRNAs, small RNA was isolated using mir-Vana™ miRNA Isolation Kit (cat. no. AM1560, Ambion). The reverse transcription (RT) of miRNAs was performed using miScript II RT Kit (cat. no. 218161, Qiagen). The RT reaction was prepared by mixing and incubating the template small RNAs, miScript Reverse Transcriptase Mix, 10× miScript Nucleics Mix, and 5× miScript HiSpec Buffer, as specified by the manufacturer. Using miScript II RT Kit, mature miRNAs were polyadenylated by poly(A) polymerase and reverse transcribed into cDNA using oligo-dT primers. The oligo-dT primers have a 3′ degenerate anchor and a universal tag sequence on the 5′ end that allows amplification of mature miRNA during the real-time PCR step. cDNA was normalized and used for real time PCR by miScript SYBR Green PCR Kit (cat. no. 218073, Qiagen), which contains the miScript Universal Primer (reverse primer) and QuantiTect SYBR Green PCR Master Mix. The real time PCR reaction was prepared by mixing 2× QuantiTect SYBR Green PCR Master Mix, 10× miScript Universal Reverse Primer, miRNA specific forward primer (Supplementary Table S1 at JXB online), and template cDNA, as specified by the kit manufacturer. U6 was used as internal control for miRNA quantification (Ding et al., 2011). The quantitative real-time PCR (qRT-PCR) temperature profile was followed as described previously (Sailaja et al., 2014).
For quantification of target genes, mRNA was isolated from the same plant samples using mirVana™ miRNA Isolation Kit as the kit provides an option for isolating small RNAs and long RNAs from the same preparation by adjusting the ethanol concentration. cDNA synthesis of mRNAs was performed using the Improm-II reverse transcription system, as specified by the manufacturer (Promega, cat. no. A3800). For cDNA synthesis, oligo dT primers were used. cDNA was treated with RNase and normalized for equal concentration. qRT-PCR was performed using SYBR Premix Ex-Taq (Takara, cat. no. RR820A) in an ABI GeneAmp 7500 Sequence Detection System. OsActin1 (Lee et al., 2011) was used as the internal control for quantification of genes. Forward and reverse primers used for qRT-PCR of target genes are listed in Supplementary Table  S1. Amplification conditions were followed as described previously (Sailaja et al., 2014).
With three biological replications, each reaction was run in duplicate and the melting curves were constructed using Dissociation Curves Software (Applied Biosystems) to ensure that only a single specific product was amplified. The comparative threshold cycle (C T ) method was used to quantify the relative expression levels of miRNAs and genes in real-time PCR. ΔC T was calculated by C T target-C T reference (internal control). Further ΔΔC T values were calculated by ΔC T of the stress sample − ΔC T of the control sample, and then the fold difference was calculated from 2 −ΔΔCT . Similarly, the ΔC T standard deviation was calculated as given at http://www. appliedbiosystems.com/absite/us/en/home/support/tutorials/realtime-pcr-trouble-shooting-guide.html.

Results
We used the massively parallel high throughput sequencer to investigate the genome-wide identification and expression profiles of miRNAs in Vandana and N22 rice cultivars, particularly for the heat stress-responsive miRNAs. Sixteen small RNA libraries were constructed by the use of total RNAs isolated from control shoot (CS), control root (CR), short duration stress (SDS) shoot (SS), SDS root (SR), long duration stress (LDS) shoot (LS), LDS root (LR), recovery shoot (RS), and recovery root (RR) tissues of both the cultivars. Small-RNA sequencing of 16 libraries yielded a total of 319 260 462 raw reads including of 271 353 814 high-quality raw sequences ( Table 1). The high-quality (HQ) sequence reads of 16 samples of O. sativa extracted from the shoot and root in different conditions (control, stress, prolong stress and recovery) were mapped to the reference genome of O. sativa Nipponbare using bowtie (Langmead et al., 2009) (Table 1). Further, the HQ score reads were mapped to miRNAs of O. sativa using bowtie (Langmead et al., 2009). As shown in Fig. 2, the highest abundance of miRNAs was found in the shoot (control) library of Vandana. In the N22 libraries, RS showed the maximum number of miRNA reads. The length distribution of sRNA reads revealed that heat stress treatments affect small RNA metabolism extensively. The sRNA class of 24 nt was the most abundant group in the control sample of each library followed by recovery samples ( Supplementary Fig. S1). Stressed tissues showed variation in terms of sequence length.

Genome-wide expression analysis of heat-responsive miRNAs in Oryza sativa
MicroRNAs identified in each library were compared in 36 pairwise combinations to decipher differentially expressed miRNAs. The expression value of each rice miRNA was analysed and miRNAs showing >2.0-fold expression are listed (Supplementary Table S5). The target prediction of differentially expressed miRNAs was performed and the list of target genes is provided (Supplementary Table S6). Heat-treated shoot tissue of N22 was compared with control shoot. Here, osa-miR1423a-5p, osa-miR1427, osa-miR2055, osa-miR1863a, and osa-miR5072 showed up-regulation, while osa-miR166n-5p, osa-miR2863b, and osa-miR396f-3p were down-regulated after SDS, LDS, and REC. The common and differentially expressed miRNAs of N22 shoot after heat stress treatments and recovery are shown in Fig. 3. In a similar way, heat-treated root tissue of N22 was compared with control root. Here, osa-miR1427 was up-regulated while osa-miR1878 was down-regulated after SDS, LDS, and REC. Interestingly, the majority of the miRNAs showed down-regulation in root of N22 after LDS and REC. MicroRNAs of N22 root showing common and differential expression after SDS, LDS, and REC are shown in Fig. 3.
High temperature-treated shoot tissue of the susceptible cultivar Vandana was also compared with control shoot. Here, osa-miR1879, osa-miR394, osa-miR3979-3p, osa-miR408-5p, osa-miR444f, osa-miR531a, osa-miR440, and osa-miR444a-5p showed up-regulation after SDS, LDS, and REC. The common and differentially expressed miRNAs of Vandana shoot after heat stress treatments and recovery are shown in Fig. 4. Likewise, heat-treated root tissue of Vandana was compared with control root. Here, osa-miR396c-5p, osa-miR5072, osa-miR5082, osa-miR528-5p, and osa-miR169i-5p were up-regulated, while osa-miR1878 showed down-regulation after SDS, LDS, and REC. MicroRNAs of Vandana root showing common and differential expression after SDS, LDS, and REC are shown in Fig. 4. The majority of miRNAs showed up-regulation in root and shoot of Vandana after SDS, LDS, and REC. In addition to differentially expressed miRNAs, there were many miRNAs that showed expression preferentially in one tissue or treatment and were absent in the other; these were also identified (Supplementary Table S7).

Expression validation of miRNAs and target genes
In order to validate the deep sequencing data, expression analysis of 17 miRNAs was performed in 22 different combinations of heat treatment versus control and N22 versus Vandana (Fig. 6). The results showed a similar trend of expression of miRNAs in NGS and qPCR, though the degree of expression varied. Further, expression correlation of miR-NAs and their target genes was analysed and suggested a negative correlation between miRNAs and gene targets (Fig. 7).

Function analysis
For each library, target transcripts representing genes with a known function were categorized into biological process, cellular component and molecular function according to the  Table S8). The target transcripts of miRNAs in the molecular function category were related to binding, protein binding, catalytic activity, hydrolase, kinase, transferase, nuclease, transporter, signal transducer, enzyme regulator activity, motor activity, receptor activity, and structural molecule activity. Most of the miRNA target genes were assigned to the binding category whose present sequences appear to be involved in nucleic acid binding, protein binding and ion binding. Since these target sequences encode transcription factors, it is in accordance with previous reports that a large proportion of the miRNA targets encode transcription factors.
The number of miRNA targets falling under different categories of molecular function was also analysed for different libraries. Significant difference was not observed in the shoot library of Vandana. However, root tissue of Vandana showed a distinct difference in the number of miRNA targets falling under catalytic activity, hydrolase activity, kinase activity, nucleotide binding, protein binding, and RNA binding. In particular, 11 miRNA targets showed kinase activity in control and stress root tissue and only three in recovery root tissue. Similarly, 12, 13, and 5 miRNA targets were identified to play a role in nucleotide binding in control, stress, and recovery root tissue, respectively. Similar to Vandana, shoot libraries of N22 also did not show much difference in number of miRNA targets under different molecular function categories; however, a significant difference was observed in root of N22. A distinct difference in root libraries was observed in the number of miRNA targets associated with binding, catalytic activity, hydrolase activity, kinase activity, nucleotide binding, and protein binding; 10, 13, and 9 miRNA targets showed kinase activity in control, stress, and recovery root of N22, respectively.
In the biological processes category, the total number of miRNA targets in each library fell into multiple classes. However, it is notable that some of these transcripts in the biological process category were related to stress response. In the shoot library of Vandana, difference in miRNA targets falling in different categories was not significant when compared among control, stress, and recovery. However, the difference was more evident in root tissue of Vandana. Here, 23, 24, and 15 target transcripts were found in control, stress, and recovery tissue, respectively. A significant difference was also observed in miRNA targets associated with biosynthetic processes, metabolic processes, protein modification processes, signal transduction, and transport. In particular, the target transcripts associated with protein modification processes were 14, 13, and 4, while metabolic processes showed 30, 25, and 19 transcripts in control, stress, and recovery root, respectively.
Similar to Vandana, shoot tissue of N22 did not show much variation in terms of number of miRNA targets falling under different categories of biological processes. However, root tissue of N22 showed significant variation. Control, stress, and recovery root tissue showed, respectively, 27, 37, and 24 target transcripts under biological processes; 23, 27, and 17 target transcripts under biosynthetic processes; 33, 37, and 27 target transcripts under cellular processes; 30, 37, and 27 target transcripts under metabolic processes; 12, 16, and 11 target transcripts under protein modification processes; 12, 20, and 12 target transcripts under stress response; 4, 8, and 4 target  transcripts under response to endogenous stimuli; and 9, 11, and 7 target transcripts under signal transduction.
Cellular components of miRNA targets were also analysed for control, stress, and recovery libraries of N22 and Vandana. In terms of number of miRNA targets showing different cellular locations, Vandana shoot tissue did not show much variation in different libraries. However, root tissue showed variation under cellular component: cytosol, intracellular, membrane, mitochondrion, nucleus, and plasma membrane. In particular transcripts localizing in plasma membrane were 11 in control, 12 in stress and 5 in recovery of Vandana root. N22 root showed 13, 15, and 10 transcripts localizing in plasma membrane in control, stress, and recovery, respectively. In N22 shoot, a significant difference in miRNA targets localizing in nucleus was observed. The number of transcripts localizing in nucleus was 16 (control), 11 (stress), and 18 (recovery). In general, target genes localizing in different cellular components were lesser in stress as compared with control. In comparison with shoot, root tissue of N22 showed a higher number of transcripts localizing in different cellular components under stress condition (Supplementary Table S8).

Discussion
This study is specifically important in addressing three major issues of heat stress response which have not been previously reported in miRNA/transcriptome studies on high temperature stress. First, we have studied the miRNA regulation during short-as well as long-duration heat stress. Second, while treating the plants with high temperature, the day/night difference of temperature was maintained. Third, microRNA regulation after recovery was also studied.
In all, 162 miRNA families were identified in two rice genotypes, which included monocot preferential miRNA families-miR437, miR444, miR528, miR1318, miR1432, and miR1436 Pandey et al., 2014). Similar to other reports of deep sequencing of microRNAs in rice, the expression of conserved miRNA families was significantly higher than that of non-conserved miRNAs (Sunkar et al., 2008;Zhu et al., 2008;Li et al., 2011;Chen et al., 2013). MicroRNA166 was most abundantly present across the 16 rice libraries sequenced in this study. Based on this study, we classified miRNA families into preferential expression in cultivars (N22 and Vandana), tissue (shoot and root), and treatments (control and heat stressed). However, this classification needs to be further validated as small RNA sequencing may not be saturated to cover all the small RNAs in a given sample. The majority of miRNA families were common. N22 preferential miRNA families included osa-miR1427, osa-miR1848, osa-miR1865, osa-miR2872, osa-miR5155, and osa-miR5484, while Vandana preferential miRNA families were osa-miR414 and osa-miR531. These miRNAs were suggested as potential candidates for selection in domestication, novel in cultivated rice when compared with wild rice . In addition, microRNA families osa-miR2120, osa-miR1846, osa-miR5079, osa-miR1318, osa-miR1863, osa-miR2275, osa-miR5072, osa-miR5160, osa-miR5162, osa-miR5505, and osa-miR5513 were suggested to be involved in cultivated rice domestication and were present in both Vandana and N22. While a previous study by Wang et al., (2012) was carried out with wild rice (Oryza rufipogon) to determine the role of miRNA genes in rice domestication, this study supports those findings for miRNAs identified in cultivated rice (Oryza sativa). O. sativa has been arranged into five distinct groups, corresponding to indica, aus, aromatic, temperate japonica, and tropical japonica rices (Garris et al., 2005), and this study provided an opportunity to investigate miRNAs of the indica (Vandana) and aus (N22) groups.
Expression of osa-miR1436 was recorded in both genotypes but preferentially in stressed tissues. It was found putatively to target a large number of genes such as CRAL/TRIO domain-containing protein, scramblase, transmembrane protein, tyrosine protein kinase domain-containing protein, glycoprotein 3-α-L-fucosyltransferase A, cytochrome P450, serine/threonine-protein kinase, WRKY34, etc., suggesting its general response to stress. This miRNA was shown to be associated with arsenic stress response in Brassica juncea (Srivastava et al., 2013). The expression of osa-miR1436 was decreased while the expression of its target gene was increased upon glyphosate treatments (Unver et al., 2010). Sunkar et al., (2008) suggested that miR1436 may be involved in drought stress in rice. Osa-miR5076 was also identified preferentially in the stressed library of both the genotypes. The predicted target of osa-miR5076 is ferredoxin-dependent glutamate synthase involved in nitrogen assimilation in rice. High temperature affects nitrogen metabolism of plants significantly (Xu and Zhou, 2006).
The miRNAs expression profile was further evaluated with a primary focus on the heat-tolerant N22 genotype of rice. In comparison with control, the SDS shoot of N22 showed upregulation of osa-miR5788, osa-miR1427, and osa-miR2055, and down-regulation of osa-miR166n-5p, osa-miR1425, and osa-miR399. The target of osa-miR5788 is U-box domain containing heat shock protein. Overexpression of AtCHIP, a U-box protein, in Arabidopsis rendered plants more sensitive to both low-and high-temperature treatments (Yan et al., 2003). miR399 controls phosphorus homeostasis by regulating the expression of a ubiquitin-conjugating E2 enzyme in Arabidopsis (Chiou et al., 2006) and participates in the regulation of multiple nutrient starvation responses (Hu et al., 2015).
The up-regulated microRNAs in N22 shoot after LDS (in comparison with control) were osa-miR5083, osa-miR5504, osa-miR5160, osa-miR5072, osa-miR2055, and osa-miR1427, and the down-regulated miRNAs were osa-miR166, osa-miR5073, osa-miR156, and osa-miR390. The miR156 targets squamosa promoter binding protein-like (SPL) and dihydroflavonol-4-reductase (DFR) pathways to coordinate the relationship between development and abiotic stress tolerance in plants (Cui et al., 2014). Recently it was shown that miR156 isoforms are highly induced after heat stress and are functionally important for heat stress memory. miR156 promotes sustained expression of stress-responsive genes and is critical after high temperature stress adaptation to recurring heat stress in Arabidopsis (Stief et al., 2014). This was not differentially expressed in SDS but appeared preferentially in LDS. miR390 was down-regulated during cadmium stress in rice (Ding et al., 2011). It contributes to coordinate multiple phytohormone signals for the regulation of cell growth and senescence in plants (Curaba et al., 2014). This miRNA targets STRUBBELIG-RECEPTOR FAMILY (SRF) 6 precursor, which was strongly induced in plants exposed for a prolonged time to heat stress (Eyüboglu et al., 2007). Our results indicate that miR156 and miR390 may be involved in acclimation or adaptation to long duration heat stress in rice, which needs to be confirmed in further studies.
The up-regulated miRNAs in shoot of N22 after REC (in comparison with control) were osa-miR1863a, osa-miR1427, osa-miR2055, osa-miR168a-3p, osa-miR167a-3p, osa-miR5083, osa-miR5160, and osa-miR5072. miR1863 is required for silencing heterochromatin by methylation in rice (Wu et al., 2010). Recent studies show the involvement of the DNA methylation process in high temperature response in plants (Sanchez and Paszkowski, 2014;Naydenov et al., 2015). Imprinting of genes may help in recovery of stressed plants after a heat stress episode. osa-miR168a-3p was shown to be an arsenic-responsive miRNA in the shoots of rice (Yu et al., 2012). Further, miR168a-3p and miR167a-3p were identified as heat stress (37 °C for 6 h)-responsive microRNAs in Arabidopsis ( Barciszewska-Pacak et al., 2015). osa-miR5083, osa-miR5160, and osa-miR5072 were up-regulated after LDS and REC indicating their significant role in response and recovery from long duration high temperature stress in rice. It is interesting to note that osa-miR2055 and osa-miR1427 were up-regulated after SDS, LDS, and REC of N22 suggesting a common function of these miRNAs in regulating tolerance to high temperature. The miRNAs showing downregulation after REC were osa-miR528-5p, osa-miR528, osa-miR166, osa-miR5073, osa-miR2863, osa-miR6246, and osa-miR1863c. miR528 was down-regulated in rice seedlings treated with hydrogen peroxide (H 2 O 2 ) to induce oxidative stress . osa-miR166 was down-regulated after SDS, LDS, and REC while osa-miR5073 showed down-regulation only after LDS and REC of N22 shoot indicating some miRNAs are recruited only when the duration of heat stress is long, and they continue their role during recovery too.
In comparison with control, the root tissue of N22 showed up-regulation of osa-miR5504, osa-miR169, osa-miR1427, osa-miR390, osa-miR166b-5p, osa-miR5082, and osa-miR319a-3p after SDS. Interestingly, osa-miR1427 showed up-regulation in shoot also after SDS, LDS, and REC. In a recent study, osa-miR169i-5p.2, osa-miR390-3p, and osa-miR5082 showed high abundance in root tips suggesting their important role in root development of rice (Ma et al., 2013). The up-regulation of osa-miR169, osa-miR390, and osa-miR5082 in N22 root after SDS suggests that these miRNAs may be involved in root development of N22 under high temperature stress. osa-miR319 was suggested to be an oxidative stress-responsive miRNA family in rice . We observed the down-regulation of osa-miR2275d and osa-miR1850.1 in N22 root after SDS. osa-miR2275 was shown to be regulated by drought and cold stress in rice (Barrera-Figueroa et al., 2012). It is important to note that osa-miR2275 is involved in regulation of biogenesis of secondary small interfering RNAs (siRNAs) of either 21 or 24 nucleotides in a phased manner (Song et al., 2012), which may have a significant function in abiotic stress response.
The significantly up-regulated microRNAs in N22 root after LDS (in comparison with control) were osa-miR5504, osa-miR1427, osa-miR5083, osa-miR169e, osa-miR156, and osa-miR6250. The induced expression of osa-miR5504 and osa-miR1427 was observed in roots after SDS also. miR6250 showed high abundance in rice roots in response to arsenite treatment (Liu, 2012). Significantly down-regulated miRNAs were osa-miR169a, osa-miR1863, osa-miR2878, osa-miR1425, osa-miR1876, osa-miR171, osa-miR1874, and osa-miR159. Interestingly, osa-miR169e was up-regulated while osa-miR169a was down-regulated in root of N22 after LDS. Tissue differential expression of members in the same family of miRNAs has been reported previously (Jeong et al., 2011). Numerous miRNAs were down-regulated in root tissue of N22 after REC suggesting the up-regulation of target genes. Expression analysis of 13 genes and 9 miRNAs suggested root as more responsive in N22 after recovery from heat stress treatments (Sailaja et al., 2014).
The miRNA composition and expression pattern in the heat susceptible cultivar Vandana was largely different form heat tolerant N22 when analysed after SDS, LDS, and REC (Fig. 4, and Tables S4 and S5). The miRNAs expressed preferentially in the tolerant rice genotype N22 at high temperature were osa-miR1439, osa-miR1848, osa-miR2096, osa-miR2106, osa-miR2875, osa-miR3981, osa-miR5079, osa-miR5151, osa-miR5484, osa-miR5792, and osa-miR5812. These miR-NAs may have a more specific role in heat stress tolerance. MicroRNA1439 was suggested to be involved in the homeostasis of metal ions and plays a role in plant response to salt stress . Recently, osa-miR1848 was shown to target the obtusifoliol 14α-demethylase gene (OsCYP51G3) and mediates the biosynthesis of phytosterols and brassinosteroids during development and stress response (Xia et al., 2015). osa-miR2106 was found to be up-regulated in arsenic stress in rice (Yu et al., 2012). Beside the specific miRNAs expressed preferentially in heat tolerant N22 as discussed above, there were a number of miRNAs that showed differential expression when compared with Vandana. Higher expression of osa-miR1859, osa-miR5504, and osa-miR5513 was observed in N22 than in Vandana shoot after SDS, LDS, and REC. osa-miR408-3p was down-regulated after SDS, LDS, and REC of N22 shoot when compared with corresponding samples of Vandana shoot. Interestingly, osa-miR408 was shown to be up-regulated in N22 and Vandana during drought stress in young seedlings, as well as in flag leaf and spikelets of adult plants (Mutum et al., 2013). Both these genotypes are drought tolerant, however, contradicting their response to high temperature. It would be interesting to analyse the differential role of the osa-miR408 family in drought and heat stress response in rice. miR408 is suggested to target the plastocyanin-like protein family and helps in maintaining the cellular redox state (Mutum et al., 2013).
It was interesting to observe that a larger number of miR-NAs were up-regulated in N22 SR than in Vandana SR suggesting that root was more responsive in N22. We had made a similar observation in an earlier study with a limited number of miRNAs and transcripts (Sailaja et al., 2014). The role of osa-miR5150-3p, osa-miR528-5p, osa-miR398b, and osa-miR408-3p in recovery response of the heat tolerant genotype needs to be further characterized as these miRNAs showed differential expression in N22 RR when compared with Vandana RR. Overall, this study suggests that N22 has distinct miRNA-mediated gene regulation machinery in comparison with Vandana after various heat stress treatments.
The functional analysis of target transcripts of miRNAs identified in different rice libraries suggests that root tissue is more responsive under heat stress treatments and recovery. The molecular function and biological process category of miRNA targets of N22 and Vandana showed significant difference in various GO categories after SDS, LDS, and REC in root, while shoot did not show much difference. Further, functional analysis of target transcripts suggested that N22 root has a more efficient recovery process as compared with Vandana. Our previous study based on expression analysis of a set of genes and miRNAs indicated that root is more responsive than shoot during heat stress in rice and N22 has a more efficient recovery mechanism to cope with high temperature stress than Vandana (Sailaja et al., 2014).
In conclusion, we identified miRNAs from 16 rice samples, their target genes and molecular pathways that are regulated by elevated temperature in rice. This is the first report of whole-genome miRNAs expression at a wider scale that includes profiling of miRNA expressed in root and shoot of heat tolerant and susceptible rice cultivars during control, short and long durations of heat stress treatment and also recovery after LDS. Further, unlike other reports where heat stress was given for a certain duration, generally a few hours, and then returned to optimal temperature, we attempted to map the miRNA expression by maintaining the day/night temperature difference in control and elevated temperature treatments, which is a more realistic situation. These results have opened a baseline for understanding miRNA-mediated regulation of rice response to high temperature. Even though we see contrasting response to heat stress in shoot growth of N22 and Vandana, the miRNA expression shows more differential activity in response to heat in root than in shoot. This study suggests that heat susceptible and tolerant rice cultivars have a different landscape of miRNA expression and regulation at elevated temperature, which might be associated with the molecular events leading to contrasting responses of rice genotypes to heat stress. The study provides the first comprehensive and specific miRNA profile of heat-induced miRNAs in shoot and root of rice. The identification of genotype and tissue preferential miRNAs will help in widening the knowledge base of miRNA-mediated gene regulation in rice.

Supplementary data
Supplementary data are available at JXB online. Fig. S1. The length distribution of small RNA reads. Table S1. List of miRNAs and target gene primers used in this study. Table S2. List of miRNA families identified in different rice libraries. Table S3. Predicted target genes of miRNAs expressed preferentially in control or stressed tissue of N22 and Vandana. Table S4. Mature miRNAs along with details of read sequence and read count. Table S5. List of differentially expressed miRNAs in 36 pairwise combinations. Table S6. List of candidate target genes of differentially expressed miRNAs. Table S7. List of miRNAs showing expression preferentially in one tissue or treatment and absent in others. Table S8. Categorization of miRNA target transcripts into biological process, cellular component and molecular function according to the ontological definitions of the GO terms.