Comparative Transcriptomes Profiling of Photoperiod-sensitive Male Sterile Rice Nongken 58S During the Male Sterility Transition between Short-day and Long-day

Background Photoperiod-sensitive genic male sterile (PGMS) rice, Nongken 58S, was discovered in 1973. It has been widely used for the production of hybrid rice, and great achievements have been made in improving rice yields. However, the mechanism of the male sterility transition in PGMS rice remains to be determined. Results To investigate the transcriptome during the male sterility transition in PGMS rice, the transcriptome of Nongken 58S under short-day (SD) and long-day (LD) at the glume primordium differentiation and pistil/stamen primordium forming stages was compared. Seventy-three and 128 differentially expressed genes (DEGs) were identified at the glume primordium differentiation and pistil/stamen primordium forming stages, respectively. Five and 22 genes were markedly up-regulated (≥ 5-fold), and two and five genes were considerably down-regulated (≥ 5-fold) under SD during the male sterility transition. Gene ontology annotation and pathway analysis revealed that four biological processes and the circadian rhythms and the flowering pathways coordinately regulated the male sterility transition. Further quantitative PCR analysis demonstrated that the circadian rhythms of OsPRR1, OsPRR37, OsGI, Hd1, OsLHY and OsDof in leaves were obviously different between Nongken 58S and Nongken 58 under LD conditions. Moreover, both OsPRR37 and Hd1 in the inflorescence displayed differences between Nongken 58S and Nongken 58 under both LD and SD conditions. Conclusion The results presented here indicate that the transcriptome in Nongken 58S was significantly suppressed under LD conditions. Among these DEGs, the circadian rhythm and the flowering pathway were involved in the male sterility transition. Furthermore, these pathways were coordinately involved in the male sterility transition in PGMS rice.


Background
Photoperiod-sensitive genic male sterile (PGMS) rice, Nongken 58S, was discovered as a spontaneous mutant in the japonica rice cultivar Nongken 58 (Oryza sativa ssp. japonica) grown in Hubei Province, China in 1973 [1]. It has since been used for production of hybrid rice, and during the past two decades, great achievements have been made in improving rice yields in China using two-line hybrid rice. Several important features of Nongken 58S have been characterised. Its fertility is highly regulated by day length at specific inflorescence developmental stages. Complete male sterility can be induced when the day length is greater than 14 h from the glume primordium differentiation stage to the pollen mother cell forming stage. However, male fertility returns gradually when the day length is shorter than 14 h [1]. This is called the second photoperiod phenomenon [2], referring to the short day photoperiod that rice requires for the transition from vegetative growth to reproductive growth [2]. This agronomic trait is genetically controlled by a recessive locus within the nuclear genome. However, the genetic loci vary depending on the genetic background of the recipient parents. That is, when crossed with a less photoperiod-sensitive cultivar, the genes for the trait appear as two or three loci [3][4][5][6][7]. The photoperiod sensitivity can also be affected by temperature when the pms gene is crossed into a temperature-sensitive cultivar. Previous studies have indicated that phytochromes and cryptochromes are involved in the male sterility transition [2]. Because male sterility is highly and coordinately regulated by day length and temperature, this leads to difficulty in the accurate identification of the male sterile phenotype in the segregation populations between japonica and indica under natural conditions. Although it would be of great value to breeders, this makes mapping the pms gene challenging [3][4][5]8,9]. It is important to understand the molecular mechanism of the male sterility transition and to explore whether the light or circadian rhythm signal transduction pathway is involved in this process. It may also be useful to understand how light or the circadian rhythm regulates microsporocyte development. This information could assist breeders in selecting recipient parents for molecular breeding programs.
Circadian rhythms, in general, take the form of sinusoidal waves that can be described using mathematical terms such as period, phase and amplitude [10][11][12]. The expression of several genes has been found to be associated with the circadian rhythm in Arabidopsis. These include TIMING OF CAB EXPRESSION1 (TOC1) [13], LATE ELONGATED HYPOCOTYL (LHY) [14], CIRCA-DIAN AND CLOCK ASSOCIATED1 (CCA1) [15] and GIGANTEA (GI). TOC1, LHY, CCA1 and unknown factor Y [16,17] comprise interlocked transcriptional feedback loops. These feedback loops play important roles in the plant central clock. These loops integrate environmental factors, such as light and temperature, into the central clock through the input signaling pathway and import the rhythm signal into downstream signaling pathways through output signaling pathways.
This paper reports the comparison of the transcriptomes of Nongken 58S under SD and LD conditions during the male sterility transition. The repressive expression profile under LD conditions identified 183 differentially expressed genes (DEGs). Gene ontology (GO) and pathway analysis of the DEGs revealed that the circadian rhythm and the flowering pathways were involved in the male sterility transition. In further analysis, qPCR results indicated that the circadian rhythms of OsPRR1, OsPRR37, OsGI, Hd1, OsLHY and OsDof in Nongken 58S were significantly different from those in Nongken 58. This suggests that the circadian rhythm and the flowering pathways were coordinately involved in regulation of the male sterility transition.

Plant materials
The varieties Nongken 58 and the mutant Nongken 58S (Oryza sativa L. ssp. japonica) were used for this study. All plants were grown under natural conditions at Wuhan University Campus. SD (10-h light/14-h dark) treatment was conducted when the seedlings had greater than five leaves. Differentiation from vegetative growth to reproductive growth was promoted after 10-12 days of SD treatment. When the inflorescences had developed to the secondary branch differentiation stage, the rice plants were divided into two groups for LD (15-h light/9-h dark) and SD treatment, respectively.

RNA isolation and cDNA synthesis
For total RNA isolation, the second fully expanded leaf was harvested at two developmental stages: glume primordium differentiation and pistil/stamen primordium formation. For diurnal expression profile analysis, the leaves were harvested at 4 h intervals over 1 day, frozen immediately in liquid nitrogen, and stored at -80°C until use [39]. The samples are annotated in Additional file 1, Table S1.
A TRIzol Reagent Kit (Invitrogen, Carslbad, USA) was used for total RNA isolation following the manufacturer's instructions. Total RNA was treated with RNase-free DNase I (New England Biolabs, Hitchin, UK) to remove DNA contamination before cDNA synthesis. Two micrograms of total RNA and Oligo (dT) were used as template and primers for first strand cDNA synthesis by M-MLV reverse transcriptase (Promega, Madison, USA).

Microarray analysis procedure and comparison strategy
Two micrograms of total RNA were used for doublestranded cDNA synthesis, and biotin-tagged cRNA was prepared using a MessageAmp™ II cRNA Amplification Kit according to the manufacturer's instructions. The resulting bio-tagged cRNA was fragmented into strands of 35-200 bases according to Affymetrix protocols. The fragmented cRNA was hybridised to an Affymetrix Gen-eChip Rice Genome Array containing 48,564 and 1,260 transcripts representing the japonica and indica cultivars, respectively. The microarray and data analysis were contracted to CapitalBio Corporation (Beijing, China). The hybridisation was performed at 45°C with rotation for 16 h (GeneChip Hybridization Oven 640, Affymetrix). The GeneChip arrays were washed and then stained (streptavidin-phycoerythrin) on an Affymetrix Fluidics Station 450. Scanning was conducted using a GeneChip Scanner 3000. All chip data represented biological triplicates.
The scanned images were examined by visual inspection and then processed to generate raw data using the default setting of GeneChip Operating Software (GCOS 1.4). dChip software was used to perform invariant-set normalization to normalize the arrays according to dChip users' manual. 5% of Perfect-Match (PM) probe signals were used as background removal. Annotation of the samples is shown in Additional file 1, Table S1. The samples were compared as follows: G-SSD2 vs. G-SLD2 and P-SSD2 vs. P-SLD2. The criteria for determining up-and down-regulated genes were fold changes (FCs) of ≥ 2 and ≤ 0.5, respectively. The complete microarray data have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO series accession number GSE29820.

Cluster, gene ontology annotation and pathway analysis
All microarray data were analysed using Significant Analysis of Microarray (SAM) 3.02 software. Genes with FC ≥ 2 or FC ≤ 0.5 were chosen for the t-test, and genes with P values < 0.05 were chosen for further analysis. Cluster analysis was performed using Cluster 3.0 software. GO annotation was performed using the GeneOntology Enrichment Analysis Software Toolkit (http:// omicslab.genetics.ac.cn/GOEAST/php/affymetrix.php) [40,41]. GOs with P values < 0.01 were selected. Pathway analysis was accomplished using the MAS 3.0 system (http://www.capitalbio.com/index.asp) (CapitalBio Corporation, Beijing, China) and by searching the Kyoto Encyclopedia of Genes and Genomes, BioCarta and GenMAPP databases. Pathways with P values < 0.01 were selected.

Quantitative PCR analysis
For quantitative PCR, the cycling conditions were 95°C for 10 min, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s, and then a melting curve from 60°C to 95°C in 0.3°C increments. Negative controls without the template DNA were used to ensure that primer dimers were not interfering with amplification. qPCR data was captured and analysed using StepOne software (ver. 2.0). The relative gene expression levels in each cDNA sample were obtained by normalisation to ACTIN1 using the formula 2 -(CT gene -CT actin1) [42]. The primers used in this study are listed in Additional file 2, Table S2. All qPCR data is the result of biological triplicates.

Results
Dynamic transcriptome profiling of photoperiod treatment during the male sterility transition As shown in Figure 1, pollen abortion was induced in Nongken 58S under LD conditions (15-h light/9-h dark) from the glume primordium differentiation to the pollen mother cell forming stages ( Figure 1A-C) [2,43]. Under SD conditions (10-h light/14-h dark), the pollen partially regained fertility ( Figure 1E, G). By contrast, pollen fertility in Nongken 58 was not affected ( Figure 1D, F). These results indicated that the fertility of PGMS could not be completely restored using normal photoperiod conditions, implying that some defective signaling pathway may be constitutively affected. To investigate the transcriptome changes during the male sterility transition process in PGMS rice, we chose the leaves, light signal receiver, as experimental material at photoperiodsensitive stage. Based on preliminary experiments, the diurnal expression profiles of those genes peaked at 02:00 h (data not shown). Thus, 02:00 h was selected as the expression profile time point for the microarray. Affymetrix GeneChips were used for global transcriptome profiling analysis. The correlation coefficients between different biological repeats were calculated, and the result indicates that the profilings were reproducible (Additional file 3, Table S3). A total of 20,444 and 19,964 transcripts were detected P values < 0.05 in the comparison of G-SSD2 vs. G-SLD2 and P-SSD2 vs. P-SLD2, respectively. In comparing the SD and LD conditions at the two developmental stages, 183 DEGs with two fold changes were detected P < 0.05. The DEGs were categorised using Cluster 3.0 software (Figure 2A, B, C). The number of DEGs with FC ≥ 2 or ≤ 0.5 increased significantly as inflorescence development proceeded ( Figure 2D, E; Additional file 4, Table  S4; Additional file 5, Table S5). At the glume primordium differentiation and the pistil/stamen primordium forming stages, 73 and 128 DEGs with FC ≥ 2 or ≤ 0.5 (P < 0.05) were detected, respectively ( Figure 2D). Five genes were identified as markedly up-regulated (≥ 5fold) and two genes were noted as considerably downregulated (≥ 5-fold) under SD conditions at the glume primordium stage ( Figure 2E; Additional file 4, Table  S4). Twenty-two genes were considerably up-regulated (≥ 5-fold), and only five genes were down-regulated (≥ 5-fold) at the pistil/stamen primordium forming stage ( Figure 2E, Additional file 5, Table S5). These results indicated that more genes were suppressed than were activated under LD conditions. This suggests that transcriptome profiling in Nongken 58S was significantly suppressed, and the repression effects were markedly intensified as inflorescence development proceeded. This is consistent with previous studies in which pollen fertility was found to decrease as LD treatment increased [44].

Defining DEGs for GO annotation and pathway analysis
To annotate the DEGs, the sequences were uploaded to the GOEAST website for GO annotation analysis [40,41,45]. The results indicated that some biological processes were significantly affected during the second photoperiod under LD conditions. Four GO biological processes related to the regulation of the transition from the vegetative to the reproductive phase (GO:0048510), the regulation of SD photoperiodism (GO:0048572), the photoperiodism (GO:0009648), and inflorescence development (GO:0010229) were investigated during the second photoperiod ( Table 1). The results implied that these processes were related to the male sterility transition in PGMS rice.
To define the potential signal transduction pathways, the DEGs were uploaded to the MAS 3.0 system. Five pathways and seven pathways were identified from the transcriptome data for G-SSD2 vs. G-SLD2 and P-SSD2 vs. P-SLD2, respectively ( Table 2). The circadian rhythm pathway was top ranked (P = 6.57E-06 and 1.38E-05) in both inflorescence developmental stages ( Table 2). The results indicated that three genes, OsPRR1 (homologous to PRR1/TOC1 in Arabidopsis), OsLHY (homologous to LHY in Arabidopsis) and MYB transcription factor (LOC_Os06g51260), belong to the circadian rhythm pathway. This suggests that the circadian rhythm may be a principal signal transduction pathway in the male sterility transition in PGMS.

Circadian rhythm involvement in the male sterility transition
As mentioned previously, several clock genes [OsPRR1 (LOC_Os02g40510), OsPRR37 (LOC_Os07g49460), OsGI (LOC_Os01g08700), OsLHY (LOC_Os08g06110) and MYB transcription factor (LOC_Os06g51260)] were identified (Additional file 4, Table S4, Additional file 5, Table  S5). Among these genes, OsPRR37 and OsGI were upregulated 3.0-fold and 8.3-fold, respectively; OsLHY and the MYB transcription factor were down-regulated 2.4fold and 2.6-fold, respectively. No difference in OsPRR1 expression was detected in the G-SSD2 vs. G-SLD2 data (Additional file 4, Table S4); in the P-SSD2 vs. P-SLD2 data, OsPRR1 was up-regulated 2.0-fold and the MYB transcription factor was down-regulated 4.4-fold. No differences in OsPRR37, OsGI and OsLHY were detected at this inflorescence developmental stage. The fold changes indicate that these clock genes were activated or repressed under LD conditions in Nongken 58S, implying that they may participate in the male sterility transition.
Because gene transcription is dynamic, the microarray data reflect expression differences at 02:00 h only. Therefore, the diurnal expression profiles of these genes were analysed in both Nongken 58S and Nongken 58 at the glume primordium differentiation and pistil/stamen primordium forming stages under LD and SD conditions. OsPRR1, OsPRR37 and OsGI exhibited the similar circadian rhythm patterns under LD conditions ( Figure  3A-C). For three of these genes, the differences between Nongken 58S and Nongken 58 only occurred at 18:00 h at the glume primordium differentiation stage. The expression levels in Nongken 58S were significantly higher than in Nongken 58. No difference was found at the pistil/stamen primordium forming stage under LD and SD conditions ( Figure 3A-C, Additional file 6, Figure S1A-C). These results indicate a circadian rhythm difference in the expression of the three genes between Nongken 58S and Nongken 58 [12].
Taken together, these differences only occurred under LD conditions at the photoperiod-sensitive developmental stage. Based on the genetic background of Nongken 58S and Nongken 58, these genes were all related to the male sterility transition and the rice circadian rhythm was involved in this process.  Table S4; Additional file 5, Table S5). Ehd1, Hd3a, RFT1 are crucial genes for rice flowering [18][19][20][21]23], and OsMADS1 plays key roles in specifying floral organ and meristem identity in rice [46][47][48]. Ehd2 and Hd1 are also important components of the flowering pathway. To determine whether these genes were involved in the male sterility transition, the diurnal expression profiles of these genes were analysed. Hd1 was the only gene exhibiting different daily expression profiles between Nongken 58S and Nongken 58 under LD conditions ( Figure 3D and Additional file 7, Figure  S2). The main differences occurred at the pistil/stamen primordium forming stage. The transcription levels at 22:00 h and 06:00 h in Nongken 58S were significantly different from those in Nongken 58, indicating that the peak phase of Hd1 was at 06:00 h in Nongken 58 and it was shifted forward to 22:00 h in Nongken 58S. The phase is a mathematical term of plant circadian rhythms [12]. The phase-shift indicated that the circadian rhythm of Hd1 in Nongken 58S was notably different from that in Nongken 58, while the flowering time of Nongken 58 is similar to that of Nongken 58S under LD conditions.  This suggests that Hd1 plays a role in the male sterility transition in PGMS.
OsPRR37 and Hd1 were differentially expressed in the inflorescences of Nongken 58S and Nongken 58 Because of light signals received from the leaf, knowledge of the expression profiles of DEGs in the inflorescence is essential. Further qPCR of the expression patterns of these genes in the inflorescence showed that OsPRR1, OsGI, OsLHY, OsDof expression was similar in Nongken 58S and Nongken 58 under LD and SD conditions (Additional file 8, Figure S3). Only OsPRR37 and Hd1 were differentially expressed in the inflorescences of Nongken 58S and Nongken 58 at the pistil/stamen primordium forming stage ( Figure 4). The expression level of OsPRR37 in Nongken 58S was significantly lower than that in Nongken 58 at 17:00 h under LD conditions ( Figure 4A), and the level of OsPRR37 was significantly lower than that in Nongken 58 at 20:00 h under SD conditions ( Figure 4B). The expression level of Hd1 in Nongken 58S was significantly greater than that in Nongken 58 at 20:00 h and 02:00 h under LD conditions ( Figure 4C), but was significantly lower than that in Nongken at 20:00 h. However, the transcription level in Nongken 58S at 02:00 h was higher than that in Nongken 58. These results suggest that OsPRR37 and Hd1 are likely involved in the male sterility transition directly under LD conditions in the inflorescence. under LD and SD conditions. The y-axis represents the expression ratio relative to the endogenous reference Actin1 gene using qPCR. NLD and SLD represent Nongken 58 and Nongken 58S under LD conditions. NSD and SSD represent Nongken 58 and Nongken 58S under SD conditions. P represents pistil/stamen primordium forming stages. Asterisks show statistically significant differences (* = P < 0.05).

Figure 3
Diurnal expression profiles of OsPRR1, OsPRR37, OsGI, Hd1, OsLHY and OsDof The diagram shows the circadian rhythms of OsPRR1, OsPRR37, OsGI, Hd1, OsLHY and OsDof under LD conditions. The y-axis represents the expression ratio relative to the endogenous reference Actin1 gene using qPCR. The boxes on the x-axis represent day (blank) and night (solid). NLD and SLD represent Nongken 58 and Nongken 58S under LD conditions. G and P indicate the glume primordium differentiation and pistil/stamen primordium forming stages, respectively. Asterisks show statistically significant differences (* = P < 0.05; ** = P < 0.01).
Moreover, the differences under SD conditions suggest a potential explanation for why pollen fertility cannot be completely recovered under SD in Nongken 58S.

Discussion
Circadian rhythm and light signals have been reported to be involved in several physiological processes, including flowering, various stress responses and metabolism [26,30,31,38]. However, this is the first report suggesting a role for circadian rhythm and flowering signals in the regulation network of male sterility. In general, these genes promote the transition from vegetative growth to reproductive growth under SD conditions in SD plants and under LD conditions in LD plants [23]. Although the male sterility transition in PGMS rice is known to be regulated by LD conditions, the mechanism of regulation remains unknown. Global gene expression profiling under LD conditions revealed a suppressive trend and categorised into the circadian rhythm and flowering pathways. Both pathways regulate the process in a coordinated manner. Further studies revealed that the diurnal expression profiles of OsPRR1, OsPRR37, OsGI, OsLHY, OsDof and Hd1 under LD conditions were reprogrammed.
In this study, OsPRR1, OsPRR37 and OsGI exhibited similar expression patterns in the leaf under LD conditions, suggesting that these three genes may function together. OsLHY and OsDof also exhibit similar differentially expressed patterns, suggesting that they may function coordinately. In Arabidopsis, TOC1/PRR1, LHY and GI, key component factors in the central clock, comprise interlocked transcriptional feedback loops to regulate plant circadian rhythms [12]. In PGMS rice, OsPRR1, OsLHY and OsGI were not only involved in the circadian rhythm, but were also involved in the male sterility transition in coordination with OsPRR37 and OsDof. The LD signals maybe integrated and transmitted to the downstream genes through rice central clock in leaves. The day length signal may have been transmitted to the inflorescence via Hd1 and OsPRR37, ultimately leading to pollen abortion.
In the rice flowering pathway, OsGI is upstream of Hd1 and is positively correlated with Hd1 expression. Hd1 is activated by OsGI under LD or SD conditions. Hd1 activates Hd3a/RFT1 expression under SD conditions to promote flowering. However, Hd1 suppresses Hd3a/RFT1 expression under LD conditions to inhibit flowering [24]. In PGMS rice, Hd1 expression may also have a positive correlation with OsGI under LD conditions in Nongken 58S, as its expression was activated after OsGI expression (Figure 3 C, D). This implies that Hd1 may function downstream of OsGI in the male sterility transition in PGMS rice. In this case, Hd3a/ RFT1 may not be the downstream regulator of Hd1 since there were no rhythm differences between Nongken 58S and Nongken 58. Hd1 was activated by OsGI under LD conditions and suppressed the unknown factors. This may influence downstream gene expression and finally affect male fertility.
In leaves, the diurnal expression patterns of OsPRR1, OsPRR37, OsGI, OsLHY, OsDof and Hd1 showed significant differences between Nongken 58S and Nongken 58. In the inflorescence, only OsPRR37 and Hd1 exhibited different expression patterns between Nongken 58S and Nongken 58 under both LD and SD. Therefore, it is speculated that the clock genes in leaves may function as sensors for day length. These genes may receive day length signals and integrate and transmit the signals into the inflorescence through a series of unknown pathways. Both OsPRR37 and Hd1 may be the effectors in the inflorescence, and they are likely involved in regulation of male sterility directly.
Previous studies have indicated that phytochromes and cryptochromes are involved in the male sterility transition. However, there were no diurnal expression differences in OsPhyA, OsPhyB, OsCry1a, OsCry1b, OsCry2 and OsCry3 under LD conditions between Nongken 58S and Nongken 58 (data not shown). Using an ELISA assay, Wang, et al. found that phyA content in Nongken 58S leaves was higher than in Nongken 58 under the identical day length treatment [49]. The differences in phytochromes and cryptochromes may be at the protein level and needs further investigation.

Conclusion
In conclusion, the transcriptome of Nongken 58S was significantly suppressed and the repression effects were markedly intensified when inflorescence development proceeded. Rice circadian rhythm genes OsPRR1, OsPRR37, OsGI, OsLHY and OsDof and the flowering gene Hd1 were coordinately involved in signal transduction in leaves. Furthermore, Hd1 and OsPRR37 may be signal sensors in inflorescences to directly affect the male sterility transition. Additional file 5: Table S5 Differentially expressed genes in P-SSD2 vs. P-SLD2 data. All the DEGs with FC ≥ 2 or FC ≤ 0.5 are listed in this table.
Additional file 6: Figure S1 Diurnal expression profiles of OsPRR1, OsPRR37, OsGI, Hd1, OsLHY and OsDof genes in leaf. The diagram shows the diurnal expression profiles of OsPRR1, OsPRR37, OsGI, Hd1, OsLHY and OsDof genes under SD conditions. The y-axis represents the expression ratio relative to the endogenous reference Actin1 gene using qPCR. The boxes on the x-axis represent day (blank) and night (solid). NSD and SSD represent Nongken 58 and Nongken 58S under SD conditions. G and P indicate the glume primordium differentiation and pistil/stamen primordium forming stages, respectively. Additional file 7: Figure S2 Diurnal expression profiles of Ehd2, Ehd1, Hd3a, RFT1 and OsMADS1 under LD conditions in leaf. The diagram shows the circadian rhythm of the expression of Ehd2, Ehd1, Hd3a, RFT1 and OsMADS1 under LD conditions. The y-axis represents the expression ratio relative to the endogenous reference Actin1 gene using qPCR. The boxes on the x-axis represent day (blank) and night (solid). NLD and SLD represent Nongken 58 and Nongken 58S under LD conditions. G and P indicate the glume primordium differentiation and pistil/stamen primordium forming stages, respectively.
Additional file 8: Figure S3 Expression patterns of OsPRR1, OsGI, OsLHY and OsDof under LD and SD conditions in inflorescence. The bar chart shows the expression patterns of OsPRR1, OsGI, OsLHY and OsDof under LD and SD conditions in inflorescence. The y-axis represents the expression ratio relative to the endogenous reference Actin1 gene using qPCR. NLD and SLD represent Nongken 58 and Nongken 58S under LD conditions. NSD and SSD represent Nongken 58 and Nongken 58S under SD conditions. P indicates the pistil/stamen primordium forming stages. Asterisks show statistically significant differences (* = P < 0.05).