Comparative analysis of microRNA profiles of rice anthers between cool-sensitive and cool-tolerant cultivars under cool-temperature stress

sterility. As a first step toward understanding the molecular mechanisms of cool tolerance in developing anthers of rice, we report here a comprehensive comparative analysis of microRNAs between cool-sensitive Sasanishiki and cool-tolerant Hitomebore cultivars. High-throughput Illumina sequencing revealed 241 known and 46 novel microRNAs. Interestingly, 15 of these microRNAs accumulated differentially in the two cultivars at the uninucleate microspore stage under cool conditions. Inverse correlations between expression patterns of microRNAs and their target genes were confirmed by quantitative RT-PCR analysis, and cleavage sites of some of the target genes were determined by 5’ RNA ligase-mediated RACE experiments. Thus, our data are useful resources to elucidate microRNA-mediated mechanism(s) of cool tolerance in rice anthers at the booting stage.


INTRODUCTION
The precise development of pollen grains in anthers is achieved by cooperative regulation between gametophyte and sporophyte cells. If such finely balanced cooperation is disorganized under stress, male sterility occurs in plant reproduction (Oda et al., 2010). Male gametophyte development in Oryza sativa L. is susceptible to cool-temperature stress (Sakai, 1937(Sakai, , 1943Satake and Hayase, 1970). Cool-temperature injury at the booting stage inflicts severe damage to rice pollen. Abnormal tapetum hypertrophy at the tetrad to early uninucleate microspore (UNM) stages is a major cause of pollen sterility, leading to a reduction of the rice grain yield (Sakai, 1949). Cool temperatures in the summer of 1993 in the Tohoku region of Japan caused heavy damage to Japanese rice production: total rice grain yield decreased about 44% in the region, and about 33% in Japan, compared to the average yield for the previous 30 years (MAFF, 2012). Sakai (1949) observed differences in the occurrence of abnormal tapetum between rice cultivars, and showed that cooltolerant rice cultivars at the booting stage have less abnormal tapetum structure. We have studied two japonica rice cultivars, Hitomebore and Sasanishiki, cultivated in the Tohoku region (Oda et al., 2010). These two cultivars are genetically not distant; however, their response to cool-temperature stress at the booting stage is clearly different. In our previous study, under cooltemperature growth conditions at the booting stage, pollen sterility following incomplete tapetum degradation was observed in Sasanishiki, while anthers from Hitomebore were able to produce normal pollen grains (Oda et al., 2010).
Genetic systems for conferring cool tolerance in rice are not clear, although many quantitative trait loci (QTLs) for cool tolerance at the booting stage have been detected (Takeuchi et al., 2001;Andaya and Mackill, 2003;Ye et al., 2010;Zhou et al., 2010;Kuroi et al., 2011). Saito and his colleagues (2001Saito and his colleagues ( , 2004Saito and his colleagues ( , 2010 have explored QTLs that correlated with cold tolerance of spikelet fertility at the booting stage in rice, and they indicated that an F-box protein gene, located on chromosome 4, is one of the genes conferring cold tolerance based on results from map-based cloning. For other QTLs, however, identification of gene(s) has not been completed.
Abiotic stress from environmental changes can severely damage plants even if the stress is short-term. Rapid regulation of transcript level by small RNA is therefore a useful response to the unexpected and unavoidable consequences of abiotic stress. MicroRNAs (miRNAs) are endogenous, approximately 21-24 nucleotides (nt) long, non-coding RNAs belonging to the family of small RNAs, and are thought to play important roles in developmental processes and stress responses by regulating target mRNA expression post-transcriptionally in plants (Bartel, 2004;Kidner and Martienssen, 2005). The primary miRNA (pri-miRNA) that is transcribed by polymerase II from an miRNA gene locus is processed via precursor miRNA into mature miRNA/miRNA* duplex, which has a stem-loop structure, by Dicer-like (DCL) 1 in plants ( Kurihara and Watanabe, 2004). The miRNA that is incorporated into RISC (RNA-induced silencing complex) recognizes specific mRNA targets based on fully or partly complementary sequence and guides the RISC to repress target transcripts by translational inhibition and silencing of the target mRNAs (Voinnet, 2009). It has been shown that plants undergoing abiotic stresses regulate the expression of specific genes by miRNAs (Abdel-Ghany and Pilon, 2008;Jeong et al., 2009;Vidal et al., 2010); such regulation was also observed in rice shoots at the prophyll emergence stage (Lv et al., 2010) and in 7-dayold rice seedlings (Yang et al., 2013) under a 4 °C treatment.
Although previous papers have suggested a relationship between miRNA profiles and rice anther or pollen development (Wei et al., 2011;Peng et al., 2012;Yan et al., 2015), further investigation is needed to understand the complete picture of the molecular genetic system(s) active in this tissue, and in particular the role of miRNAs in developing rice anthers under cool stress. Here, we focused on miRNA accumulation in immature anthers of rice under cool stress related to cool-temperature injury at the booting stage, and compared the miRNA profiles between cool-sensitive Sasanishiki and cool-tolerant Hitomebore cultivars. Comprehensive analysis of anther miRNA using next-generation sequencing (NGS) showed that some miRNAs are differentially accumulated in UNM anthers under cool stress in the cool-tolerant cultivar. The data obtained in this study are useful resources to understand microRNA-mediated mechanism(s) of cool tolerance in rice anthers at the booting stage.

Plant materials
Two japonica rice cultivars, Sasanishiki and Hitomebore, were grown in an experimental field at Miyagi Prefectural Furukawa Agricultural Experimental Station, Japan under normal growth conditions. For the cool growth condition, cool-temperature stress was applied using a deep water irrigation system with cool water (19 °C, 25 cm deep) (Matsunaga, 2005) in the experimental field, according to Oda et al. (2010). Anthers at the UNM and tricellular pollen (TCP) stages were collected and stored at −80 °C until analysis.
Illumina sequencing Total RNA was extracted from rice anthers using the mirVana miRNA Isolation Kit (Ambion, Austin, TX, USA). Each cDNA library was constructed using the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA, USA) with a standard protocol, and single-read multiplex sequencing (read length: 100 bp) was conducted by Illumina Hiseq2000 sequencer (Illumina).
Among uncategorized sequences, if they were mapped to less than 10 loci in the reference genome, upstream and downstream sequences of the corresponding genomic positions were extracted (the extracted windows were shifted or increased per 10 nt from 50 to 150 bp upstream or downstream), and their secondary structures were predicted using RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/ RNAfold.cgi; Hofacker, 2003) with default parameters. Extracted segments were considered as novel miRNAs (nov-miRNAs) if their secondary structure indicated the following characteristics: (1) minimum free energy is less than -18 kcal/mol; (2) the total number of mismatches and gaps between putative miRNA and miRNA* is four or less, but the number of gaps is one or less (Meyers et al., 2008); (3) putative miRNA and miRNA* sequences are in the same library and miRNA counts exceed miRNA* counts; (4) the normalized read count value of a putative miRNA (reads per million, RPM) is over 10 in the library.
Quantitative reverse transcription (qRT)-PCR To validate the Illumina sequencing of rice miRNAs, expression levels of miRNAs or target genes were investigated by qRT-PCR. For the target genes, total RNA was extracted from anthers at each developmental stage using the mirVana miRNA Isolation Kit, followed by firststrand cDNA synthesis by Superscript III (Invitrogen, Carlsbad, CA, USA). Small RNA was extracted from the total RNA according to the small RNA isolation protocol (Invitrogen), followed by cDNA synthesis using the Mir-X miRNA first-strand synthesis kit (Clontech, Palo Alto, CA, USA). The reaction mixture (10 μl) consisted of 1×mRQ Buffer, 1.25 μl mRQ Enzyme and 3.75 μl RNA sample, and the mixture was incubated for 1 h at 37 °C and then terminated at 85 °C for 5 min to inactivate the enzyme. Reaction products were diluted with 90 μl ddH 2 O. qRT-PCR was performed in an MJ Mini Thermal Cycler and MiniOpticon Real-Time PCR System (Bio-Rad Laboratories, Hercules, CA, USA) using SYBR Advantage qPCR Premix (Clontech). The PCR reaction mixture (25 μl) consisted of 0.2 μM forward and reverse primers, 1×SYBR Advantage qPCR Premix and about 100 ng cDNA, and the PCR cycle parameters were 95 °C for 30 sec followed by 45 cycles of 95 °C for 5 sec and 60 °C for 10 sec. The miRNA forward primers were designed based on the full miRNA sequences, and the mRQ 3' primer in the Mir-X miRNA first-strand synthesis kit was used as the reverse primer. U6 snRNA was used as the reference gene for miRNA expression, and eukaryotic translation elongation factor 1a (eEF-1a), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), ubiquitin 5 (UBQ5) and U6 snRNA were used as the references for target gene expression. Primers used for qRT-PCR experiments are listed in Supplementary Table S1. The ΔΔCt method was adopted to calculate relative expression value.

5' RNA ligase-mediated (RLM) RACE
To investigate cleavage sites of predicted target genes, 5' RLM-RACE was conducted for representative target genes. Approximately 250 ng of mRNA extracted from anthers of Hitomebore at the UNM stage under cool conditions was ligated to a 5' RACE Adapter, without calf intestine alkaline phosphatase treatment, using a FirstChoice RLM-RACE Kit (Ambion), followed by reverse transcription according to the kit protocol. One microliter of reverse transcription product was used as a template for outer 5' RLM-RACE PCR using the 5' RACE Outer primer supplied in the kit and a gene-specific outer primer designed for each putative target gene (Supplementary Table  S2). Next, 1 μl Outer PCR product was used as a template for inner 5' RLM-RACE PCR using the 5' RACE Inner primer supplied in the kit and a gene-specific inner primer designed for each putative target gene-specific upstream site of the outer gene-specific primer-binding site (Supplementary Table S2). The electrophoresed inner PCR products were purified, and cloned into the pTAKN-2 vector (BioDynamics Laboratory, Tokyo, Japan). More than six independent insert-positive clones were sequenced with a 3500 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) and analyzed by Sequencing Analysis software v.5.4 (Applied Biosystems).

RESULTS AND DISCUSSION
Characteristics of small RNA sequences in rice anthers Illumina high-throughput sequencing was conducted on eight libraries of small RNAs in rice anthers collected from each of two cultivars (Sasanishiki or Hitomebore) at the UNM or TCP stage under normal or cool growth conditions (Table 1). We selected these eight samples to compare small RNA profiles between the coolsensitive Sasanishiki and cool-tolerant Hitomebore in immature anthers under cool stress. After removing low-quality data and adaptor sequences, we obtained over 670 million reads in total (Supplementary Table S3). Sequences ranging from 18 to 30 nt were extracted from all the reads and were mapped where possible to the corresponding cultivar genome; the percentage of reads that mapped perfectly to the genome was 96.3 (Supplementary  Table S3). Next, we removed rRNAs, tRNAs, snRNAs and snoRNAs from those reads that mapped perfectly to the genome, and finally obtained small RNA sequences comprising ta-siRNAs, ra-siRNAs, known miRNAs (osa-miRNAs), novel miRNAs (nov-miRNAs) and un-annotated RNA (Supplementary Table S3). From an analysis of size distributions, the frequencies of small RNAs were maximum at 21 and 24 nt in all samples (Fig. 1). Under normal growth conditions, sequences of 24 nt were a higher proportion than those of 21 nt, regardless of cultivar and developmental stage (Fig. 1A), while under cool conditions, sequences of 21 nt were a higher proportion than those of 24 nt at the UNM stage in both cultivars (Fig. 1B). Thus, cool-temperature stress appeared to promote the accumulation of 21-nt small RNAs at the UNM stage in rice. This may be a reflection of the stagespecific miRNA expression profile in anthers under cool stress, because the sizes of miRNAs and ta-siRNAs were mainly 21 nt, as determined by DCL1 and DCL4, respectively, in plants (Bartel, 2004;Voinnet, 2009) but ta-siR-NAs were rarely detected in any of our libraries compared with the level of miRNA (Supplementary Table S3).
To annotate osa-miRNAs, the reads of each library were mapped to the precursor and mature sequences from the rice miRNA database available in miRBase (release 19). As a result, a total of 3,765,844 reads were mapped   Table S3), and 239 out of 289 miRNA families available in the miRBase were detected in rice anthers in this study. Furthermore, 48 miRNA families were regarded as nov-miRNAs (Table 2), according to the criteria described in MATERIALS AND METHODS. After this analysis, we checked a more recent version of miRBase (release 21), and found that two of the 48 nov-miRNAs are now annotated as osa-miRNAs. Thus, we finally identified 241 osa-miRNAs and 46 nov-miRNAs in this study.

Characteristics of miRNAs expressed in rice anthers and their relationship with cool tolerance
The raw count of each osa-miRNA and nov-miRNA was normalized to RPM. We selected 110 miRNAs whose expression levels were over 10 RPM in at least one of the eight samples for further comparative analysis between Sasanishiki and Hitomebore. The 110 miRNAs were detected in rice anthers under both normal and cool conditions; neither normal-condition-specific nor coolcondition-specific miRNAs were included in this set of miRNAs. Under cool conditions, 89.1% of the miRNAs were common among the two stages (UNM and TCP) and the two cultivars ( Fig. 2A). At the UNM anther stage under cool conditions, most of the miRNAs were expressed at a higher level in Hitomebore than Sasanishiki (Fig.  2B), although this difference was not observed in the TCP stage (Fig. 2C). These data suggest that the differential accumulation of miRNAs between Sasanishiki and Hitomebore effects cool tolerance of rice anthers at the UNM stage, when the tapetum cells are being degraded (Oda et al., 2010). We therefore focused on the expression level of miRNAs at the UNM stage to explore miRNA(s) related to cool tolerance in rice in this study.
We selected miRNAs differentially expressed between Sasanishiki and Hitomebore at the UNM stage under the cool condition that satisfied both of the following conditions: 1) a value of more than 1 RPM was detected in both the cultivars, and 2) the RPM value differed by more than two-fold between the cultivars. We thus identified 13 miRNAs that were differentially expressed: seven osa-miRNAs (osa-miR397b, osa-miR398b, osa-miR408-3p, osa-miR528-5p, osa-miR171b/c/d/e/f, osa-miR1432-5p and osa-miR5072) and six nov-miRNAs (nov-miR017, nov-miR018, nov-miR023, nov-miR024, nov-miR030 and nov-miR031) ( Table 3). In addition, we identified two cultivar-specific miRNAs (osa-miR531a/c and nov-miR039) that accumulated to more than 10 RPM in one cultivar and less than 1 RPM in the other cultivar; osa-miR531 and nov-miR039 were specifically expressed in Sasanishiki and Hitomebore, respectively (Table 3). It is noteworthy that seven osa-miRNAs, out of the 13 differentially expressed and two cultivar-specific miRNAs, have previously been reported to be expressed in developing rice pollen (Wei et al., 2011). Thus, we considered that these 15 miRNAs were candidates for rice miRNAs related to cool tolerance in immature anthers. For the 15 candidate miRNAs, we calculated the RPM ratios between different growth conditions (cool/normal) to understand the up-or down-regulation of their expression under cool stress (Fig. 3). In Sasanishiki, four miRNAs (osa-miR398b, osa-miR1432-5p, osa-miR5072 and osa-miR531a/c) were up-regulated over two-fold by cool-temperature stress (Fig. 3A). In Hitomebore, seven miRNAs (osa-miR397b, osa-miR398b, osa-miR408-3p, osa-miR528-5p, osa-miR1432-5p, nov-miR018 and nov-miR031) were up-regulated over two-fold by cooltemperature stress (Fig. 3B). Most of the up-regulated miRNAs showed a greater change of RPM in Hitomebore than in Sasanishiki in response to cool stress (Fig. 3), indicating a greater influence of cool-temperature stress in Hitomebore than in Sasanishiki. Such differences in miRNA profiles between stress-sensitive and stress-tolerant genotypes have been observed in other plants subjected to abiotic or biotic stress (drought stress in cowpea, Barrera-Figueroa et al., 2011; salt stress in maize, Ding et al., 2009; drought or disease stress in soybean, Kulcheski et al., 2011). These examples suggest that plant miRNAs play an important role in stress tolerance.
To validate the miRNA expression profiles obtained using Illumina sequencing, qRT-PCR was conducted against five representative miRNAs (osa-miR397b, osa-miR398b, osa-miR408-3p, osa-miR528-5p and nov-miR018) (Fig. 4), which were selected as miRNAs showing a high cool/normal temperature RPM ratio in Hitomebore (Fig.  3B). The expression levels and patterns of miRNAs in the qRT-PCR results correlated almost exactly with those from Illumina sequencing, although a clear difference between cultivars was not detected in nov-miR018 expression (Fig. 4, I and J).
Prediction of target genes related to cool tolerance in rice at the booting stage Analysis of miRNAs enables computational prediction of miRNA targets because plant miRNAs have almost perfectly complementary sequences to their target mRNAs. Therefore, we predicted the target genes of the 15 candidate miRNAs using the psRNATarget web program with default parameters. In total, 64 target genes were predicted for 12 of the miRNAs, and these were categorized into 34 gene products (Table 4). No target genes were predicted for the other three miRNAs. Interestingly, the most frequent annotations were related to redox genes, with 20 target genes being redox-related: two genes (LOC_Os07g46990 and LOC_Os08g44770) targeted by osa-miR398b and osa-miR528-5p encode copper/zinc superoxide dismutase (SOD); 14 genes (e.g., LOC_Os01g44330) targeted by osa-miR397b and osa-miR528-5p encode laccase precursor; LOC_Os06g37150 targeted by osa-miR528-5p encodes Lascorbate oxidase precursor; two genes (LOC_Os01g03620 and LOC_Os01g03640) targeted by osa-miR528-5p contain a multicopper oxidase domain; and LOC_Os08g39960 targeted by osa-miR531 encodes 3-oxoacyl-reductase (Table 4). Of the 15 miRNAs, miR171c, miR398b and miR408 are highly conserved in Arabidopsis (Sunker and Zhu, 2004), and miRNA-target regulatory modules in Arabidopsis were also predicted in this study, such as miR171 and its target genes (scarecrow-like family genes) (Llave et al., 2002), miR398 and its target genes (copper/ zinc SOD genes, CSD1 and CSD2) (Sunkar et al., 2006),     and miR408 and its target gene (plastocyanin) (Abdel-Ghany and Pilon, 2008). These data suggest that our predictions of target genes show good reliability, and these miRNA-target gene properties may be conserved between plant species. We investigated expression levels and patterns of four representative predicted target genes (LOC_Os01g03640, LOC_Os02g20970, LOC_Os07g46990 and LOC_Os08g44770) by qRT-PCR, and obtained clear inverse correlations between expression levels of the target genes and those of their corresponding miRNAs (Fig. 5). Furthermore, we determined the cleavage sites of six predicted target genes using 5' RLM-RACE analysis (Supplementary Table S2; Fig. 6). Of the six target genes investigated, major cleavage sites of four target mRNAs (LOC_Os07g46990, LOC_Os07g38290, LOC_Os06g06050 and LOC_Os06g37150) were mapped between nucleotides 10 and 11 (counting from the 5' end of each miRNA) (Fig.  6). However, we obtained unexpected results for the remaining two target genes: 14 of 20 clones from LOC_Os02g44360 were cleaved between nucleotides 9 and 10 (Fig. 6E), and all 22 clones from LOC_Os01g71106 were cleaved between nucleotides 6 and 7 (Fig. 6F). Although we have no explanation for these slippages from canonical cleavage sites, Wei et al. (2011) observed a similar phenomenon in developing rice pollen, suggesting that this may be due to subtle differences in regulatory mechanisms between tissues. Thus, using qRT-PCR and 5' RLM-RACE techniques, we validated the computational prediction of target genes of miRNAs related to cool tolerance in rice at the booting stage.
Possible miRNA-mediated mechanisms of cool tolerance in rice anthers Tapetum cells in rice anthers are degraded by programmed cell death (PCD) during tetrad to vacuolated pollen stages (Li et al., 2006), and retarded degradation of the tapetum can cause pollen sterility (Li et al., 2006;Li et al., 2015). Pollen sterility in rice plants due to cool-temperature stress is mainly due to tapetum hypertrophy (Imin et al., 2006;Oda et al., 2010). Tapetum degradation, therefore, is essential for normal male gametophyte development in rice. Reactive oxygen species (ROS) accumulate in the tapetum during young microspore formation (Hu et al., 2011), and ROS accumulation is involved in tapetum degradation by PCD (Li et al., 2006). Thus, it is possible that SOD, one of the ROS scavengers (Mittler et al., 2004), plays a crucial role in the induction of male sterility under cool stress by disturbing PCD of tapetum cells in rice anthers. In this study, differential expression of osa-miR398b was observed between cool-sensitive and cool-tolerant rice cultivars, and the sequence of osa-miR398b indicates that it could target two SOD genes (LOC_Os07g46990 and LOC_Os08g44770) (Table 4). In Arabidopsis, SOD is regulated at the post-transcriptional level by miR398  (Sunkar et al., 2006;Beauclair et al., 2010); likewise, we clarified that one SOD (LOC_Os07g46990) was downregulated by osa-miR398b accumulation at the UNM stage under cool-temperature conditions in cool-tolerant Hitomebore (Fig. 5A). We also showed that this SOD (LOC_Os07g46990) transcript can be cleaved as a target of osa-miR398b (Fig. 6A). Similarly, accumulation of osa-miR528-5p may down-regulate the other SOD gene, LOC_Os08g44770 (Fig. 5C). These results suggest that inhibition of SOD by osa-miR398b and osa-miR528-5p contributes to the maintenance of normal tapetum degradation and subsequent pollen maturation in Hitomebore under cool stress. The involvement of miRNA-mediated down-regulation of SOD genes in cool tolerance in rice at the booting stage will be elucidated in the future by functional analyses of the candidate miRNAs and their targets.

Conclusion
Using NGS technology, we conducted a comprehensive comparative analysis of miRNAs in rice anthers under cool conditions from the cool-sensitive Sasanishiki and cool-tolerant Hitomebore cultivars. The miRNA profiles reported here are useful resources to deepen our understanding of the regulation of gene expression in plants under abiotic stress. The present study showed that a limited number of miRNAs are differentially expressed between the cultivars and are upregulated under cool stress at the UNM stage, indicating that small RNA-mediated silencing is involved in the coolstress response in rice anthers. These miRNAs and their targets are candidates for explaining the basis of the cooltolerance trait in Hitomebore anthers.