Profiling miRNA expression in photo-thermo-sensitive male genic sterility line (PTGMS) PA64S under high and low temperature

ABSTRACT

Peiai64S (PA64S) is a photo-thermo-sensitive genic male sterile line (PTGMS), with wide application in hybrid seed production in rice (Oryza sativa L.). Micro-RNAs are 21–24 nt, endogenously expressed small RNAs that have been characterized in various developmental stages of rice, but none have been studied with respect to the regulation of TGMS in rice. Here, we employed high-throughput sequencing to identify expression profiles of miRNAs in the anthers of PA64S at high (PA64S-H) and low temperature (PA64S-L). Two small RNA libraries from PA64S-H and PA64-L anthers were sequenced, and 263 known and 321 novel candidate miRNAs were identified. Based on the number of sequencing reads, a total of 133 known miRNAs were found to be differentially expressed between PA64S-H and PA64S-L. Target prediction showed that the target genes encode MYB and TCP transcription factors, and bHLH proteins. These target genes are related to pollen development and male sterility, suggesting that miRNA/targets may play roles in regulating TGMS in rice. Further, starch and sucrose metabolism pathways, sphingolipid metabolism, arginine and proline metabolism, and plant hormone signal transduction pathways were enriched by KEGG pathway annotation. These findings contribute to our understanding of the role of miRNAs during anther development and TGMS occurrence in rice.

KEYWORDS:

• Mirna
• TGMS
• Oryza sativa L

Introduction

Rice (O. sativa L.) is one of the most important global grain crops and provides a staple food for almost half of the world population.¹ The development of hybrid rice is the main way to increase its yield potential, because hybrid rice varieties have a yield advantage of about 20% or higher than the improved inbred varieties.² Hybrid rice technologies are mainly based on two male sterility systems, namely environmentally sensitive genic male sterile (EGMS) and cytoplasmic male sterile systems (CMS).³ Different from the three-line system (CMS), the two-line system is based on the discovery and application of thermosensitive genic male sterile (TGMS) lines, which serve as both the male sterile lines and maintainer lines under different environmental conditions.⁴ The advantages of the two-line system compared to the three-line system include the wider germplasm resources used for breeding parents, better grain quality, higher yields, and simpler breeding and hybrid seed production processes.⁵

PA64S is one of the most important indica rice genic male sterile (GMS) lines in the two-line system, which was first identified as the photoperiod-sensitive GMS line (PGMS) by Mingsong Shi in Hubei, China in 1973,⁶ where the development of its anthers are sensitive to length of photoperiod, with abortive pollen grains generated under long photoperiod condition, and with fertile pollen grains occurred under short photoperiod condition,during the meiosis period of anther development. PA64S was subsequently found to be thermal-sensitive in pollen development, and is also a TGMS line, where at high temperature, it is male-sterile and produces sterile pollen, but at low temperature it converts to male-fertile plants that produce fertile pollen.⁷ Compared to other TGMS lines, PA64S is a indica-type (O. sativa indica) TGMS line with an japonica-type NK58S-derived sterility gene,⁸ and has wide compatibility and good agronomic traits. With its lower critical point temperature of fertility alteration (~23°C), PA64S is a valuable male sterile germplasm source and has become the most widely used female parent for two-line hybrid rice breeding.⁸

Extreme temperatures are a key factor limiting global rice distribution. Super hybrid rice cultivars grown in tropical and subtropical climates are highly productive, but they are often harmed by low and/or high temperatures.⁹ Heat and cold stress seriously affect the normal growth of rice throughout its life span.¹⁰ During the rice life cycle, the reproductive phase is the most sensitive to temperature stress, and even small changes in temperature can cause significant losses to food crops during flowering.¹¹ Cold and heat tolerance are complex traits regulated by multiple genes, and these traits are highly genotype dependent. As temperature stress and climatic changes become more prevalent, it is urgent to identify genes associated with temperature stress tolerance and to understand their regulatory mechanisms to develop crops with enhanced temperature stress tolerance through genetic manipulation.¹²

Recent studies have documented that MicroRNAs (miRNAs) are fundamental sequence-specific regulatory elements (~24 nucleotides) and play key roles in the regulation of growth and development of eukaryotic organisms.^13,14 In plants, mature miRNAs are processed from stem-loop regions of long primary transcripts by a Dicer-like 1 (DCL1) enzyme, then incorporated into the RNA-induced silencing complexes (RISC), guiding RISC to target mRNA, where they silence the expression of specific miRNAs by directing either mRNAs cleavage or translational repression.^14–18 The miRNAs are non-coding, single-stranded RNAs encoded by endogenous genes involved in post-transcriptional regulation of target gene expression in plants and animals, and most miRNA target genes are transcription factors.^19–21 Both miRNA and transcription factors can regulate the expression of downstream genes. Most of the miRNAs down-regulate the expression of target genes, and transcription factors activate or inhibit genes to determine the fate of target genes.^22,23 In short, miRNAs not only regulate transcription factor expression, but also interact with transcription factors, all of them forming a complex regulatory network that participates in plant response mechanisms under developmental and stress conditions.²⁴ Accordingly, a number of miRNA functions have also been identified in various biological processes, including developmental regulation, metabolism, and response to biotic and abiotic stresses.^25–27

As a class of negative regulators, miRNA biogenesis and activity is well characterized in sporophytic tissues,²⁸ and play a broad regulatory role during male gametophyte development in plants.²⁹ Studies have reported that miR390, miR166, miR159, miR530, miR167, miR164, miR168 and two new miRNAs are differentially expressed in photoperiod-sensitive male-sterile mutant 7B-1 and its wild type. The predicted targets of these miRNAs are potentially involved in anther development and regulation of male-sterility in tomatoes.³⁰ Recent transcriptomics studies identified 24 conserved miRNAs that are differentially expressed in the pollen of a cytoplasmic male-sterile line of rice and its maintainer line. Targets of these miRNAs include the MYB family protein, PPR-containing proteins, and kinases, which are involved in carbohydrate metabolic pathways, signal transduction, and reproductive processes.³¹

There is increasing evidence showing that miRNAs play key roles in the regulation of anther development. In plants, miR159 generally regulates the expression of GAMYB or GAMYB-like genes, which are involved in the regulation of anther development and microsporogensis.^32–34 In rice, the gamyb deletion mutant is aborted due to the inability of the tapetum to perform PCD (programmed cell death) and degradation.^35,36 Overexpression of miR159 in Arabidopsis thaliana and cereal crops results in sterile anthers.^37–39 Thus, the role of GAMYB in anthers is highly conserved and regulated by miR159. In addition, miR159 interacts with miR319 and miR167. The expression of the target genes GAMYB and TCP4 are inhibited by miR159 and miR319, respectively. They can regulate the expression of miR167a alone or in the form of a protein complex to regulate the activity of ARF6/8. Therefore, the miR159-miR167-miR319 network is very important for the development of reproductive organs, including the ripening of anthers.^40,41

To our knowledge, there has been no report comparing miRNA expression profiles of TGMS during reproductive development under high and low temperatures in rice. To do this, we used high-throughput sequencing technology to investigate the differential miRNA expression and targets from rice anthers at the early uninucleate stage of PA64S at high (PA64S-H) and low temperature (PA64S-L). Differential expression patterns of miRNAs were analyzed between PA64S-H and PA64S-L. Targets were predicted, and their expression profiles were selectively validated. These results may shed light on the regulatory roles of miRNAs related to temperature during pollen development as well as the occurrence of TGMS in rice.

Materials and methods

The TGMS line, PA64S, was used in this study. Plants were grown in the field until the fertility sensitive period at Hunan Normal University (latitude 28°19ʹ’ N; longitude 112°95ʹ’ E) under standard growing conditions. At jointing-booting stage, rice (PA64S-H) were transferred to a growth chamber at 30°C for 2 weeks of high temperature treatment, whereas rice (PA64S-L) were placed in a cooling pool at 22°C for 2 weeks of low temperature treatment. In the early mononuclear stage anthers with a floret length of ~5 mm were manually collected from plants of both temperature treatments and stored at −80°C for RNA isolation.⁴² The resulting product was reverse transcribed and amplified by polymerase chain reaction (PCR) to produce a cDNA sequencing library. The cDNA libraries were sent to Kang Cheng Biotechnology Co (Shanghai China). Solexa sequencing was performed on the HiSeq 2,000 platform (BGI).

After the Illumina sequencing, the raw sequences were first obtained in fastq format. Subsequently, the raw sequence passes through a filtering process to produce high quality readings. Unique sequences represented the sum of small RNA sequences of different base compositions, each unique small RNA contains several small RNA tags of identical sequence. All unique sequences were aligned against the rice genome to analyze their expression in the genome. Unique sequences recognized as tRNA, rRNA, snRNA, and snoRNA by searching through GenBank from the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/) and Rfam (11.0) where RNA family databases were dislodged.

The sequence is paired with the miRNA sequence of the species in miRBase (http://mirbase.org/), to obtain known miRNAs already present in the sample, and the base distribution. To identify miRNAs with differential expression between PA64S-H and PA64S-L, the expression of miRNAs was first quantified in transcripts per million. P value and fold change were calculated to confirm the significance of expression difference. The missing criteria for differential expression of miRNA was that | fold – change | > 2 and P value < .05. The sequences of known miRNA、exist miRNA and novel miRNA were used to Patmatch software to predict the potential target genes. After obtaining the target gene of miRNA expressed in a single sample, we performed GO (Gene Ontology) function analysis and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis on the part of the gene. In addition, we processed target gene enrichment analysis on the known miRNA、exist miRNA and novel miRNA.

Results

Sequencing and date analysis

To identify the differentially expressed miRNAs in an early stage of microspore development, two sRNA libraries were constructed and sequenced from PA64S-H and PA64S-L anthers at the early uninucleate stage, which produced about 14,030,132 and 16,844,884 raw reads, respectively. After removing sequences of low quality, adaptor contaminants, insert- null, RNAs smaller than 18 nucleotides, and polyA, 6,788,401 and 12,484,929 clear reads remained from PA64S-H and PA64S-L, respectively (Table 1). The clean reads were sequence-matched with the GenBank and Rfam databases, and the ncRNAs such as rRNA, tRNA, snRNA, and snoRNA, and the repeated sequences were filtered to obtain unannotated reads containing the miRNA (Table S1-S4).

Table 1. Statistics of small RNA sequences from the PA64S-H and PA64S-L libraries.

	PA64S-H		PA64S-L
Type	Count	Percentage (%)	Count	Percentage (%)
Total-reads	14030132	100%	16844884	100%
High-quality	12537263	89.3596%	15112119	89.7134%
3ʹadapter-null	136949	1.0923%	1052179	6.9625%
Insert-null	104738	0.8354%	10196	0.0675%
5ʹadapter-contaminants	127002	1.0130%	24331	0.1610%
Smaller-than-18 nt	5379902	42.9113%	1540181	10.1917%
PolyA	271	0.0022%	303	0.0020%
Clean-reads	6788401	54.1458%	12484929	82.6153%

The proportion of several non-coding miRNAs was low in PA64S-H and PA64S-L, among which rRNA accounted for the highest proportion and snRNA accounted for the least (Figure 1; Figure S1). The unannotated reads were aligned with the reference rice genome using miRDeep2 software to obtain positional information on the reference genome (Figure S2). In PA64S-H and PA64S-L, the number of miRNAs that could be aligned was 4,289,438 and 7,743,506, respectively, which accounted for about 60% of the total miRNAs (Table S5).The majority of sRNAs were 21–24 nt in both PA64S-H and PA64S-L libraries, which formed the typical size range for Dicer-derived products.⁴³ Among these, sRNAs with a length of 24 nt were the most abundant (30.37% in PA64S-H and 39.97% in PA64S-L), followed by those of 21 nt (16.53% in PA64S-H and 15.60% in PA64S-L) (Figure 2). The size distribution of sRNAs from the TGMS was quite similar to that of other plant species, such as A. thaliana⁴⁴ and Phaseolus vulgaris,⁴⁵ in which the 24 nt class dominates the sRNA transcriptome. However, there are some exceptions, as in Brassica juncea, where the 21-nt class dominates the sRNA in total and unique reads.⁴⁶

Figure 1. Summary of unique and total of small RNAs in non-coding RNAs in GenBank between PA64S-H and PA64S-L. a,b. Unique and total sRNAs in PA64S-H. c.d. Unique and total sRNAs in PA64S-L.

Figure 2. Length distribution of small RNAs in the different libraries. a. PA64-H. b. PA64-L.

Identification of known and new miRNAs

In the present study, a total of 263 known miRNAs (Table S6) were identified in the anthers of PA64S-H and PA64S-L at the early uninucleate stage. Among these known miRNAs, 151 were found in both sequence libraries, and 56 were specific to PA64S-H or PA64S-L (Figure 3(a)). The composition of miRNA families was slightly different in the two libraries, and there was a clear PA64S-H or PA64S-L-specific profile, for example, miR172 and miR158 was only expressed in PA64S-H, while miR169, miR4399, and miR9473 were only expressed in PA64S-L (Table S6). Further, sequencing frequencies of conserved miRNAs were also detected. Expression levels varied greatly between distinct miRNA families. Some of the miRNA family members, such as miR159, miR166, miR319, and miR168, were abundantly expressed in both libraries (Table S6), suggesting their conserved and essential roles for rice anther development. Compared to known reference data, these miRNAs are associated with pollen development. For example, miR159 regulates GAMYB or GAMYB-like expression, the gamyb deletion mutant is aborted due to the inability of the tapetum to perform PCD and degradation in rice;^35,36 miR159 and miR319 regulate genital development by regulating the activity of ARF6/8, including sepals, petals and anthers.^40,41 In contrast, other miRNA families, such as miR437, miR5509, miR5796, and miR171, were expressed at lower levels in both rice lines, suggesting the functional divergence among these miRNA families.

Figure 3. Venn diagram analysis of expressed genes of PA64S-H and PA64S-L. a. known miRNA. b. novel miRNA.

Generally, the first nucleotide of miRNAs is thought to be related with recruitment of the AGO protein. In Arabidopsis, AGO1 tends to recruit the miRNAs beginning with a uracil (U), while AGO2 and AGO5 show preferences for miRNAs initiating with adenosine (A) and cytosine (C), respectively.⁴⁷ In this study, the first nucleotide bias of known miRNAs was “U”, followed by “A” (Figure 4). The first nucleotide bias of known miRNAs from TGMS was quite similar to that of the cytoplasmic male sterile line and its maintainer line.³¹ This indicated that the identified miRNAs in rice anthers are conserved, and could be recruited by different AGO proteins and enter different regulatory pathways. After excluding sRNA reads homologous to known miRNAs and other noncoding RNAs from GenBank and the Rfam databases, many unannotated sequences that could not match any of the above databases were analyzed using miRDeep (http://deepbase.sysu.edu.cn/miRDeep.php) to predict novel miRNA candidates. A total of 321 putative new miRNAs were identified in the two libraries (Table S7), where 144 were present only in PA64-H, 114 only in PA64-L, and 63 in both sequencing libraries (Figure 3(b)). These novel miRNAs formed typical stem-loop hairpin structures (Figure S5), and their corresponding new sequences were mainly derived from the same precursor gene when aligned with the rice genome. This suggests that the unidentified small fragments of miRNAs may be rice-specific miRNAs. The first nucleotide bias and nucleotide bias at each position in the mature sequence of novel miRNAs was further analyzed (Figure S4) and was similar to that of known miRNAs.

Figure 4. The analysis of miRNA first nucleotide bias and nucleotide bias at each position in known miRNAs.

a. The miRNA first nucleotide bias analysis of known miRNAs in PA64S-H.b. The miRNA first nucleotide bias analysis of known miRNAs in PA64S-L.c. The miRNA nucleotide bias at each position analysis of known miRNA in PA64S-H.d. The miRNA nucleotide bias at each position analysis of known miRNA in PA64S-L.

Differential expression of miRNAs during polloen development of PA64S-H and PA64S-L

To identify miRNAs related to anther development and regulation of male-sterility in PA64S under high and low temperature, a total of 263 miRNAs were used for differential expression analysis between

PA64S-H and PA64S-L (Table S6), and the normalized expression values of these miRNAs were used to draw the scatter plot (Figure 5(a)). The selection threshold were fold-change |(log2)| ≥2 and P-value < 0.05. Considering that extremely low abundances might be false results, the miRNAs with less than 15 norm reads were removed from expression. Among these differentially expressed, known miRNAs, 59 were up-regulated and 74 were down-regulated in PA64S-L compared to those in PA64S-H (Table S15), while other miRNAs were equally expressed (Figure 5(a)). Studies have shown that miR390 and miR164 are significantly differentially expressed in the photosensitive male sterile mutant 7B-1 and its wild type,³⁰ and low temperature stress leads to the induction of miR171 and miR172 in Brachypodium serrata and Arabidopsis,⁴⁸ respectively. To validate the differential expression data obtained from the high-throughput sequencing, we further tested the expression patterns of six selected, known miRNAs by quantitative real-time RT-PCR analysis. The expression levels of MIR531, MIR6300, and MIR9773 decreased, whereas MIR1118, MIR2275, and MIR9473 were up-regulated in PA64S-H compared with that of PA64S-L (Figure 5(b)). The qRT-PCR results were consistent with the high-throughput sequencing, implying that the high throughput sequencing results were reliable.

Figure 5. Comparison and validation of the differential expression of known miRNA in the different libraries. a. Scatter plot of expression rations for known in Pa64S-H and PA64S-L. Red points represent miRNAs with fold change>2, showing upregulated miRNAs in PA64S-L; blue points represent miRNAs with 1/2≤ fold change ≤2, indicating equally expressed miRNAs in both libraries; green points represent miRNAs with fold change<1/2, meaning downregulated miRNAs in PA64S-L. Fold change = normalized expression in PA64S-L/normalized expression in PA64S-H. b. Comparison of the miRNA expression levels between PA64S-H and PA64S-L by qPCR for MIR531, MIR6300, MIR9773, MIR1118, MIR2275, MIR9473. The experiments were performed with three biological replicates, and the error bars represent the standard error.

Target prediction and function analysis

Given that these miRNAs found during the reproductive stage may have functions associated with anther development, subsequent analyzes focused on their targets. A total of 5944 target genes were identified in 769 miRNAs (Table S10), including known miRNAs corresponding to 2226 target genes (Table S11), and novel miRNAs correspond to 233 target genes (Table S12). The number of target genes corresponding to novel miRNAs was significantly lower than that of conserved miRNAs. Most novel miRNAs did not predict any target gene, and some new miRNAs corresponded to the same target gene as the known miRNAs. For example, the target gene corresponding to novel-m0280-5p is ABH1 transcription factor, and the target gene corresponding to novel-m0098-5p is F-box domain containing protein (Table S12). This reflects, to some extent, that the functions of these novel miRNAs may be similar to those of conserved miRNAs, and regulate the development of rice flowers by regulating similar target genes.

Among the prediction targets, some miRNA targets related to anther development are worth mentioning. For example, growth-regulating factors (OSGRF1/2/3/6/7/9/10/12), kinases, MYB transcription factors, and TCP transcription factors have potential roles in anther development. API5, receptor-like kinase, and bHLH proteins are involved in tapetum development/degeneration^{31,36,41,49–52} (Table S10). In addition, 12 putative auxin response factors (ARFs) were predicted to be targeted by miRNAs (Table S10), which are quite similar to that mentioned in a previous study,²⁹ and indicate that the miRNA-mediated silencing of the ARF gene has a potential role in pollen development. This result suggests that miRNAs can directly regulate gene expression at both transcription and translation level.

To gain more information on miRNA targets for TGMS at high and low temperature, the targets of differentially expressed miRNAs were further analyzed. Interestingly, many differentially expressed miRNA targets have been shown to be associated with cold or heat stress responses, as well as with anther development. MIR5809-y, which specifically expressed in PA64S-L, is predicted to target transcripts encoding a receptor-like protein kinase (Os01G0155500) (Table S15; Table S10). Mutations in receptor-like protein kinase have been reported to result in male sterility in Arabidopsis.⁵³ The target gene bHLH (Os02g0116600), corresponding to MIR3440-y, is responsive to the regulation of cold stress in rice and associated with anther development^10,54 (Table S15; Table S10). In addition, many heat shock proteins, which are the targets of MIR164-y, MIR419-x, and MIR8011-x, are responsible for protein folding, assembly, and translocation, and play pivotal roles in protecting plants against stress by stabilizing proteins and membranes^55,56 (Table S15; Table S10). These significantly differentially expressed miRNAs are specifically expressed in PA64S-H (Table S15), and their target gene expression was downregulated under heat stress (Figure R2).⁵⁷

GO terms and KEGG pathway enrichment analysis of targets

To illuminate the function of the differentially expressed miRNAs, gene ontology (GO) analysis and KEGG pathway annotation were used to enrich and classify the given target genes into specific informative groups.⁵⁸ Here, the most frequent “biological process” terms were ‘‘metabolic process’’ and ‘‘cellular process’’ (Figure S3). “Cell”, “cell part”, and “organelle” were the top three “cellular components” (Figure S3), and the most abundant terms in ‘‘molecular function’’ were ‘‘binding’’, ‘‘catalytic activity’’, ‘‘transcription regulator activity’’, and ‘‘transcription activity’’ (Figure 6, Figure S3). Further, the targets of differential expression of known miRNAs were analyzed. The result showed that the targets were classified into 22 reference terms of “biological processes”, with the most frequent term being ‘‘metabolic process’’, followed by ‘‘cellular process’’, ‘‘response to stimulus’’, and “biological regulation” (Figure 6(a)). The results indicate a correlation between pollen development and low temperature. Comparatively, the results for ‘‘biological process’’ were classified into 18 reference terms in novel miRNAs, with the most common terms being ‘‘cellular process’’, ‘‘metabolic process’’, ‘‘response to stimulus’’, and ‘‘localization’’ (Figure 6(b)).

Figure 6. Gene ontology analysis for the targets of the differential expression of known and novel miRNAs based on the terms of cellular component, biological process and molecular function in PA64S-H and PA64S-L. a. known miRNAs. b. novel miRNAs.

To further illustrate our findings, the KEGG enrichment networks were also used for comparative analysis between differentially expressed target genes in PA64S-H and PA64S-L (Table S8, Table S9). As a result, most of the target genes corresponding to differentially expressed miRNAs were involved in metabolic pathways and biosynthesis of secondary metabolites (Table 2). Among them, several pathways are especially noteworthy. First, 13 target genes were enriched in the starch and sucrose metabolism pathway (ko00500) (Table 2, Table S8, Table S9). Starch and sucrose metabolism provide raw materials for all stages of plant growth, and the accumulation of soluble sugar is directly related to pollen fertility.⁵⁹ Second, differentially expressed miRNA target genes were also heavily involved in sphingolipid metabolism (ko00600) (Table 2, Table S8, Table S9). Sphingolipids are key cellular membranes and signaling molecules involved in several cellular activities.^60–62 Sphingolipid homeostasis is also reported to be involved in responses to ABA and cold temperatures,^63,64 as well as regulation of PCD and pollen fertility.^65,66 This indicates that miRNAs affect pollen fertility under low temperatures indirectly by silencing sphingolipid metabolism-related genes. Third, arginine and proline metabolism pathways were also annotated (ko00330) (Table S8). Retardation of glutamic and aspartic acid conversion to proline is considered to be a cause of pollen sterility in the male sterile line.^67,68 In addition, the plant hormone signal transduction pathway was also enriched (ko04075) (Table S8), suggesting that miRNAs are involved in the regulation of plant growth and development and respond to environmental signals by affecting important transcription factors in the hormone signaling pathway.⁶⁹

Table 2. Some KEGG enrichment pathways analysis.

KEGG Pathway	Target genes with pathway annotation (272)	All genes of the species with pathway annotation (5086)	P value	Q value	Pathway ID
Metabolic pathways	106 (38.97%)	1979 (38.91%)	0.5155038	0.9952613	ko01100
Biosynthesis of secondary metabolites	52 (19.12%)	1050 (20.64%)	0.7611412	0.9952613	ko01110
Pyruvate metabolism	14 (5.15%)	125 (2.46%)	0.006478281	0.105095012	ko00620
Phenylpropanoid biosynthesis	14 (5.15%)	157 (3.09%)	0.04023822	0.22131021	ko00940
Starch and sucrose metabolism	13 (4.78%)	196 (3.85%)	0.2485438	0.781137657	ko00500
Sphingolipid metabolism	7 (2.57%)	32 (0.63%)	0.001227664	0.036011477	ko00600
Plant hormone signal transduction	6 (2.21%)	254 (4.99%)	0.9946631	0.9952613	ko04075
Endocytosis	12 (4.41%)	140 (2.75%)	0.07010644	0.308468336	ko04144
Purine metabolism	12 (4.41%)	165 (3.24%)	0.9946631	0.56361197	ko00230
Arginine and proline metabolism	1 (0.37%)	91 (1.79%)	0.993576	0.9952613	ko00330

Discussion

Peiai64S (O. sativa L) is a photoperiod temperature-sensitive genic male sterile line (PTGMS). Its fertility alternates between sterile and fertile at a certain day length and temperature, and has been widely used in hybrid seed production.⁶ The photosensitivity of PA64S is regulated by a long-chain non-coding RNA LDMAR. A SNP (G→C) mutation on LDMAR alters the secondary structure of LDMAR, reducing the transcriptional level of LDMAR under long-day conditions, which leads to premature PCD in pollen and results in photosensitive male sterility in rice.⁷⁰ Zhou revealed the TGMS mechanism of RNase Z^S1-mediated UbL40 mRNA processing in Annon S-1.⁷¹ So far, what confers the TGMS traits in PA64S is still unknown.

It has been reported that the “miRNA-target gene” regulatory pathway regulates the fertility of plants. For example, in Arabidopsis, miR167 regulates its target ARF transcription factor, miR159 regulates the target MYB transcription factor, and miR319 regulates the TCP transcription factor to participate in anther development and maturation.^41,72,73 However, the relationship between miRNA and TGMS in rice has not been reported. In this study, we used high throughout sequencing technology to profile miRNAs expression and target gene function in PA64S anther under high and low temperatures. The result showed that a total of 263 known miRNAs was identified in the anthers of PA64S-H and PA64S-L at the early uninucleate stage (Table S6). As a result, many known miRNAs are abundantly expressed in two libraries (Table S6). For example, MiR159 and miR319 regulate the development of sepals, petals and anthers by modulating the activity of the target gene ARF6/8.^40,41 The results suggested that abundantly expressed known miRNAs are conserved and play important roles in rice anther development. Among the 263 known miRNAs, 133 differentially expressed miRNAs were identified in PA64S-L compared with PA64S-H, of which 59 were up-regulated and 74 were down-regulated (Table S15). Interestingly, many male fertility-related miRNAs, including miR390 and miR164, were significantly, differentially expressed.^30,74 Furthermore, miR172, which was only expressed in PA64S-H, is not only induced by low temperature stress, but also regulates male and female differentiation of floral organs^48,75 (Table S15).

Targets of these differentially expressed miRNAs were also analyzed. Growth-regulating factors (GRFs), kinases, MYB transcription factors, and TCP transcription factors were characterized in this study (Table S10). OsGRF1 not only regulates leaf growth, but also controls flowering time.⁵⁰ Overexpression of AtGRF1 and AtGRF2 in Arabidopsis also delayed flowering compared to the wild type.⁵¹ Kinases act as major regulators of various pathways. Mutations in receptor-like protein kinases have been reported to result in male sterility in Arabidopsis.⁵³ MIR5809-y was specifically expressed in PA64S-H (Table S15) and is predicted to target a transcript encoding a receptor-like protein kinase (Os01G0155500) (Table S10). In Arabidopsis, the target MYB and TCP transcription factors corresponding to miR159 and miR319, respectively, are involved in the development and maturation of anthers.^41,73 To further illuminate the functional role of the targets, the predicted targets were analyzed with GO and KEGG. Some metabolic and biosynthetic pathways that may be associated with anther development, include biosynthesis of secondary metabolites, the phenylpropanoid biosynthesis pathway, starch, sucrose, and sphingolipid metabolism. In addition, the plant hormone signal transduction pathway is also enriched, which is associated with anther development and stress responses.^32,35,37,76 In addition, we used CRISPR-Cas9 technology to knock out target genes for some clearly differentially expressed miRNAs. Mutation phenotypic analysis showed that pollen fertility was transformed under high temperature and low temperature (data not show). These results suggest that the miRNA-target gene pathway may be involved in regulating pollen fertility at high or low temperatures.

Conclusions

We employed high-throughput sequencing to identify expression profiles of anther miRNAs from PA64S-H and PA64S-L, and found some clues that may be involved in pollen fertility under low temperature in PA64S. These findings may shed light on the possible mechanism of miRNAs regulating TGMS occurrence in rice during anther development. However, further research is needed to verify the miRNA/target gene function role by using the CRISPR-Cas9 technology to knock out for validation analysis. In addition, the functions of a large number of differentially expressed miRNAs/target genes are still unknown, which is an obstacle to clarify the fertility mechanism in PA64S.

Disclosure of potential conflict of interest

The authors declare no conflict of interest.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

This work was supported by grants of NSFC [31570316; 31371244], Projects of Hunan Provincial Education Department [17C0974], and Hunan Provincial Innovation Foundation for Postgraduate [CX2017B230, CX2018B300].