Small RNA Profiles of the Rice PTGMS Line Wuxiang S Reveal miRNAs Involved in Fertility Transition

MicroRNAs (miRNAs) play key roles in the regulation of plant growth and developmental processes. In this study, RNA-seq was used to examine the expression profiles of miRNAs in a novel, photo-thermosensitive genic male sterile (PTGMS) rice line, Wuxiang S (WXS), during fertility transition. A total of 497 known miRNAs and 273 novel miRNAs were identified. In a differential expression analysis, 26 miRNAs exhibited significant differential expression between WXS (Sterile, S) and WXS (Fertile, F). Some of these miRNAs were validated by quantitative real-time PCR. Among these miRNAs, 11 showed decreased expression levels, and 15 showed increased expression levels in WXS (S) compared to WXS (F). Some of these miRNAs, such as osa-miR156a-j, osa-miR164d, and osa-miR528, were shown to be negatively correlated with their targets. These targets have previously been reported to be related to pollen development and male sterility, suggesting that these miRNAs may be involved in the regulation of pollen development in the rice PTGMS line WXS. Furthermore, miRNA-mediated editing events were also observed. In this study, a possible model for the control of signaling pathways during the process of fertility transition in the rice PTGMS line WXS by miRNAs was developed. These findings contribute to our understanding of the roles of miRNAs during anther development in PTGMS lines in rice.


INTRODUCTION
MicroRNAs (miRNAs) are endogenous, small (20-24 nt), non-coding RNAs that play important roles in various biological and metabolic processes, including development, signal transduction, and biotic and abiotic stress responses (Bartel, 2004(Bartel, , 2009. In plants, primary miRNA transcripts (pri-miRNA) are mainly transcribed by RNA polymerase II and subsequently cleaved by Dicer-like1 (DCL1), resulting in mature miRNA sequences (Voinnet, 2009). Mature miRNAs then enter into the RNA-induced silencing complex (RISC) and negatively regulate gene expression at the post-transcription or translation level by degrading target mRNAs or by translational repression (Bartel, 2004).
A deep-sequencing study of the transcriptome detected more than 100 miRNAs in plant male gametophytes during development (Le et al., 2011). Recently, differential expression patterns of miRNAs between the cytoplasmic male sterility (CMS) line and its maintainer line have been reported in maize (Shen et al., 2011), cotton (Wei et al., 2013), Brassica juncea (Yang et al., 2013), cabbage (Wei et al., 2015), and rice (Yan et al., 2015). Additionally, it has been found that miRNA172 can control flowering time by down-regulating APETALA2-like target genes in Arabidopsis (Chen, 2004). Moreover, in Arabidopsis, miR167 overexpression has been reported to lead to male fertility defects (Sire et al., 2009), whereas miR159a overexpression results in decreased expression of MYB33, leading to male sterility and flowering time delays in Arabidopsis (Achard et al., 2004). Other studies have also reported that MYB33/MYB65 facilitates Arabidopsis anther development. These mutants were conditionally sterile but exhibited increased fertility under greater light levels or lower temperatures (Millar and Gubler, 2005). In photoperiod-and thermo-sensitive genic male sterile (PTGMS) rice lines such as PA64S and NK58S, a 21-nt small RNA (sRNA, osa-smR5864w) is produced by a non-coding RNA precursor p/tms12-1. A C-to-G substitution in the small RNA leads to loss-of-function, which gives rise to the male sterility found in PA64S and NK58S (Ding et al., 2012;Zhou et al., 2012Zhou et al., , 2014. However, our knowledge of the involvement of miRNAs in rice PTGMS lines is still incomplete. Hybrid rice breeding has made a tremendous contribution to food security in China (Cheng et al., 2007;Peng et al., 2008), and the utilization of male sterility in hybrid rice technologies is mainly based on three-line and two-line systems. The threeline system includes the CMS line, a maintainer line that maintains the sterility of the CMS line, and a restorer line that restores the fertility of the CMS line in hybrid rice. However, in the two-line system, the fertility of the male sterile line is influenced by the alteration of environmental conditions (daylength and temperature); thus, the male sterile line is called an environmentally sensitive genic male sterile (EGMS) or a PTGMS line. These male sterile lines can be used not only as male-sterile lines but also as maintainer lines depending on environmental factors. Therefore, the application of a two-line breeding system is simple, inexpensive, labor-saving, and effective and overcomes the limitations of the three-line system in hybrid rice (Liu et al., 2001;Yang et al., 2007;Xu et al., 2011;Zhou et al., 2012). Recently, we developed a novel male sterile line Wuxiang S (WXS) for use in two-line systems. This line was derived from a mutant tms5 locus in indica rice. In this study, we used RNA-seq to identify conserved and novel mircoRNAs that may be involved in fertility transition during pollen development in WXS. The interaction network between these miRNAs and their targets was also investigated using GO and KEGG analysis, and we attempted to elucidate the potential regulatory mechanism of pollen development during fertility transition in the rice PTGMS line WXS.

Rice Materials
The rice PTGMS line WXS was generated by our laboratory. Between May and August 2014, WXS plants were grown in the natural ecological paddy field of Huazhong Agricultural University (30 • 28 ′ N, 114 • 21 ′ E), Wuhan, Hubei province of China. From 20 July to 3 August 2014, natural ecological conditions were daily average temperatures between 25.5 • C and 34.6 • C, with approximate day lengths of 14 h light and 10 h dark in Wuhan (Wuhan Weather Bureau), which could induce WXS sterility, and no pollen was produced in male sterile rice, a condition designated as WXS (S). When the young panicle length was approximately 1 cm, we transferred 20 of these WXS plants from the paddy field into a cooling pond to treat them with lower temperatures (approximately 21 • C) and shorter day lengths (approximately 12 h light/12 h dark) for 2 weeks. These plants were converted into male-fertile plants, designated WXS (F). Young panicles were separately collected from the natural ecological paddy field and from the cooling pond at the pollen mother cell (PMC) formation stage (P2) and the meiosis stage of PMC (P3) and were correspondingly named SP2, SP3, FP2, and FP3. Collected samples were frozen in liquid nitrogen and then stored at −80 • C for future use. To analyze pollen fertility, mature anthers and pollen grains were also collected and stained with 1% potassium iodide solution (I 2 -KI). PMCs were stained using carbol fuchsin dye for cytological observations. Rice anther staging was performed in this study according to the method of Itoh et al. (2005).

RNA Isolation, Small RNA Library Construction, Sequencing, and Data Analysis
Total RNA was extracted from the young rice panicles at the PMC and meiosis stages using Trizol reagent (TaKaRa, Dalian, China) according to the manufacturer's instructions. RNA quality and quantity were measured using a NanoDrop2000 and Agilent 2100 bioanalyzer. Then, small RNA libraries were constructed and analyzed as described by Molly and Thomas (2010). Briefly, small RNAs ranging from 18 to 30 nt were size fractionated using 15% polyacrylamide gel electrophoresis (PAGE), and suitable isolated fragments were extracted from the gel and ligated with 5 ′ and 3 ′ RNA adaptors using T 4 ligase. Subsequently, the resulting products were reverse transcribed and amplified using 15 cycles of PCR to produce sequencing libraries. Finally, flow cell sequencing was performed on the HiSeq 2500 platform. Automated base calling of the raw sequences and vector removal were performed with the PHRED and CROSS MATCH programs.

Differential Expression Analysis of miRNAs
To identify miRNAs with differential expression in the four samples (SP2 and FP2, SP3 and FP3), miRNA read counts were normalized to transcripts per million (TPM) using the following formula: normalized expression = (miRNA count/total count of clean reads) ×10 6 . If one miRNA was not expressed in one of two samples, the normalized read count of this miRNA was arbitrarily set as 0.01 for further calculation (Chen et al., 2012). When miRNA expression levels were less than 1 TPM in both samples, differential expression analysis was not performed. Fold changes and P values were calculated to determine the significance of expression differences between SP2 and FP2 and between SP3 and FP3. Hierarchical clustering and k-means clustering of the expression patterns were performed in Mutiexperimental Viewer v4.7 (Saeed et al., 2006).

Prediction of the miRNA Targets and Functional Analysis
The sequences of all known and novel miRNAs were aligned to the annotated rice genome (MSU-v7.0) to predict potential target genes. The criteria used for target prediction in this study were as described by Hu et al. (2014). The potential targets of known and novel miRNAs that were significantly differentially expressed between two samples (SP2-FP2 or SP3-FP3) were searched against the Gene Ontology database and KEGG databases (Kanehisa and Goto, 2000). The biological process, molecular function and cellular component of the targets were obtained, and then GO functional enrichment analysis was performed using AgriGO with Oryza sativa MSU-v7.0, and graphical results depicting overrepresented GO terms were generated via singular enrichment analysis to fully understand the function of these targets.

miRNA Detection by Stem-Loop RT-PCR and Validation by qPCR
After total RNA was extracted from the young panicles at the PMC and meiosis stages, RNA-free DNaseI (Fermentas, USA) was used to remove DNA contamination for 30 min at 37 • C. Approximately 2 µg of total RNA was reverse transcribed using miRNA-specific stem-loop primers and a Fermentas RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA), as described previously (Varkonyi et al., 2007). The reactions were incubated for 30 min at 16 • C; followed by 60 cycles of pulsed reverse transcription at 30 • C for 30 s, 42 • C for 30 s and 50 • C for 1 s; and finally terminated by incubating at 70 • C for 5 min. For miRNA targets, cDNA templates were reverse transcribed using the Oligo dT 18 primer with a Fermentas RevertAid First Strand cDNA Synthesis Kit (Fermentas, USA).
Expression analysis of the known miRNAs and their potential targets was performed using the ABI Step One Plus TM Real Time PCR System (Applied Biosystems, USA) and a SYBR Green Master Mix (Roche, Germany). Quantitative real-time PCR (qPCR) was performed using the following parameters: 10 min at 95 • C, followed by 40 cycles of 15 s at 95 • C and 60 s at 60 • C. U6 snRNA and OsActin were chosen as endogenous controls for miRNAs and targets, respectively. The reactions were performed with three biological replicates, and a melting curve analysis was carried out to verify that only one specific amplification occurred. Comparative expression levels were calculated in the four different samples using the 2 −△△CT method (Livak and Schmittgen, 2001). Novel miRNAs were detected by electrophoresis analysis after stem-loop RT-PCR. All of the primers used in this study are listed in Table S1.

Identification of Potential miRNA Editing Sites and Validation of These Sites
The RNA editing level was calculated by the ratio of reads supporting the mismatch in a site to the total reads detected at that site. MiRNA editing sites located on known miRNA sequences were identified according to the methods described by Luciano et al. (2004) and Yang et al. (2013). Genomic DNA and total RNA were isolated from young panicles to validate miRNA editing. Precursor miRNA sequences were amplified from rice DNA, and the corresponding mature miRNA sequences were also amplified from cDNA, which was reverse-transcribed using stem-loop RT-PCR. All primers used in this study are listed in Table S1.

Cytological Observation of WXS under Fertility Transition
WXS was grown in the paddy field until the panicle length was approximately 1 cm, and then 20 plants were selected and transferred to a cooling pond to treat them with low temperatures for 2 weeks. Under the low temperature treatment (∼ 21 • C), WXS was male-fertile with normal anthers ( Figure 1A). Pollen grains from these samples were stained using a 1% I 2 -KI solution ( Figure 1B). Fertile pollen grains maintained the normal morphology from pollen mother cell formation (Figures 1E,a) to mature pollen grains ( Figures 1E,b-d). In contrast, under natural conditions, the WXS plants had thin anthers ( Figure 1C) and a complete absence of I 2 -KI-stained pollen ( Figure 1D), exhibited sterility with abnormal pollen mother cells (Figures 1E,e,f) and aberrant dyads and tetrads during meiosis (Figures 1E,g-k), and eventually produced cracked pollens (Figures 1E,l). These results suggested that the male sterile line WXS can display two different phenotypes under two different environmental conditions; these were named WXS (S, Sterile) and WXS (F, Fertile).

Overview of Small RNA Library Sequencing Data
To identify the roles of miRNAs in the process of fertility transition during anther development in WXS, four small RNA libraries (SP2, SP3, FP2, and FP3) were constructed for deep sequencing. A total of 9,114,807, 6,775,838, 22,076,091, and 9,560,446 raw reads were generated via Illumina sequencing from the SP2, SP3, FP2, and FP3 libraries, respectively (Table S2). After . Carbol fuchsin dyeing of the WXS (S) (E,a-d) and WXS (F) (E,e-l) microsporocytes from the pollen mother cell formation stage to the pollen meiosis stage. a, e, f were considered the pollen mother cell formation stage; b, c, g, h, i, j, k were considered the meiosis stage; d, l were considered the mononuclear stage.
When the common and unique reads of these small RNAs were compared between the four libraries, more than 64% of the total sRNAs were shared when any two libraries were compared: 64.58% were shared between SP2 and SP3, 67.14% were shared between SP3 and FP3, 68.48% were shared between SP2 and FP2, and 69.14% were shared between FP2 and FP3 (Figure 2A). Correspondingly, more than 11% of unique small RNAs were found when any two libraries were compared, suggesting that there was a less abundant but much more diverse pool of small RNAs that could be assumed to represent induced-specific small RNAs (Figure 2A). These data emphasize the differences and complexities in the assemblages of small RNAs between the different environments in the process of fertility transition in WXS.

Identification of Known miRNAs
In plants, the first nucleotide of a miRNA is important for the construction of the RNA induced silencing complex (RISC). In Arabidopsis, different Argonaute (AGO) proteins are preferentially enriched for the different initiating nucleotides of miRNAs; for example, AGO1 prefers miRNAs beginning with a U, whereas AGO2 and AGO4 recruit miRNAs with a 5 ′ terminal A (Mi et al., 2008). In our study, we found that 48% of the first nucleotides of miRNAs were U, followed by C at 26%, A at 18%, and G at 8% (Figure 2B). These results were similar to those of a previous study in which most of the miRNAs in rice began with a U (Yan et al., 2015). To identify known miRNAs, unique small RNA sequences were aligned against mature miRNAs from rice deposited in the miRBase database (Release 20) (Kozomara and Griffiths, 2011). Currently, approximately 713 mature miRNAs derived from 592 pre-miRNAs are published in the miRBase (June 2013). In our study, the length distribution of small RNAs was similar between the four libraries, with fragment sizes mainly between 21 and 24 nt in length ( Figure 2C). A total of 497 known miRNAs were identified, of which 411 miRNAs were in SP2, 446 in FP2, 383 in SP3, and 410 in FP3 (Table S3). Among these miRNAs, 339 known miRNAs were shared between all four RNA libraries, accounting for 68.21% of the identified miRNAs ( Figure 2D). Moreover, a total of 64 miRNAs, of which 13 were in SP2, 31 were in FP2, 6 were in SP3, and 14 were in FP3, were detected in only one of the four samples ( Figure 2D). This indicates that those miRNAs are likely to be specifically expressed in the different anther development stages in WXS.

Identification of Novel miRNAs
A large number of small RNAs generated by high-throughput sequencing were used to identify novel miRNAs. In this study, a total of 373 novel miRNAs were found, many of which had low expression levels (Table S4). Of these miRNAs, 201 miRNAs were located on the 5 ′ arm of the miRNA precursor, and 160 miRNAs were located on the 3 ′ arm; 12 miRNA * sequences were also obtained ( Table S4). The miRNA* sequence plays an important role in miRNA prediction. However, because most miRNA * s are degraded soon after their separation from their corresponding miRNAs, it is thought that miRNA * s are usually rather low in abundance or are even undetectable (Meyers et al., 2008).
Fragment size distribution analysis indicated that most of the novel miRNAs displayed a nucleotide length of 21 nt ( Figure S1). In plants, the determination of which strand of the miRNA:miRNA* duplex is incorporated into the RISC is largely based on the identity of the first nucleotide. In the present study, the first nucleotide analysis of these novel miRNAs showed that the base rates of U, C, A, and G were 49, 26, 17, and 9%, respectively ( Figure S1). These novel miRNAs exhibited differential expression profiles. To validate the predicted novel miRNAs, five novel miRNAs, including novel-miR-32, novel-miR-169, novel-miR-221, novel-miR-259, and novel-miR-368  ( Table S4), were selected for confirmation by stem-loop RT-PCR. PCR indicated that the DNA fragments were approximately 60-bp long. Moreover, the secondary structures of these novel miRNAs were predicted using MFOLD (http://mfold.rna.albany. edu/) and manually checked. The results showed that the five novel miRNAs had perfect secondary structures ( Figure 3A). Both the PCR results and the predicted secondary structures showed that these novel miRNAs are real and reliable and revealed that these five novel miRNAs were expressed in the young panicles of WXS ( Figure 3B).

Comparative miRNA Expression Profiles of WXS during Fertility Transition
Differentially expressed miRNAs were identified in this study based on the normalization of the read counts. The relative abundances of 72 known miRNAs were found to be significantly different ( * P < 0.05; * * P < 0.01) during fertility transition (SP2-FP2 (WXS (S)-WXS (F)) or SP3-FP3(WXS (S)-WXS (F)) ( Table S5). Fifty-five and 51 miRNAs were significantly differently expressed in SP2-FP2 and SP3-FP3, respectively ( Figure S2). Among them, expression levels of 26 miRNAs were significantly different between WXS (S) and WXS (F) ( Table 2). In this study, 10 different clusters were obtained from Hierarchical Clustering (HCL) analysis of the expression patterns using the software MeV ( Figure 4A). The results show that many miRNAs were only up-regulated in SP2 or SP3 of WXS (S) (Clusters 3, 7, 8, and 9), though some miRNAs were up-regulated in both SP2 and SP3 of WXS (S) (Clusters 2 and 10). In contrast, some miRNAs appeared to be downregulated only in FP3 of WXS (F) (Cluster 4). The differential expression of these miRNAs (Cluster 4) suggests that they may be related to male sterility during fertility transition. Their high expression levels during the meiosis stage (SP3 and FP3) might inhibit some genes with essential roles in STD represent normalised expression level of miRNAs. Normalised expression (TPM) = count of miRNA/total count of clean sRNAs 9 × 10 6 . " ** " and " * " represnent significant differential expression and significant differential expression, respectively.  anther development. Eight known miRNAs (miRNA156a-j, miRNA3979, miRNA159c/d/e, miRNA171a/c/e/i, miRNA398b, miRNA164d, miRNA528, and miRNA408) were selected to authenticate expression profiles. The results showed that in WXS (S), expression levels of miRNA156a-j and miRNA171a/c/e/i were down-regulated ( Figure 4B), whereas the rest of these miRNAs were up-regulated. These results were consistent with the high-throughput sequencing data, indicating that the results of this study were reliable.
Our results also showed that the MYB transcription factor (LOC_Os01g59660) and No apical meristem protein (LOC_Os02g36880) were both predicted to be targets of osa-miR164d. However, the expression patterns of the two targets were different during the process of fertility transition. Transcripts of LOC_Os01g03720 were positively correlated with the expression of miR164d (Figure 5C), whereas the expression level of LOC_Os02g36880 was negatively correlated (at nearsignificant levels) with miR164d ( Figure 5D). The results suggest that the expression levels of different targets of the same miRNA might be regulated differently. Some negative correlations were also found between expression levels of target genes and their corresponding miRNAs in WXS (S) and WXS (F), indicating that miRNA-mediated mRNA silencing occurred during anther development in the fertility transition.

miRNA Target Prediction and Functional Analysis
In plants, miRNAs regulate gene expression by interacting with their targets. Identifying the candidate genes targeted by the miRNAs would contribute to our understanding of the biological functions of miRNAs. Using psRNATarget, we identified 5192 targets for 435 known miRNAs, for an average of 11.9 targets per miRNA in this study (Table S6). We also identified multiple targets for the miRNAs at the 3 ′ terminal, such as Nramp6 for miR156c-3p and OsFBDUF58 for miR160a-3p, etc. As previously described, most of the predicted miRNA targets are particular transcription factors. We also identified some targets other than these conserved targets for conserved miRNAs. For example, besides NAC factors, osa-miR164 also targets the No apical meristem protein (Table S6). Predicted targets consist mainly of transposon proteins, retrotransposon proteins, growth-regulating factors, MYB family transcription factors, F-box domain-containing proteins, MADS-box family proteins, and SBP-box gene family members ( Table S6). Targets of novel miRNAs were also predicted and show a much broader range of potential functions. In addition to transcription factors, various enzymes, transposon proteins, and unidentified expressed proteins were also targeted by novel miRNAs (Table S7). Based on the GO annotations (Figure 6) from this study, targets were enriched for the "metabolic process, " "cellular process, " and "single-organism process" terms of the "biological process" cluster. "Cell", "cell part, " and "organelle" were the three most abundant "cellular component" terms. In the "molecular function" cluster, the top three terms were "binding, " "catalytic activity, " and "nucleic acid binding transcription factor activity" (Figure 6).

Observation of the miRNA Editing Events
RNA editing is another post-transcriptional modification that generates divergence between RNAs and their genomic DNA sequences. Previous studies have documented that miRNA editing occurs in rice during the grain-filling stages (Yi et al., 2013). In this study, a large number of miRNA editing events were found during fertility transition (Table S9). These miRNA editing events mainly occurred at nucleotide positions 11 and 18 from the 5 ′ terminus (Figure 7A). The most abundant nucleotide substitution was U to A, which accounted for over 30% of substitutions in the SP2 and FP3 libraries; this was followed by the U to C type substitution, which comprised more than 25% (Figure 7B). These results were similar to those reported by Yan et al. (2015).
To validate the occurrence of miRNA editing, three editing types (U to A, U to C, and A to G) were examined in osa-miR166m, osa-miR2118h, and osa-miR399a, respectively. Using PCR amplification, precursor miRNA sequences from genomic DNA and mature miRNA sequences from cDNA were cloned to characterize nucleotide substitutions. This sequence comparison analysis further confirmed that miRNA editing events were reliably detectable during miRNA generation ( Figure 7C).

DISCUSSION
Previous reports have studied the relationship between miRNAs and CMS in maize (Shen et al., 2011), cotton (Wei et al., 2013), Brassica juncea (Yang et al., 2013), and rice (Yan et al., 2015). However, few studies have investigated miRNAs from rice PTGMS lines. Therefore, characterizing the role of miRNAs in rice PTGMS lines would be extremely useful and could contribute to an improved understanding of the molecular functions of miRNAs during fertility transition in the male sterile rice lines used in two-line systems. In this study, we examined miRNA profiles determined by high-throughput sequencing and qPCR to investigate the expression of miRNAs in the rice PTGMS line WXS. These analyses revealed that miRNAs targeted many transcription factors associated with metabolism and signal transduction and play important roles in anther development. In this study, we found a total of 26 conserved miRNAs that were differentially expressed between WXS (S) and WXS (F) during fertility transition ( Table 2). Of these, 11 miRNAs were downregulated and 15 miRNAs were up-regulated. Of these, several members of the miR156 family showed significant differential expression levels in the four libraries ( Table 2). It was apparent that expression levels of miRNAs were higher in WXS (F) than in WXS (S) during fertility transition.
In rice, it is well known that miR156 targets SBP-box gene family (SPL) proteins, which play important roles in the proper development of sporogenic tissues. Previous studies have reported that miR156 regulates the timing of flower formation through its target SPL3, which activates the expression of APETALA1 (Yamaguchi et al., 2009). In the present study, the low expression levels of miRNA156 found in WXS (S) might affect flower formation during pollen development. Many studies have reported that the formation of pollen exine is important for pollen development (Qin et al., 2013). In loblolly pine, laccase genes were reported to be expressed and involved in plant pollen development (Turlapati et al., 2011). The osa-miR3979 target laccase was found to participate in regulating the pollen abortion process in a maize CMS line (Shen et al., 2011). However, in the present study, osa-miR528 was predicted to target the laccase precursor protein (LOC_Os01g62600). Furthermore, our qPCR analysis confirmed that the expression of osa-miR528 is negatively correlated with that of its target laccase (LOC_Os01g62600) (Figure 5H). These results indicated that osa-miR528, which targets laccase, might be involved in the fertility transition of WXS. In all eukaryotes, the ubiquitous, multifunctional calcium sensor calmodulin (CaM) mediates calcium action by regulating the activity and function of many proteins. Calcium and CaM play crucial roles in pollen germination and pollen tube growth. In Arabidopsis, a calmodulin-binding protein was reported to be essential for pollen germination (Golovkin and Reddy, 2003). Previous studies have also documented that osa-miR1432 and osa-miR812d may be involved in Ca 2+ -mediated signaling pathways by targeting genes that encode EF-hand family proteins and CaM/Cadependent protein kinase, respectively (Yan et al., 2015). In our study, osa-miR5976 was predicted to target a gene coding for a calmodulin binding protein (LOC_Os12g36940), and expression of osa-miR5976 was up-regulated in WXS (S) ( Table 2).
The conserved targets of miR159 are MYB transcription factors, which have been reported to be involved in flower development and are essential for fertility (Jones et al., 2006;Tsuji et al., 2006). In rice, mutations in OsGAMYB have resulted in defects in anther and pollen development, and overexpression of osa-miR159 leads to male sterility (Kaneko et al., 2004;Tsuji et al., 2006). In this study, three members of osa-miR159 (osa-miR159c/d/e) were found to be up-regulated in WXS (S) (SP2 and SP3). The expression levels of osa-miR159 and its target, GAMYB transcription factor (LOC_Os01g59660), were confirmed by qPCR, which showed that these miRNAs and their targets were negatively correlated ( Figure 5B). Therefore, it is highly likely that osa-miR159 silenced the expression of MYB proteins, which affect anther development in WXS (S). Moreover, relatively high expression levels of osa-miR164d were found in the meiosis stage of WXS (S) (SP3). The gene coding for the transcription factor NAC is considered to be the target of osa-miR164. In many plant species, NAC transcription factors are involved in plant tolerance to biotic or abiotic stresses. The SP2 to SP3 stages in WXS (S) are important for fertility transition during anther development. Up-regulated expression of osa-miR164d in WXS (S) (SP3) indicated that it may be involved in fertility transition. Moreover, osa-miR444f was predicted to target the gene encoding a pentatricopeptide repeat (PPR) protein (LOC_Os04g14450). PPR proteins are a large family with tandem repeats of a degenerate 35 amino acid motif. In higher plants, PPR proteins can suppress the expression of mitochondrial genes associated with cytoplasmic male sterility to restore fertility (Desloire et al., 2003;Wang et al., 2004). In our study, the higher expression of osa-miR444f in WXS (S) may silence the expression of PPR proteins, ultimately affecting their functions in anther development. Furthermore, previous reports have also shown that accumulation of two Arabidopsis miRNAs (miR171, miR398) oscillates during the diurnal cycle: increasing during the light period and decreasing in darkness (Christelle et al., 2009). The miRNA171c targets SCL6-II, SCL6-III, and SCL6-IV regulate shoot branching in Arabidopsis (Wang et al., 2010;Manavella et al., 2013). In this study, miRNA398b was down-regulated in WXS (S) ( Table 2). However, miRNA171, which showed higher expression in WXS (F), may target an ankyrin repeat domain-containing protein (ARDCP, LOC_Os03g04300).
Based on the comprehensive analysis described above, we proposed a possible model by which miRNAs control the signaling pathways involved in fertility transition of the rice PTGMS line WXS (Figure 8). MiR171, miRNA156, and miRNA3979 interact with their respective target genes (ARDCP, SPL, and GRMP) to modulate phase transition. These interactions then result in the subsequent formation of organs and the meristem and the accumulation of biomass, and finally result in morphological changes to the rice anther. Moreover, the interactions of miRNA444, miRNA159 and miRNA164 with their corresponding target genes (PPR, GAMYB, and NAC) modulate the expression of developmental genes, GA/ABArelated genes, and auxin-responsive genes to promote the transmission of signals that enhance developmental processes and maintain energy supply. These properties affect development and metabolic processes. Fertility transition in the rice PTGMS line WXS may occur due to combined action of the regulatory interactions described above. Although fertility transition in rice PTGMS lines is a highly complex process, it is certain that miRNAs must play very important roles in the regulation of rice anther development. This study also showed that miRNAs may play a critical role in the fertility transition stage, which adapts to variations in temperature and photoperiod. Our functional description and target analysis of miRNAs provide more clues for understanding the different ways gene expression is regulated during fertility transitions. These miRNA profiles also revealed an important component of the gene regulatory circuit and may provide insights for further investigations of PTGMS rice lines.

CONCLUSION
In this study, we performed high-throughput sequencing to identify the miRNA expression profiles of the rice PTGMS line WXS during fertility transition. A total of 497 known and 273 novel miRNAs were identified, and some of these miRNAs were validated by stem-loop RT-PCR analysis. MiRNA editing events were also observed and validated during anther development. Target prediction analysis indicated that some miRNAs are involved in anther development and male sterility in the WXS line. The characterization and comparative expression profiling of the miRNA transcriptome in this study lays the foundation for understanding the complex miRNA-mediated regulatory networks in rice anther development. A possible model for the control that miRNAs exert on signaling pathways during fertility transition in the rice PTGMS line WXS was proposed. Further functional studies on differentially expressed miRNAs will provide a better understanding of miRNA-mediated regulation mechanisms during fertility transitions in rice PTMGS lines.

AUTHOR CONTRIBUTIONS
YD and HZ conceived and designed the experiments. HZ, QQ, HC, JJ performed the experiments. HZ and JH analyzed the data and drafted the manuscript. YD contributed reagents/materials/analysis tools and modified the manuscript. All of the authors carefully checked and approved this version of the manuscript.

DATA ACCESS
The miRNA raw data released have been deposited at NCBI in the Gene Expression Omnibus (GEO) database under the accession number: GSE74003 (http://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc=GSE74003).

FUNDING
This study is supported by National Natural Science Foundation of China (31471464), and "973" Program of China (2013CB126900).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 00514 Table S1 | Primers used in this study.     Figure S1 | Fragment size distributions and first nucleotides of the novel miRNAs. The fragment size distribution analysis indicated that most of the miRNAs have nucleotide lengths of 21 nt (A). The first nucleotide analysis of these novel miRNAs showed that the bases rates of U, C, A, and G appeared at frequencies of 49, 26, 17, and 9%, respectively (B).