Transcriptome Analysis Reveals Photoperiod-Associated Genes Expressed in Rice Anthers

Environmental conditions, such as photoperiod and temperature, can affect male fertility in plants. While this feature is heavily exploited in rice to generate male-sterile lines for hybrid breeding, the underlying molecular mechanisms remain largely unknown. In this study, we use a transcriptomics approach to identify key genes and regulatory networks affecting pollen maturation in rice anthers in response to different day lengths. A total of 11,726 differentially expressed genes (DEGs) were revealed, of which 177 were differentially expressed at six time points over a 24-h period. GO enrichment analysis revealed that genes at all time points were enriched in transport, carbohydrate, and lipid metabolic processes, and signaling pathways, particularly phytohormone signaling. In addition, co-expression network analysis revealed four modules strongly correlated with photoperiod. Within these four modules, 496 hub genes were identified with a high degree of connectivity to other photoperiod-sensitive DEGs, including two previously reported photoperiod- and temperature-sensitive genes affecting male fertility, Carbon Starved Anther and UDP-glucose pyrophosphorylase, respectively. This work provides a new understanding on photoperiod-sensitive pollen development in rice, and our gene expression data will provide a new, comprehensive resource to identify new environmentally sensitive genes regulating male fertility for use in crop improvement.


INTRODUCTION
Hybrid breeding has made a great contribution to the yield increase of rice (Oryza sativa) relative to inbred varieties, with improvements of up to 20% (Khush, 2001;Michel et al., 2001;Yuan, 2004;Cheng et al., 2007). Originally, after the discovery of a wild abortive cytoplasmic male sterile (CMS) line in the 1970s, a three-line system containing the male-sterile line, a maintainer line, and a restorer line was used (Yuan, 1984). Due to the time-consuming and germplasm-limited nature of using separate maintainer and restorer lines (Chen et al., 2010), a more modern two-line system was developed in 1973, based on the discovery and application of environmentally sensitive genic male sterile lines (EGMS) (Shi, 1981). Photoperiod-sensitive (PGMS) lines and thermosensitive (TGMS) lines have been successfully developed and applied widely in rice breeding (Mou, 2016).
Since photoperiod is more stable and predictable than temperature, breeders have focused on PGMS lines. Several genes, such as photoperiod-sensitive male sterility-1 (pms1) (Zhou et al., 2011;Fan et al., 2016), pms2 (Zhang et al., 1994), pms3 (Ding et al., 2012a,b;Zhou et al., 2012), Programmed Cell Death 5 (OsPDCD5) , and Carbon Starved Anther (CSA) , have been reported to control PGMS. Pms1 encodes a long non-coding RNA (lncRNA) that acts as a phased small-interfering RNA (phasiRNA)-generating locus to produce the PMS1T transcript (21-nt phasiRNA). PMS1T is highly expressed during pollen formation in wildtype plants under long-day (LD) conditions, and a mutation in PMS1T that causes an accumulation of phasiRNAs causes PGMS (Fan et al., 2016). A single-nucleotide polymorphism (SNP) in pms3, another lncRNA required for normal pollen development, reduces pms3 transcript levels under LD conditions, leading to PGMS in a japonica cultivar, NK58S (Ding et al., 2012a,b). Subsequently, Zhou et al. (2012) identified the same SNP in NK58S-derived cultivars and reported the role of this mutation in producing one mutated small RNA, leading to PGMS and TGMS in japonica and indica backgrounds, respectively (Zhou et al., 2012). OsPDCD5 is an ortholog of mammalian programmed cell death 5, whose decreased expression induces photoperiodsensitive male sterility in rice under LD photoperiods (≥13.5 h sunlight) . CSA encodes an R2R3 MYB transcription factor that regulates sugar accumulation in rice anther sink tissue; its knockout mutant accumulates insufficient sugar and starch levels in anthers under short-day (SD) conditions, resulting in male sterility . A TGMS line currently used for hybrid seed production, thermosensitive genic male sterile 5 (tms5), has been shown to contain a mutation in a conserved ribonuclease, RNase Z (Zhou et al., 2014). These findings reveal the complexity of the interaction between genetic components and environmental signals in determining male fertility in plants (Kim and Zhang, 2018). However, the global understanding of male reproduction in plants in response to different photoperiods remains largely unknown.
To comprehensively analyze molecular mechanisms controlling plant male fertility in response to photoperiod, we performed a global RNA sequencing (RNA-seq) analysis on anthers grown in LD and SD conditions. Genes associated with carbohydrate transport and metabolism, and phytohormone signaling, were highly correlated with day length during late pollen developmental. This work provides insights into the molecular mechanisms of photoperiod-dependent male fertility, which can lead to identification of new PGMS genes for hybrid breeding.

Plant Materials
Rice variety O. sativa L. ssp. japonica 9522 was grown in paddy fields around Shanghai (China) during the normal rice growing season (June-August). Photoperiod treatment was started after panicle initiation and continued until all samples were collected: LD was the natural daylight condition, with a photoperiod of ∼14 h; SD conditions were simulated by covering plants with black cloth for the last 2 h of natural daylight to create a 12h photoperiod. Duplicate samples of anthers were collected at mitosis I (stage 11) based on anther length (about 2.2 mm), which was described by Zhang et al. (2011) at six time points (00:00, 04:00, 08:00, 12:00, 16:00, and 20:00; Supplementary Figure 1). All samples were frozen immediately in liquid nitrogen and stored at −80 • C.

RNA Sequencing
RNA isolation, library construction, and sequencing were performed by Novelbio Ltd (Shanghai, China). Briefly, TRIzol R reagents (Invitrogen, Carlsbad, CA, United States) were used to isolate RNA from each of two biological samples from each experimental condition. RNA quality was assessed by 1% agarose gel electrophoresis, Nanodrop (ThermoFisher), and Bioanalyzer (Agilent 2100), with A260/A280 = 1.8-2.2, and RNA integrity number (RIN) > 6.5. Libraries were synthesized using the standard Ion proton RNA-Seq Kit v2.0, according to the manufacturer's instructions. Library sequencing was conducted by Ion Proton TM (Life Technologies, Carlsbad, CA, United States).

Sequence Data Processing, Mapping, and Annotation
The genome and annotation of the rice reference cultivar Nipponbare was downloaded from the MSU Rice Genome Annotation Project 1 (Kawahara et al., 2013). FastQC (Version 0.11.1) was used to remove reads with 20% of the base quality lower than 13. Clean reads were mapped to the reference genome using Mapsplice v2.1.8 and then read counts were calculated (Wang K. et al., 2010), which were normalized to reads per kilobase per million (RPKM) to obtain relative levels of transcript expression. Counts and RPKM for each gene are given in Supplementary Data Sheet 1. A summary of sequencing data for each sample is given in Supplementary  Table 1 and Supplementary Figure 2. The sequencing data were also uploaded in NCBI (GSE163030). Expression data for genes of interest in other vegetative and reproductive tissues (Supplementary Figure 5B) were sourced from the Rice Expression Profile Database (RiceXPro) (Sato et al., 2012).
Quantitative Reverse Transcription (qRT)-PCR Validation cDNA was synthesized using the Fast RT Kit (Tiangen Biotech, Beijing), following the manufacturer's instructions. qRT-PCR analysis was performed on the LightCycler 96 (Roche). The OsActin gene was used as internal control, and the relative expression level of target genes was calculated based on the 2 − Ct method (Livak and Schmittgen, 2001). Primer sequences for qRT-PCR were found in qPrimerDB (Supplementary Data Sheet 9) (Lu et al., 2017).

Construction of Gene Co-Expression Networks
Weighted gene co-expression network analysis R package (version 3.6) (Langfelder and Horvath, 2008) was used to construct co-expression networks of all expressed genes with automatic, one-step network construction and module detection. First, the soft thresholding power β was chosen based on the lowest power for which the scale-free topology fit index reached a high value. The function "blockwiseModules" was then used to construct topological overlap matrix (TOM) and module detection. The associations between modules and trait were estimated using Pearson correlation coefficient between module eigengenes and phenotype where 0/1 were used to define SD/LD or light/dark conditions. The expression patterns of four traitrelated modules were plotted by "ggplot2" (Wickham, 2016). Hub genes were identified by both gene significance and module membership. The module membership (MM) is the correlation of gene expression profile with module eigengenes, and gene significance (GS) is the association of individual genes with that trait. Hub genes were set to a threshold of GS and MM > 0.8 (Horvath and Dong, 2008).

Visualization of Co-expression Genes
Cytoscape (V3.7.2) (Shannon et al., 2003) was used to visualize co-expression networks. Input data from turquoise module were chosen based on the weight in terms of correlation value from the TOM (Topological Overlap Matrix Value), so that a higher value refers to a strong co-expression. The top 100 highest correlations with CSA and Ugp1 genes were chosen and highlighted the top 1% correlation value of genes.

DEGs in SD and LD Anthers
A total of 11,726 unique DEGs were identified between SD and LD conditions (Supplementary Data Sheet 2). Six genes 2 http://systemsbiology.cau.edu.cn/agriGOv2/ from the RNA-seq data were selected for confirmation by quantitative reverse-transcription PCR (qRT-PCR), and patterns of expression with the two methods were found to be consistent (Supplementary Figure 3).
The highest numbers of DEGs were expressed at 12:00, with 4552 genes up-regulated and 4662 genes down-regulated in SD anthers compared with LD anthers (Figure 1A and . Gene numbers in each cluster given in brackets. DEGs assembled to seven clusters: three clusters up-regulated in SD (2, 3, and 6); three clusters up-regulated in LD (4, 5, and 7); and Cluster 1, with higher daytime expression under SD conditions, and nighttime expression under LD conditions. The value of center point of seven clusters represents relative expression. (C) Venn diagram showing all DEGs at six time points. A total of 177 genes showed differential expression between SD and LD at all time points. Figure 4). Interestingly, the lowest number of DEGs was observed 4 h later at 16:00, with only 824 up-regulated and 546 down-regulated genes in SD anthers. All 11,726 DEGs assembled into just seven clusters, confirming that the largest differences in expression between LD and SD appeared at 12:00, followed by 08:00 and 04:00 ( Figure 1B). Additionally, the 12:00 time point saw expression of the largest number of unique DEGs (4502), whereas the least numbers of unique DEGs occurred at 20:00 (55) and 00:00 (64) (Figure 1C). These results suggest that photoperiod dramatically affects gene expression patterns in anthers, with the largest changes occurring in the morning (04:00, 08:00, and 12:00), even before dawn.

Anther DEGs at All Time Points
A total of 177 genes were differentially expressed at all time points in SD and LD anthers ( Figure 1C). Most of these genes (138) were up-regulated in LD anthers, with low or undetectable expression in SD anthers (Supplementary Figure 5A). A large proportion of these DEGs (51 genes, 29%) exhibit significantly higher expression in the anther than in other rice tissues (Supplementary Figure 5B), including previously reported anther-related genes such as ATP Binding Cassette G15 (OsABCG15) (Wu et al., 2014), LESS ADHESIVE POLLEN 6 (OsLAP6) (Zou et al., 2017), and HOTHEAD-Like1 (HTH1) (Liu et al., 2016). These 177 DEGs were found in GO categories primarily associated with cell wall synthesis and metabolism, including hydrolytic enzyme expression and regulation, indicating that the anther cell wall structure and composition may be photoperiod-sensitive (Supplementary Figure 5C). These genes may play key roles in anther development and pollen formation, and can be used as markers to distinguish SD and LD conditions.

Functional Analysis of DEGs in SD and LD Anthers
Gene ontology and kyoto encyclopedia of genes and genomes functional analysis of anther DEGs indicated that genes involved in transport, lipid and carbohydrate metabolism, and metabolic regulation were most affected by photoperiod (Figure 2 and Supplementary Figure 6). GO enrichment analysis revealed that DEGs at all time points were enriched in GO "biological process" categories of "transport, " "localization, " "establishment of localization" (transport is a sub-item of these two terms), "carbohydrate metabolic process, " and "lipid metabolic process" (Figure 2A and Supplementary Data Sheet 3). DEGs involved in regulatory processes were generally enriched at 12:00 and 04:00, while metabolic processes were enriched at 08:00 and 12:00 (Figure 2A). In "molecular function" categories, the GO terms of "transport activity" and "enzyme regular activity" were enriched at most time points except 16:00 and 04:00, respectively, while "catalytic activity" was enriched at 08:00 and 00:00 ( Figure 2B). For the "cell component" GO category, DEGs associated with "membrane" and "cell wall" were seen at all time points, while specific cellular components were most affected at 12:00 and 16:00 ( Figure 2C).
KEGG analysis showed that genes encoding proteins involved in "pentose and glucuronate interconversions" in carbohydrate metabolism were enriched at 08:00, 12:00, 20:00, and 00:00 (Supplementary Data Sheet 4 and Supplementary Figure 6). UDP-glucose pyrophosphorylase 2 (OsUgp2), a reported TGMS gene required for starch accumulation in pollen (Mu et al., 2009), was present in this pathway. In addition, genes associated with "starch and sucrose metabolism" and "carbon metabolism" were specifically enriched at 12:00, indicating a major difference in carbohydrate metabolism between SD and LD anthers at this time. Genes encoding "photosynthesis-antenna proteins" were enriched at all time points except 08:00. Rice chlorophyll a/b binding protein gene, OsCAB1R, a circadian clock-controlled gene (Sugiyama et al., 2001), is part of this pathway, indicating that rhythm genes may play photoperiod-sensitive roles in the anther. Calmodulin-related genes were mainly enriched in "plant-pathogen interaction" at 08:00, 20:00, and 04:00, such as calcium-dependent protein kinase 9 (OsCPK9), OsCPK21, and OsCPK25/26. OsCPK9 overexpression has been shown to improve pollen fertility under drought stress and shows sensitivity to ABA (Wei et al., 2014). Specific enrichment of mitogen-activated protein kinase (MAPK) pathway genes involved in defense against biotic stress was observed at 04:00.

DEGs Involved in Transport and Carbohydrate and Lipid Metabolism
Genes involved in transport were differentially expressed in anthers at all time points under the two photoperiod conditions, including 189 genes encoding proteins involved in ion transport, 82 genes associated with protein transport, and 54 genes associated with lipid transport ( Figure 3A and Supplementary Data Sheet 5). These genes could be grouped in two main clusters that were generally up-regulated in either LD or SD conditions. Genes up-regulated in SD anthers were generally most highly expressed at 12:00, but genes up-regulated in LD anthers had high expression at the 08:00 time point as well. For ion transport, DEGs predicted to encode potassium transporters, including high-affinity K + transporter (HKT) (Jabnoune et al., 2009) and HAK family (Bañuelos et al., 2002) proteins, as well as sodium, calcium, hydrogen (proton), iron, boron, and zinc transporters were represented at all time points. For protein transport, genes encoding 15 rat sarcoma-related transporter and 17 mitochondrial translocase proteins were differentially expressed. A small GTPase protein, OsRacD, which has been shown to participate in light signal transduction and the fertility switch of NK58S (Ye et al., 2004), was more highly expressed under SD conditions at 08:00 and 12:00. For lipid transport, DEGs encoding 44 lipid transfer protein (LTP) family members were detected. Among these was Photoperiod-sensitive dwarf 1 (Psd1), a non-specific LTP gene whose mutation can cause dwarfism under LD and low-temperature conditions (Deng W. et al., 2020). Psd1 was up-regulated under SD conditions at 12:00. Other transport-associated DEGs included sugar transporters, such as monosaccharide transporters (MSTs) and sucrose transporters (SUTs). For example, MST7 was up-regulated under SD conditions at all time points, while MST5 and SUT1 were up-regulated under LD conditions at all times except 16:00.
In the GO category "carbohydrate metabolism, " DEGs encoding glycosyl hydrolase and transferase proteins were highly  represented, including 250 hydrolases and 118 transferases ( Figure 3B and Supplementary Data Sheet 5). These genes again clustered into two groups, each up-regulated under one of the two photoperiod conditions, and genes up-regulated under SD conditions were again most highly expressed at 12:00 ( Figure 3B). DEGs up-regulated under LD conditions exhibited two patterns of expression: a larger group most highly expressed at 08:00 and 12:00 pm, while a second, smaller group exhibited peak expression at 00:00 and 04:00 ( Figure 3B). Hexokinase (HXK) family genes including HXK5 (Cho et al., 2009), HXK7 (Kim et al., 2016), and HXK10 (Xu et al., 2008) as well as four invertases (INV) (Oliver et al., 2005;Deng X. et al., 2020) were identified, and these were spread across all three major patterns of expression. DEGs encoding proteins involved in lipid metabolism related genes also exhibited three distinct patterns of expression ( Figure 3C and Supplementary Data Sheet 5). Genes generally up-regulated in SD anthers had maximum expression at 12:00; genes up-regulated in LD anthers generally peaked at 08:00 and 12:00, while a third group of genes exhibited maximum expression in SD anthers at 12:00 and in LD anthers at 04:00. Two genes associated with male sterility, Wax-deficient anther 1 (Wda1) and HUMIDITY-SENSITIVE GENIC MALE STERILITY 1 (HMS1) Chen et al., 2020), were among the genes up-regulated in SD anthers. Wda1 is associated with biosynthesis of very-long-chain fatty acids for the formation of anther epicuticular wax crystals, whose mutation causes complete male sterility . HMS1 encodes a β-ketoacyl-CoA synthase, and its mutant is humidity-sensitive genic male sterile .

DEGs Involved in Circadian Rhythm
Based on our finding that OsCAB1R, a circadian clock-controlled gene, was a DEG in the "photosynthesis-antenna proteins" GO class (Supplementary Data Sheet 4), we analyzed, and found, differential expression of other genes known to be affected by circadian rhythm (Supplementary Figure 7 and Supplementary Data Sheet 5). Rice CIRCADIAN CLOCK ASSOCIATED 1 (OsCCA1), homologous to Arabidopsis AtCCA1 (Izawa et al., 2002), showed a higher expression in the early morning, with almost twofold higher expression under LD compared with SD conditions at 04:00; ZEITLUPE 2 (ZTL2) had a similar pattern but peaked under LD at 08:00. At 12:00, OsCCA1 was expressed nearly seven times higher in SD compared with LD anthers. Rice PSEUDO-RESPONSE REGULATOR (OsPRR), OsPRR95, homologous to Arabidopsis circadian oscillator AtPRR5 and AtPRR9 (Murakami et al., 2005), showed fivefold expression higher at 12:00 under SD compared with that of LD conditions. OsPRR37 shows up-regulation at night under LD (20:00, 00:00, and 04:00), four times higher at 00:00 and 04:00, while OsPRR73, Early flowering 3 (OsELF3), and ZTL1 showed little difference between SD and LD conditions. Rice LUX ARRHYTHMO (OsLUX), homologous to a nighttime repressor of circadian gene expression AtLUX (Murakami et al., 2007), was also upregulated at 16:00 in SD conditions. Other rice circadian clock or flowering genes FLAVIN BINDING, KELCH REPEAT, F-BOX1 (OsFKF1), TIMING OF CAB EXPRESSION 1 (OsTOC1), and Gigantea (OsGI) (Hayama et al., 2003;Murakami et al., 2003) were expressed more highly at 16:00 in LD anthers than in SD anthers, although by less than a twofold difference. These results indicated that rhythm genes in the anther are sensitive to photoperiod, particularly under SD conditions.

Construction of Co-expression Gene Networks
To examine the regulatory networks of genes under SD and LD conditions, 39,966 genes (99.8%) were categorized into 34 modules based on scale-free topology model (β = 7; Supplementary Figure 8). The remaining 78 genes that did not assemble into existing modules were grouped into a 35th "module" (gray). Standard modules ranged in size from 177 to 5852 genes, e.g., black contained 1496 genes that generally expressed highly at LD 12:00; dark red had 555 genes that generally expressed highly at SD 16:00 (Supplementary Table 2 and Supplementary Figure 9B).

Identification of LD-and SD-Related Modules
Variation in gene expression was higher between SD and LD conditions than between dark/light conditions (Figure 4), suggesting that photoperiod has a significant effect on anther development. Expression of 36% of genes correlated with SD or LD conditions, suggesting that a large proportion of genes responded directly or indirectly to day length. Gene expression in four modules in particular was observed to correlate significantly (p < 0.05) with photoperiod: genes in turquoise (5852 genes, r = 0.76) and magenta (1132 genes, r = 0.71) modules were upregulated in LD, while genes in yellow (4929 genes, r = −0.71), and blue (2,698 genes, r = −0.71) modules were up-regulated in SD (Supplementary Table 2, Supplementary Figure 9, and Supplementary Data Sheet 6).

Functional Specific Enrichment Analysis of Related Modules
Gene ontology analysis revealed enrichment of different functional categories of genes among modules (Figure 5). The turquoise, magenta, and blue modules were enriched for biological processes including "transport, " "localization, " "establishment of localization, " and "carbohydrate metabolic process, " and contained genes encoding various potassium (HKT and HAK family) and sugar (MST, SUT) transporters, and carbohydrate metabolic proteins such as sucrose synthase (SUS) and invertases (Supplementary Data Sheet 7).
Genes associated with "lipid metabolic process" were enriched in turquoise and blue modules (Figure 5A), which is also consistent with the GO enrichment of DEGs (Figure 2A). The blue module was enriched in more GO categories than other modules, including "protein metabolic process" and "macromolecule modification" (Figure 5A). In the "protein metabolic process" category were TGMS genes, such as Temperature-sensitive Male Sterile 10 (TMS10) and TMS10-Like gene (TMS10L) (Yu et al., 2017) and ubiquitin fusion ribosomal protein L40 (UbL404), which acts downstream of TMS5 (Zhou et al., 2014), as well as genes encoding light-responsive proteins, such as photoreceptors Phytochrome B(PHYB) and PHYC, and FIGURE 4 | Module-trait associations. Correlation coefficients between different modules and traits show in matrix. Each cell contains correlation coefficient with p value in bracket. Photoperiod, short-day, and long-day condition; LD acts as positive control, in which red color shows positive correlation with LD; SD is the opposite; Light and Dark: daytime vs. night conditions, daytime is defined as positive blank including samples at 08:00, 12:00, and 16:00; night includes 20:00, 00:00, and 04:00.
KEGG enrichment analysis revealed genes associated with "oxidative phosphorylation, " the "TCA (tricarboxylic acid) cycle, " and "carbon metabolism pathway" in the turquoise module (Supplementary Figure 10). The magenta module was enriched in diterpenoid biosynthesis, while yellow module was enriched in categories associated with DNA replication, RNA transport, and ribosome biosynthesis. These results indicate that protein metabolism and transcription/translation processes are also affected by day length.

Identification of Hub Genes
Hub genes are those with a high degree of connectivity and co-relationships with genes in modules of interest for a given phenotype (Langfelder and Horvath, 2008). Based on the module membership and gene significance, 239, 62, 139, and 56 hub genes were identified in turquoise, magenta, blue, and yellow modules, respectively (Supplementary Data Sheet 6). As identified in our earlier GO and KEGG analyses, these genes were involved a broad range of processes, including macromolecule, carbohydrate, and lipid metabolisim and transport. Of particular interest, six hub genes have previously been reported to be involved in male sterlity or low fertility, including LOC_Os01g16810 (CSA) , LOC_Os09g38030 (UDPglucose pyrophosphorylase gene 1, OsUgp1) (Chen et al., 2007), LOC_Os01g54620 (cellulose synthase catalytic subunit genes 4, OsCesA4) (Tanaka et al., 2003), LOC_Os02g04840 (collapsed abnormal pollen 1, CAP1) (Ueda et al., 2013), LOC_Os04g39980 (dioxygenase for auxin oxidation gene, DAO) (Zhao et al., 2013), and LOC_Os06g04090 (secondary wall NAC 1, OsSWN1) (Chai et al., 2015). The discovery of known key regulators of pollen sterility indicates the robustness of our co-expression and hub gene construction method to mine for key anther regulator genes. CSA encodes a key R2R3 MYB transcription factor that controls male fertility in a photoperiod-dependent manner . OsUgp1 is required for callose depostion in pollen mother cells in meiosis; knockdown of OsUgp1 causes a thermosensitive male sterile defect (Chen et al., 2007). OsCesA4 is required for cellulose synthesis, and its mutant shows brittle culm and low fertility (Tanaka et al., 2003). CAP1 encodes an arabinokinase-like protein required for nucleus, cytoplasm, and intine cell wall development in rice pollen. DAO encodes a 2-oxoglutarate-dependent-Fe (II) dioxygenase essential for anther formation and dehiscence, pollen fertility, and seed initiation (Zhao et al., 2013). OsSWN1 is a NAC transcription factor gene involved in secondary cell wall biosynthesis, whose down-regulation causes a decrease in lignin and an increase in polysaccharide content in cell walls, resulting in anther dehiscence defects (Chai et al., 2015). The identification of these genes in our analyses suggests that these previously characterized pollen development regulators could also have a day lengthdependent role.
An additional 30 hub genes encoding transcription factors were detected, including 5 previously unreported MYB proteins and 10 zinc-finger proteins ( Table 1). The function of these unreported genes in regulating anther development needs to be further characterized, but it seems likely that they will be key regulators of rice anther development particularly affected by photoperiod.

Genes Highly Co-expressed With CSA and Ugp1
Two hub genes, CSA and Ugp1, have been reported to be involved in rice photoperiod-sensitive and thermosensitive male fertility, respectively. We visualized genes with high correlation to CSA and Ugp1 (Figure 6 and Supplementary Data Sheet 8). Notably, CSA and Ugp1 shared 74 overlapped correlation genes, such as MATRILINEAL (OsMATL) and OsHXK5 (Cho et al., 2009;Yao et al., 2018), which were up-regulated in LD anthers, unlike that of CSA and Ugp1-as undirected network relationships. Four co-expressed genes were also selected for qRT-PCR analysis and confirmed this prediction outcome (Supplementary Figure 11). OsMATL encodes a pollen-specific phospholipase that is involved in haploid induction, whose knockout mutant displays reduced seed set (Yao et al., 2018). OsHXK5 encodes a hexokinase that functions as a glucose sensor. Plants overexpressing OsHXK5 showed a hypersensitive and slow growth phenotype (Cho et al., 2009). Six TFs such as MYB, Zinc Finger, and MADS-Box FIGURE 6 | High co-expressed genes of CSA and Ugp1. Blue indicates the top 1% TOM value genes; rectangle represents hub genes. OsGA20ox3 (gibberellin 20 oxidase 3, LOC_Os07g07420); OsMATL (MATRILINEAL, LOC_Os03g27610); CSA (Carbon Starved Anther, LOC_Os01g16810); Ugp1 (UDP-glucose pyrophosphorylase, LOC_Os09g38030); Invertase (LOC_Os12g37480, LOC_Os02g01310); MYB-TFs (LOC_Os01g51260, LOC_Os06g46560); MADS-box genes (LOC_Os06g11970); Zinc Finger genes (LOC_Os02g44120, LOC_Os01g68900, and LOC_Os12g08070); F-box (LOC_Os04g19800); Box 1, lipid metabolic process-related genes; Box 2, carbohydrate metabolic process-related genes; Box 3, Transport-related genes. proteins were in network; 5 lipid metabolic-associated genes, 11 carbohydrate metabolic genes, and 11 transport genes were also detected. In addition, there were 38 hub genes highly coexpressing with CSA and Ugp1. These results indicate that CSA and Ugp1 are likely to share regulatory pathways in controlling environmentally mediated male development in rice.

DISCUSSION
Transcript profiling of anthers grown under different photoperiod conditions provides information to help understand the molecular mechanisms that affect male fertility in response to day length and thus improve our ability to create new PGMS lines for rice breeding. In this study, we have conducted the first genome-wide transcriptomic and gene network analysis of anthers grown under short and long day lengths and revealed functional and regulatory genes that respond to photoperiod.

Photoperiod Changes Gene Expression in Anthers
Approximately one-third of expressed genes (11,726 out of ∼30,000 in different samples) were differentially expressed in anthers grown under SD and LD conditions, with gene expression generally down-regulated under SD conditions (Figure 1A, Supplementary Table 1

, and Supplementary Data Sheet 2).
A total of 177 genes were differentially expressed at all time points (Figure 1C), of which 51 were specifically expressed in the anther and other reproductive tissues (Supplementary Figure 5B and Supplementary Data Sheet 2). Several genes whose mutation is known to cause male sterility were found in this set of 51 genes. OsABCG15 plays a role in the formation of the anther cuticle and pollen exine and is involved in the export of lipid precursors from the tapetum to anther locules (Wu et al., 2014;Zhao et al., 2015). OsLAP6 functions in sporopollenin metabolism and affects bacula elongation to regulate pollen exine formation (Zou et al., 2017). HOTHEAD 1 (HTH1) encodes a glucose-methanol-choline (GMC) oxidoreductase; its RNAi line and mutant both display defective anther walls and aborted pollen with decreased amounts of long-chain fatty acids and cutin . Other genes in this set of 177 DEGs, especially genes that specifically expressed in anthers, may also affect male fertility, particularly in a photoperiodsensitive manner.

DEGs Associated With Transport and Metabolism
Gene ontology enrichment analysis of DEGs at different time points revealed a high proportion of genes associated with transport and metabolism, particularly of ions, carbohydrates, and lipids (Figure 2). Generally, these genes were up-regulated at 12:00 under SD and at 08:00 and 12:00 under LD conditions, indicating that these pathways may be associated with the growth and physiological adaption of plants in response to photoperiod. In stage 11 anthers, mitosis and accumulation of nutrients are the main developmental processes, consistent with our GO enrichment results. Correct ion homeostasis is essential for normal plant growth; accordingly, many genes encoding ion transporters were detected (Supplementary Data Sheet 3). Several potassium transporter family genes were enriched, like HKT, HAK, and AKT family genes (Rodríguez-Navarro and Rubio, 2006). HKT1;4 can also exclude Na + in productive stage of rice (Suzuki et al., 2016). Cation/H + exchanger 14 (OsCHX14), involved in K + homeostasis, may play an important role during rice flowering stages (Chen et al., 2016). OsBOR1 and OsBOR4 DEGs encode boron transporters; pollen from osbor4 mutants have defects in pollen tube elongation (Nakagawa et al., 2007;Tanaka et al., 2013). OsFRDL1 is a MATE (multidrug and toxic compound extrusion) family gene, a citrate transporter of iron, whose mutant exhibits a decrease in pollen viability and grain fertility (Yokosho et al., 2009(Yokosho et al., , 2016. SULTR-like phosphorus distribution transporter (SPDT) encodes a plasma-membrane-localized phosphorus transporter that preferentially allocates phosphorus to grains (Yamaji et al., 2017). Psd1 encodes a non-specific lipid transfer protein (nsLTP) that regulates cell division and elongation; psd1 mutants show photoperiod-thermosensitive dwarfism (Deng W. et al., 2020).
Many DEGs were also found to encode sugar transporters and catabolic proteins, such as monosaccharide transporters (MST1 to MST8), sucrose transporters (SUT1, SUT2, SUT3, and SUT5), and invertases (INV1,INV2,INV3,and INV4). Proteins from these three family genes cover sucrose transport from source to sink. Sucrose is loaded into the phloem by SUTs and unloaded in sink tissues, like anthers or seeds, where sucrose is converted by invertases into hexoses and imported into cells by MSTs (Lim et al., 2006). OsMST4 has been shown to play a role in monosaccharide supply during grain filling , and OsMST5 is expressed in panicles before pollination and may be involved in pollen development (Ngampanya et al., 2003). MST8 is predicted to transport monosaccharides from anther wall layers to the microspore, affecting the unloading process of the sourcesink tissue during male reproduction . SUT1 is predicted to play a role in the uptake of sucrose from the apoplast, and its lack during anther development suppresses pollen germination (Hirose et al., 2010). INV4 is expressed specifically in the anther, down-regulated by cold stress, resulting in abnormal accumulation of sucrose (Oliver et al., 2005).
Other DEGs associated with carbohydrate metabolism could also affect anther viability. Glucan synthase-like 5 (GSL5) encodes a callose synthase, and its mutant showed decreased fertility due to defects in the pollen callose wall . Oryza sativa glucanase 1 (Osg1) encodes a plant β-1,3-glucanase essential for callose degradation during tetrad dissolution; RNAi plants exhibited delayed release of young microspores into anther locules, causing male sterility (Wan et al., 2011). Hexokinase (HXK) family proteins are involved in sugar sensing and signaling , and genes encoding HXK5, HXK7, and HXK10 were among DEGs detected in our analysis. HXK10 can also phosphorylate hexose sugars; RNAi knockdown lines showed non-dehiscent anthers (Xu et al., 2008). Genes encoding disproportionating enzymes DPE1 and DPE2 were also among our DEGs; these proteins are involved in starch metabolism and maturation of starch granules (Akdogan et al., 2011;Hwang et al., 2016).
Differentially expressed genes encoding transport and metabolic genes clustered into two, similarly sized groups that were up-regulated in either SD or LD anthers, suggesting that these biological processes may occur, or be regulated by, different pathways depending on photoperiod to maintain nutrient balances essential for normal development (Figures 3A-C). While previously reported genes affecting male sterility were found among our photoperiod-sensitive DEGs, their mutations generally caused total male sterility rather than PGMS (i.e., Osg1, OsFRDL1, and OsBOR4). Discovery of new potential PGMS regulators will require careful analysis of differential expression of candidate genes under SD and LD conditions.

Differential Expression of Circadian Genes
We analyzed the expression pattern of genes associated with the control of circadian rhythm and found that OsLUX, OsFKF1, OsGI, and OsTOC1 show higher expression in LD over SD anthers at16:00, OsPRR37 was up-regulated under night at LD, while OsPRR95 has higher expression at 12:00. OsCCA1 was upregulated at 12:00 and 16:00 in SD compared with LD anthers (Supplementary Figure 7). These genes and proteins interact with each other in complex ways. In Arabidopsis, CCA1 positively regulates PRR7 and PRR9 expression, which in turn affect feedback regulation of CCA1 (Farré et al., 2005). CCA1 can also interact with TOC1 to bind the GI promoter (Hsu and Harmer, 2014). LUX is a MYB transcript factor required for normal circadian rhythm that can repress PRR9 expression (Chow et al., 2012). FKF1 is an F-box protein that forms a complex with GI to regulate the transition to flowering (Sawa et al., 2007). In rice, overexpression of CCA1 or GI delays flowering (Hayama et al., 2003). The presence of these important circadian genes in our DEGs suggests that they may also play a critical role in anther development in response to photoperiod and warrant further investigation as potential candidates for PGMS line creation.

Protein Metabolism and Phytohormones Play a Key Role in Photoperiod-Sensitive Anther Development
Weighted gene co-expression network analysis (WGCNA), which evaluates potential interaction between all expressed genes, revealed large gene clusters and regulatory networks accociated with photoperiod length and confirmed that gene expression differs significantly between SD and LD anthers. Four modules were identified with a high correlation to photoperiod: genes in turquoise (5852 genes) and magenta (1132 genes) modules were up-regulated in LD anthers, while genes in blue (4929 genes) and yellow (2698 genes) modules were up-regulated in SD anthers (Figure 4 and Supplementary Figure 8).
These modules were generally enriched with the same GO terms as DEGs (Figures 2, 5), but the blue module was also enriched in genes from new categories, i.e., protein-related and macromolecule-related processes, which contained TMS10 (thermosensitive male sterility 10), TMS10L (thermosensitive male sterility 10-like), and UbL404 (ubiquitin fusion ribosomal protein L40) genes. TMS10 and TMS10L are two thermosensitive genes that redundantly regulate post-meiotic tapetal development and pollen development in low and high temperatures (Yu et al., 2017). UbL404 is a target of TMS5, which encodes a conserved RNase Z protein. UbL404 is degraded by RNase Z in wild-type plants, but accumulates in tms5 plants at high temperatures, causing temperature-sensitive male sterility (Zhou et al., 2014). A subsequent mutation in GATA10, which decreases UbL404 levels, can restore fertility at high temperatures (Jin et al., 2019). Our modules thus contain known PGMS and TGMS effectors and are likely to be a good source to uncover as-yet-unknown genes affecting PGMS.
The "protein process" GO category also contained genes associated with light response and anther-specific genes ( Figure 5 and Supplementary Data Sheet 7). Light response genes included PHYB, PHYC, and Hd6. PHYB and PHYC are photoreceptors in rice, single mutants of which induce early flowering; a triple mutant with PHYA causes incomplete sterility (Takano et al., 2005). Hd6 encode a casein kinase II that delays heading time in LD conditions (Takahashi et al., 2001). Anther-specific genes included Pollen Mitosis Relative (PMR), essential for correct pollen mitosis , and DWARF AND RUNTISH SPIKELET (DRUS1) and DRUS2, which regulate anther development and pollen maturation redundantly through sugar utilization or conversion (Pu et al., 2017). Genes involved in the MAPK (mitogen-activated protein kinase) and Ca 2+related signaling pathway were also found, such as Calciumdependent protein kinase 2 (OsCDPK2), which is repressed by light and disrupts seed development (Frattini et al., 1999;Morello et al., 2000). Genes associated with phytohormone signaling included brassinosteroid-signaling kinase (BSK) family members and Ethylene Receptor 2 (ETR2). Thus, while mitosis and sugar transport are key biological processes occurring in stage 11 anthers, hormone, and other signaling pathways are also very active during anther late development.
Our WGCNA analysis revealed 496 hub genes in these four modules (Supplementary Data Sheet 6). A known PGMS gene, CSA, and TGMS gene, Ugp1, were among these hub genes, both known to affect male sterility through carbohydrate metabolism (Chen et al., 2007;Zhang et al., 2013). Other genes affecting carbohydrate metabolism were also among the hub genes, including CAP1, DPE2, and OsSWN1. These results provide further evidence that carbohydrate metabolism is a key characteristic of late anther development affected by photoperiod. Multiple phytohormone-related genes were also present in the hubs: ABA 8'-Hydroxylase gene 3 (OsABA8ox3) is a key gene in abscisic acid (ABA) catabolism (Zhu et al., 2009); OsNAP can be induced by ABA and, in turn, affects ABA biosynthesis (Liang et al., 2014); Gibberellin 20-oxidase 3 (OsGA20ox3) and Gibberellin 2-oxidase 4 (OsGA2ox4) affect gibberellic acid biosynthesis, GA20ox3, in a temperature-dependent manner (Lo et al., 2008;Sakata et al., 2014); Dioxygenase of auxin oxidation (DAO) affects the auxin pathway (Zhao et al., 2013); and SK3/SHAGGY-like kinase (GSK2) acts as a negative regulator of brassinosteroid signaling (Tong et al., 2012). Phytohormones thus appear to be key regulators mediating the adaption of anther development to changing day length.

TGMS and PGMS Share Regulatory Network Elements
Given the presence of CSA and Ugp1 as hub genes, we searched for, and found, other hub genes regulating environment-sensitive male sterility (Table 1 and Supplementary Data Sheet 6). TMS5 was clustered in the red module, which did not show significant correlation with photoperiod. OsPDCD5, Ugp2, and Osgata10 clustered in the turquoise module, while TMS10 clustered in the blue module. Both OsPDCD5 and TMS10 are involved in PCD, which, in rice anthers, begins before stage 11; had younger anthers been used in our analysis, WGCNA significance scores for these genes may have been higher.
Analysis of interaction networks of key hub genes revealed further interactions between PGMS gene CSA and TGMS gene Ugp1 that two genes share large part of co-expressed genes. Although these genes were up-and down-regulated in an opposite pattern CSA and Ugp1, these genes are likely to be common targets of Ugp1 and CSA regulation and have similar roles in co-regulating photoperiod-or temperature-sensitive anther development. This common list included genes encoding zinc finger, MADS-Box, and MYB transcript factors, which may play a key relational role in mediating photoperiod-associated male reproduction. Moreover, carbohydrate metabolic genes and transport genes like invertase were also detected, and these genes may be associated with the sugar transport and metabolism under SD and LD. Further analysis on 38 hub genes highly coexpressed with CSA and Ugp1 may reveal key genes regulating photoperiod-dependent male fertility.
In summary, we provided a transcriptomic view in anther development at different photoperiods and narrowed down candidates not only for future discovery of new PGMS but also for genes that interacted with reported EGMS, CSA and Ugp1, based on DEG analysis and WGCNA. Further investigations to confirm the biological function of these genes will reveal new genetic and molecular controls of PGMS in rice and possible other crops.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, GSE163030.

AUTHOR CONTRIBUTIONS
DZ designed the study. DW, JL, YL, and SS performed experiments. SS analyzed the RNA-seq data and drafted the manuscript. WC, GL, XZ, and DZ revised the manuscript. All authors contributed to the article and approved the submitted version.