The Rice Phytochrome Genes, PHYA and PHYB, Have Synergistic Effects on Anther Development and Pollen Viability

Phytochromes are the main plant photoreceptors regulating multiple developmental processes. However, the regulatory network of phytochrome-mediated plant reproduction has remained largely unexplored. There are three phytochromes in rice, phyA, phyB and phyC. No changes in fertility are observed in the single mutants, whereas the seed-setting rate of the phyA phyB double mutant is significantly reduced. Histological and cytological analyses showed that the reduced fertility of the phyA phyB mutant was due to defects in both anther and pollen development. The four anther lobes in the phyA phyB mutant were developed at different stages with fewer pollen grains, most of which were aborted. At the mature stage, more than one lobe in the double mutant was just consisted of several cell layers. To identify genes involved in phytochrome-mediated anther development, anther transcriptomes of phyA, phyB and phyA phyB mutants were compared to that of wild-type rice respectively. Analysis of 2,241 double-mutant-specific differentially expressed transcripts revealed that the metabolic profiles, especially carbohydrate metabolism, were altered greatly, and heat-shock responses were activated in the double mutant. This study firstly provides valuable insight into the complex regulatory networks underlying phytochrome-mediated anther and pollen development in plants, and offers novel clues for hybrid rice breeding.

and distorted anther lobes ( Fig. 2A-D). In addition, both the scanning electron microscopic and the semi-thin section observations of anther structures from the phyA and phyB single mutants, phyA phyB double mutant and WT showed that the structures of the four anther lobes of the double mutant were different from one another ( Fig. 2E-L), indicating that their developmental synchronization was interrupted. The four lobes remained at different stages, and one or two lobes even could not produce any visible pollen grains (Fig. 2H,L). Before anthesis, the middle layer and tapetum largely degenerated, leaving only the epidermis, which was covered by a cuticular membrane, and the endothecium, which was deposited by U-shaped ligno-cellulosic thickenings adjacent to two splits in the WT and single mutants (Fig. 2M-O). The deposition of secondary thickenings were also observed in the endothecium of some lobes in the anthers of the phyA phyB mutant, but the cell layers of the anther walls were more numerous, at least one lobe was just consisted of several cell layers (Fig. 2H,L and P). The extra cell layers of the anther wall in the phyA phyB mutant were consisted of enlarged cells with low cytoplasmic density, which were distinguished from those of the middle layer and tapetum layer cells in single mutants and WT (Fig. 2P).
Pollen development was severely retarded and pollen viability decreased significantly in the phyA phyB mutant. As shown in Fig. 2H,L, pollen development was affected severely in the phyA phyB mutant, resulting in multiple shriveled pollen grains, most of which were arrested at stage 11. To observe pollen defects in the phyA phyB mutant more clearly, we used transmission electron microscopy to compare pollen grains among WT, single and double mutants. At the mature stage, pollen grains in WT and single mutants were round and enriched with storage materials, such as starch granules and lipid bodies (Fig. 3A-F), which are necessary for pollen viability and function. However, in the phyA phyB mutant, most of the pollen grains were irregular in shape and filled with materials of high or low electronic density, no starch granule were visible ( Fig. 3G-J). In addition, the transmission electron micrographs of pollen wall structures of the single mutants, double mutant and WT showed no significant difference in the exine structures. However, the double mutant formed a much thinner intine than those of WT and the single mutants (Fig. 3H,J). Considering that pectin and cellulose are the major constituents of intine, we hypothesize that phyA phyB mutant pollens contain much less pectin and cellulose.
To examine the vitality of pollen grains in the phyA phyB mutant, Alexander's staining was used to distinguish the viable and aborted pollen grains. As shown in Fig. 4A-C, almost all of the pollen grains in the WT, and phyA and phyB single mutants were stained purple, indicating that they were active. In contrast, approximately 80% pollen grains in the phyA phyB mutant were stained green, which suggested that a large number of pollen grains were aborted (Fig. 4D).
Thus, the reduced fertility of the phyA phyB mutant was determined by two factors, one is the abnormal anther development and failure of anther dehiscence, the other one is reduced pollen vitality.
Overview of anther transcription profiles of phyA and phyB single and double mutants and WT rice. To identify genes responsible for abnormal anther development in the phyA phyB mutant, 12 cDNA libraries used for digital gene expression analysis were constructed using anthers of phyA and phyB single and double mutants and WT at stage 12. Three biological replicates were established for each material. Over 130 million raw reads were generated from the 12 libraries, with phyA mutant libraries producing the fewest raw reads and phyA phyB mutant libraries producing the most (Table 1). After removing low quality reads, as well as reads with adaptors or more than one N base, the number of clean reads varied greatly among different libraries, from 8,326,096 in the phyA library to 14,102,775 in the WT library ( Table 1). The clean reads were then mapped to the rice genome, with no more than one mismatch allowed in the alignment. In total, over 36,000 genes were found to be expressed in rice anther ( Table 1). The number of genes expressed in phyA and phyB single and double mutant anthers was similar to that expressed in WT anthers.

Identification of candidate genes responsible for abnormal anther and pollen development in
the phyA phyB double mutant. Compared with WT anthers, significantly more genes were differentially expressed in the phyA phyB double mutant than in either single mutant. In total, 461 (374 genes), 530 (423 genes), and 2,539 (1,951 genes) differentially expressed transcripts were identified in phyA, phyB and phyA phyB mutants, respectively (Fig. 5A, Table S1). These results indicate that gene expression profiles changes dramatically in the double mutant. Among the 2,539 differentially expressed transcripts between phyA phyB double mutant and WT libraries, a small percentage of differentially expressed transcripts (298, 11.74%) were common to phyA/WT or phyB/WT datasets. The remaining 2,241 differentially expressed transcripts (1,717 genes) were specific to the phyA phyB/WT dataset (Fig. 5A). Because a single mutation of PHYA or PHYB gene has no impact on the fertility, the 2,241 differentially expressed transcripts specific to the phyA phyB/WT dataset were considered to be the candidate genes responsible for the abnormal anther and pollen development of the phyA phyB double mutant. These phyA phyB/WT-specific differentially expressed transcripts were used in subsequent analyses.
Among the candidate transcripts responsible for fertility control in the phyA phyB mutant, 1,111 transcripts (847 genes) were up-regulated and 1,130 transcripts (870 genes) were down-regulated (Table S2). Gene Ontology (GO) analysis revealed that genes involved in stimulus responses, carbohydrate metabolism and transport were over-represented in the up-regulated transcripts, while genes involved in photosynthesis were enriched in the down-regulated ones (Fig. S1). Half of genes involved in stimulus responses in the up-regulated dataset encode  heat-shock factors (Hsfs) and proteins (Hsps) ( Table S2). Among genes involved in carbohydrate metabolism in the up-regulated dataset, genes participate in cellulose degradation take the largest proportion, followed by genes involved in starch metabolism, chitin hydrolysis, galactose catabolism and callose hydrolysis (Table S2).
To confirm gene expression patterns detected by RNA sequencing, 25 genes differentially expressed between the phyA phyB mutant and the WT were examined by quantitative reverse transcription-PCR (qRT-PCR). As shown in Fig. 5B, most gene expression patterns detected by qRT-PCR were similar to those revealed by RNA-Seq. Correlation coefficient between the two methods reached to 0.933. The qRT-PCR results indicate that the RNA-Seq results were quite reliable.
The metabolic profile changed significantly in the phyA phyB mutant. GO analysis of the differentially expressed transcripts that were specific to the phyA phyB double mutant revealed that genes involved in metabolic processes, including carbohydrate metabolism and photosynthesis, were over-represented in the phyA phyB mutant (Fig. S1). To overview the metabolic profile changes in the phyA phyB mutant, we carried out the MapMan analysis tool. As shown in Fig. 6, both major and minor CHO metabolisms were significantly affected in the phyA phyB mutant. Starch and sucrose metabolism-and trehalose biosynthesis-related transcripts were mainly down-regulated, which indicated that these sugars changed significantly in the phyA phyB mutant. Based on the list of genes up-regulated in the phyA phyB mutant, it was found that a large proportion of these genes (35/68, 51.5%) were encoding glycoside hydrolases and β-galactosidases. Members of eight glycoside hydrolase families, 3, 9, 10, 16, 17, 18, 31 and 38, were up-regulated in the phyA phyB double mutant. The up-regulated β-galactosidases belonged to families 11, 17 and 44. Considering that enzymes of glycoside hydrolase families 3, 9, 10 and β-galactosidases can hydrolyze cellulose, cellulose hydrolysis should be accelerated in the phyA phyB mutant. Two genes encoding cellulase, which mainly contributes to cellulose decomposition in plants, were also up-regulated. This result was in accordance with the abnormal anther walls and the thinner intines of phyA phyB mutant pollen grains. In addition, genes involved in two steps of photosynthesis, light reactions and the calvin cycle, were down-regulated significantly (Fig. 6), indicating the phyA phyB mutant produced much less organic materials such as sugars, lipids and proteins, than WT. Most of the differentially expressed genes that participated in fatty acid synthesis, elongation and desaturation, and phospholipid and glycolipid synthesis were down-regulated (Fig. 6). In accordance with this result, the phyA phyB mutant accumulated lower levels of storage materials in the mature pollen grains (Fig. 3H). In contrast, most of the differentially expressed transcripts related to nucleotide metabolism, amino acid metabolism, wax, terpene and flavonoid metabolisms, were significantly up-regulated (Fig. 6). The expression patterns of cell wall-related transcripts were also changed significantly in the phyA phyB mutant. Cell wall modification enzymes, including expansin, xyloglucan endotransglycosylase and endoxyloglucan transferase, were up-regulated, while most of the genes encoding pectinesterase were down-regulated (Fig. 6). All these four enzymes have cell wall remodeling activities, therefore, changes in the expression levels of genes encoding these enzymes may lead to altered anther and pollen wall structures and components in the phyA phyB mutant.
Hsps and Hsfs may play import roles in the anther and pollen development in the phyA phyB mutant. Hsps are a class of universal molecular chaperons that play important roles in various biological activities, such as protein folding, assembly and translocation, and degradation of truncated and misfolded proteins. Hsps can also protect plants against various abiotic stresses by assisting plants to reestablish cellular homeostasis [25][26][27] . Hsps from four major families have been identified in rice, including 29 small Hsp (sHsp), 26 Hsp70, 9 Hsp90 and 10 Hsp100 proteins 28 .
According to the list of up-regulated genes in the phyA phyB double mutant, Hsps formed the largest proportion of genes involved in abiotic stress (Table S2). As shown in Table 2, 38 transcripts (25 genes) were Hsps related to abiotic stress responses. Four categories of Hsps were identified in the up-regulated gene set, including 8 sHsp, 5 Hsp40, 5 Hsp70, 6 Hsp90 and 1 Hsp100 (Table 2). In contrast, fewer Hsps were down-regulated in the phyA phyB double mutant, except for four genes encoding Hsp40 (Table S2). To verify the Hsp expression patterns in rice anthers, qRT-PCR was adopted. As shown in Fig. 5B, all six Hsps were up-regulated significantly in the phyA phyB double mutant, which was consistent with the RNA-Seq results.
As we know, the expression of Hsps is activated by Hsfs which interact with heat shock elements in the promoter regions of Hsps 29 . Consistent with the up-regulation of Hsps in the phyA phyB mutant, five Hsfs, namely OsHsfA2c, OsHsfA2d, OsHsfB2b and OsHsfB2c, were also up-regulated significantly in the double mutant ( Table 2). In addition, only one Hsf, OsHsfC2b, was down-regulated in the phyA phyB double mutant (Table S2). The qRT-PCR results confirmed the expression of the three Hsfs (Fig. 5B). The up-regulation of both Hsfs and Hsps indicate that they might be the key regulators of anther development in the phyA phyB mutant.

Discussion
Phytochromes are a class of photoreceptors responsible for perceiving red/far-red light in plants. In rice, the mutations of phytochromes not only lead to abnormal photomorphogenesis, but also result in reproductive defects, such as small panicles, infrequent anther dehiscence and low seed-setting 23,24 . Previous studies have demonstrated that failure of fertilization in the phyA phyB phyC mutant was caused by impaired dehiscence of the anther wall 23,24 . This study firstly reported that the phyA phyB double mutant also had defect in male fertility, which is similar to that of triple mutant. To identify genes involved in the phytochrome-mediated reproduction process, a series of cytological observations and RNA-Seq investigations of phyA and phyB single and double mutants, and WT were carried out in this study. Cytological examination of slides of developing anthers and Alexander's staining of pollen grains revealed that phyA and phyB single mutants showed similar anther structures and pollen viability as those of WT (Figs 2, 3 and 4). However, the phyA phyB double mutant exhibited hysteretic anther differentiation, more cell layers in at least one lobe, and fewer pollen grains-most arrested at stage 11 in the relatively normal lobes (Fig. 2H,L). These results suggest the existence of novel signaling pathways regulating anther development through phytochromes.
PIFs are a class of basic helix-loop-helix transcription factors that bind directly to phytochromes in the biologically active Pfr form 30 . They usually localize to the nucleus and bind specifically to downstream genes at the light-regulated cis-element G-box and PBE-box in their promoter regions [31][32][33][34][35] . Thus, identifying PIFs and their downstream genes under the regulation of phytochromes is an efficient way to establish novel phytochrome-mediated signaling pathways. In this study, two PIFs were specifically differentially expressed in the phyA phyB mutant (Table S2), which indicates they function downstream of PHYA and PHYB in rice. The two PIFs showed opposite expression patterns; PIF16 was up-regulated significantly in the phyA phyB mutant, whereas PIF14 were down-regulated (Table S2). To identify the candidate genes targeted by the two PIFs, we analyzed the distribution of putative PIF-binding motifs (including G-box and PBE-box) in the 2 kb upstream regions of the 2,241 differentially expressed transcripts (1,717 genes) specific to the phyA phyB double mutant. The result showed that more than half of the differentially expressed genes had at least one G-box or PBE-box in their promoter regions (Fig. S2). In detail, 566 (32.7%) and 912 (53.1%) differentially expressed genes had at least one G-box and PBE-box in their promoter regions, respectively. Among them, 319 genes had both the two cis-elements in their promoter regions. The large number of differentially expressed genes with G-box and/or PBE-box in their promoter regions might act as direct targets of the two PIFs in the signaling pathways of PHYAand PHYB-mediated rice growth and development.
Environmental stresses, such as heat, drought and cold, have adverse impacts on anther development, especially during meiosis and mitosis [36][37][38][39] . Wheat exposed to elevated temperatures during pollen mitosis 1 and 2 triggers two types of abnormal anther development 40 . In one type, tapetum degradation occurs during meiosis, but the first microspore mitosis could not finish, resulting in immature pollen grains with normal exine but no cytoplasm. In the other type, the first mitotic division of the microspores proceeds to completion; however, few microspores undergo complete second mitosis, with mostly sterile pollen grains as a consequence 41 . Prolonged cold treatment of rice at the young microspore stage leads to small anthers and severe male sterility, which results from the hypertrophy and break down of tapetum cells 42 . In major cereals, such as wheat, maize, and rice, water deficit at the meiotic stage leads to partial or complete male sterility through inhibition of the future development of microspores and pollen grains 37,43,44 . The failure of male gametophyte development under stress treatments was considered to be the result of impaired carbohydrate metabolism [43][44][45] .
Carbohydrate metabolism disturbance severely affects anther development, even lead to male sterility 46,47 . Sugars are the key components of developing anthers, especially in the anther wall fraction where sugar synthesis, secretion, storage, mobilization, and regulation take place 48 . At the early stage of rice anther development, starch is the most abundant sugar in the endothecium; it then decreases after meiosis, but reaches a much higher level during the first pollen grain mitosis 44 . In addition, cellulose and lignified thickenings, which are essential for anther dehiscence, are deposited in the endothecium during microspore maturation 49 . During rice anther maturation, reducing sugar contents increase while those of non-reducing sugars decrease 44 . Water deficit has no obvious effect on starch content at the meiotic stage, but tends to have a significant inhibitory effect during later stages. After first mitosis, starch content in drought-stress-treated anthers is less than half that of control plants 44 . Following meiosis, reducing sugar contents of drought-stress-treated anthers during subsequent stages are lower than controls, whereas non-reducing sugar contents are higher 44 . Comparing the expression profiles between phyA, phyB and phyA phyB mutants and that of WT revealed that the carbohydrate metabolic profiles changed significantly in the double mutant. Major and minor CHO metabolisms were altered greatly, and photosynthesis was inhibited in the phyA phyB mutant (Fig. 6). Examinations of anthers and pollens in the phyA phyB double mutant showed similar results. In the phyA phyB mutant, more than one anther lobe did not generate tapetum which supplies nutrients, such as sugars, lipids and proteins, for pollen development and maturation 47 . Instead, these anther lobes produced several cell layers in anther locules without any visible pollen grains (Fig. 2H,L). In addition, the phyA phyB mutant produced much less pollen grains, most of which were arrested at stage 11 (Fig. 2H,L). Transmission electronic microscopy analyses revealed that the pollen intines were much thinner, and the storage materials, especially starches and lipids, were absent in the phyA phyB mutant (Fig. 3G-J). The cytological observations of the anthers and pollen grains of the double mutant indicated its metabolism changed significantly, with sharp decreases in cellulose, starch and lipid contents.
According to our study, the mutations of PHYA and PHYB seem to have impacts on anther development similar to those of stress treatments. In this study, 28 Hsps and 5 Hsfs were up-regulated in the phyA phyB mutant (Table 2), which indicated that PHYA and PHYB mutations altered rice sensitivities to environmental stress. Consistent with this hypothesis, phyBD, phyABDE and phyABCD showed sequentially enhanced resistance to abscisic acid and NaCl treatments in Arabidopsis 50 . Similarly, the phyA phyB mutant in rice exhibited less sensitivity to NaCl and low temperature stresses than WT and single mutants (Unpublished data). Previous studies revealed that Hsps and Hsfs play important roles in anther and pollen development. Many Hsps are expressed at different stages of anther development in various plant species, such as Arabidopsis, tobacco and tomato [51][52][53][54] . For example, the Hsp40 protein TMS1 is required for normal pollen tube growth under heat shock environments in Arabidopsis 55 . In tobacco, two small Hsps are up-regulated significantly in the transition from unicellular to bicellular pollen, while nine other small Hsps are down-regulated at the same stage in tobacco 54 . HsfA2 and Hsp17-CII are activated at the tomato diploid pollen mother cell stage, and their expression levels are up-regulated significantly under short and prolonged heat shock treatments in subsequent anther developmental stages 56 . Considering both the abiotic stress treatments and the absence of phytochromes could rapidly induce the expression levels of stress-responsive genes, such as HSFs and HSPs, we suggest that the phytochrome status has a great impact on plant stress responses. Previous studies revealed that male sterility induced by abiotic stress treatments were caused by impaired carbohydrate metabolisms in many plant species [43][44][45] . In this study, we also found that the carbohydrate metabolism profiles changed significantly in the phyA phyB mutant. Thus, we inferred that the growth and development of the double mutant wasted much carbohydrates and energy, and that the male sterility of the double mutant resulted from the imbalanced carbohydrate metabolism.

Methods
Plant materials and growth conditions. Four plant materials were used in this study: WT rice (Oryza sativa L. cv. Nipponbare), phyA (phyA-4), phyB (phyB-1) and phyA phyB (phyA-4 phyB-1) mutants 21,22 . The phyA mutant was obtained by screening Tos17 retrotransposon mutant panels 21,22 , which has a Tos17 insertion in the middle of the second exon of the PHYA gene. The phyB mutant was isolated from γ-ray-mutagenized rice seeds 22 , which has 1 bp insertion (C) in the the first exon. All the mutants used in this study were homozygous in the background of Nipponbare. Seeds were surface sterilized in 70% (v/v) ethanol for 30 s and then placed in 5% NaClO (v/v) for 20 min. After rinsing six times in sterile double-distilled water, the seeds were incubated in the dark at 28 °C for 3 days to induce germination. During July to August, seedlings were transplanted into soil at the farm of the Shandong Academy of Agricultural Sciences in Jinan, Shandong, China (latitude 36°40′N; longitude 117°00′E). Anthers from the four plant materials were collected before anthesis. The collected anthers were frozen in liquid nitrogen and transferred to a −80 °C freezer until RNA extraction.
Anther structural analysis and pollen viability test. Anthers from the WT, phyA, phyB and phyA phyB mutants were fixed before anthesis as described by Li et al. 10 . Then they were embedded in Spurr's resin (Sigma-Aldrich, St Louis, MO). Semi-thin sections were generated with an Ultracut E ultramicrotome (Leica Microsystems). After staining with 0.5% (w/v) toluidine blue-O (Sigma-Aldrich, St Louis, MO), the sections were observed and photographed under an Olympus BX51 microscope.
Transmission electron microscopic analysis of anthers from WT, phyA and phyB single, and phyA phyB double mutants were sampled at stage 12 and fixed for 4 h in phosphate buffer, at pH 7.0, supplemented with 2.5% glutaraldehyde. Then they were rinsed with the same buffer and post-fixed for 2 h in phosphate buffer containing 1% osmium tetroxide (Pelco International, Redding, CA, USA) at room temperature. After dehydration by a graded ethanol series (50%, 70%, 80%, 90%, 95% and 100%), the specimens were embedded in Spurr resin (SPI). Ultrathin sections (50 to 70 nm thick) collected on uncoated nickel grids (100 mesh) were stained with 2% uranyl acetate and 1% lead citrate, then examined at 80 kV using a transmission electron microscope (Hitachi 7650, Ibaraki, Japan).
The samples used for scanning electron microscopic analysis were the same as those of used in transmission electron microscopic analysis. All the samples were embedded for 4 h in 0.1 M sodium phosphate buffer solution containing 2.5% glutaraldehyde, at pH 7.0, rinsed three times in 0.1 M phosphate buffer, at pH 7.0, and further fixed in 0.1 M sodium phosphate buffer solution containing 1% osmium tetroxide for 2 h. They were rinsed again in the same buffer, dehydrated in an acetone series from 50% to 100% and then rinsed three times with 100% acetone, which was exchanged three times with isoamyl acetate (Sigma-Aldrich). They were then processed for critical-point drying using liquid CO 2 , and gold-coated (10-nm thickness) in an E-1045 sputter coater (Hitachi, Tokyo, Japan). The specimens were then examined in a TM3000 (Hitachi, Tokyo, Japan) with an accelerating voltage of 15 kV.
For the pollen viability test, anthers were stained in Alexander stain for 1~2 days according to Alexander (1969) 57 .
RNA extraction, cDNA library construction, and sequencing. Total RNA was extracted from rice anthers using the modified CTAB method as described in Xu et al. 58 . The quantity and purity of total RNA were analyzed on a Bioanalyzer 2100 instrument (Agilent, Santa Clara, CA, USA). RNAs with integrity numbers above 7.0 were used for library construction. Approximately 10 μg total RNA was used for generating cDNA according to the Illumina protocol (Illumina, Aan Diego, CA, USA). The products were single-end sequenced on an Illumina Hiseq. 2500 sequencer at LC Sciences (Hangzhou, China) following the vendor's recommended protocol. Processing of fluorescent images into sequences, base-calling and quality value calculations were carried out using the Illumina data processing pipeline (version 1.8).
Read alignment and assembly. Raw reads containing adaptor sequences, with low quality and unknown nucleotides N were filtered to obtain clean reads. The clean reads were then subjected to quality assessment, which included classification of total and distinct reads and percentage counts in libraries, saturation analysis of libraries and correlation of biological replicates. All clean reads were mapped to the rice genome (http://phytozome.jgi. doe.gov/pz/portal.html) using Bowtie (version 1.0.0), with only one base-pair mismatch allowed. The number of perfect clean reads corresponding to each gene was calculated and normalized to the number of RPKMs 25 .
Differentially expressed genes identification and functional analysis. Transcript profiles of all the sequencing data were analyzed by calculating the number of reads per kilobase of exon model per million mapped reads (RPKM) 59 . To identify candidate genes involved in anther development in the phyA phyB mutant, a series of comparisons between phyA and phyB single and double mutants, and WT anther expression profiles were carried out using DEseq. 60 . Transcripts with fold change ≥2, FDR ≤ 0.01 and read number ≥20 in at least one of the two comparison samples were considered to be differentially expressed. Gene Ontology (GO) analyses were performed using AgriGO 61 .
qRT-PCR analysis. Three biological replicates of total RNA were used for qRT-PCR analysis. Total RNA was extracted as described above. After DNase I treatment, approximately 4 μg of total RNA was reverse transcribed using M-MLV RTase cDNA syntheis kit (Takara, Dalian, China).
qRT-PCR was performed on a StepOne plus Real-Time PCR system (ABI, Carlsbad, CA, USA) using SYBR Green PCR master mix (Takara, Dalian, China). The EF-1α gene was used as an internal control. Primers used in qRT-PCR experiments were listed in Table S3.