GWAS of days to heading (DTH) under different N concentration conditions

To evaluate the effect of N fertiliser on flowering, we measured DTH in 149 Oryza sativa japonica varieties that were not highly genetically structured and interrelated (Supplementary Fig. 1a and Supplementary Data 1) to avoid excessive allelic heterogeneity and increase the sensitivity of signal detection in GWAS. We cultivated the rice population in the field, with three levels of N fertilization (low, normal and high; 0, 9 and 18 kgN/10 a) and measured DTH in 2016. A comparative analysis of DTH among the three N conditions revealed that 64.8% of plants showed delayed flowering under normal versus low conditions, and 78.3% and 82.7% of plants in high showed delayed flowering when compared with normal and low plants, respectively (Fig. 1a). Further, DTH in each variety varied in response to the N condition, suggesting that this population was suitable for the isolation of genetic factor(s) related to flowering time responses to N levels.

Next, we performed GWASs of DTH under the three conditions and observed similar pattern to the previous GWASs using DTH of Japanese rice cultivars, and some peaks were related to previously identified genes, HESO1 and Hd1 (Supplementary Fig. 1b-g). However, there was a significant peak on the terminal end of the long arm of chromosome (Chr.) 3 under the high N condition (Supplementary Fig. 1d). To further clarify the genomic region related to the response to N levels, we conducted additional GWASs considering the difference and ratio of DTH among the conditions (Fig. 1b, Supplementary Figs. 2a-e and 3a-f). We detected a prominent peak, which corresponded to that found in the GWAS under the high N condition, in the GWASs using the difference and ratio of DTH between low and high, and normal and high conditions (Fig. 1b and Supplementary Fig. 2a, c, d), but not in the results for that between low and normal (Supplementary Fig. 2b, e). These results highlighted that the candidate region was significantly correlated to N-responsive flowering regulation. We additionally performed GWASs of principal component (PC) scores derived from PCA with the DTH data (Supplementary Figs. 2f, g and 3g, h). We found a slight but prominent peak on Chr. 3 in a GWAS using PC2, which corresponded to that found in the above GWASs (Supplementary Fig. 2g), suggesting that PC2 mainly represented a response to N levels. Actually, the contribution of PC1 of the PCA was higher than 98% and that of PC2 was less than 2% and the factor loadings of the N conditions were different in PC2 but not in PC1 (Supplementary Fig. 2f).

Similar analyses for DTH data of 128 varieties grown under the same three N levels in 2017 also revealed a non-significant but prominent association between the candidate region and DTH difference and ratio under low vs. high and normal vs. high conditions (Supplementary Figs. 4a-j and 5a-i and Supplementary Data 2). The GWASs of PC scores from the analyses using the DTH data of 2017 roughly replicated those in 2016 (Supplementary Figs. 4k-m and 5j, k). These results strongly implied that the peak on Chr. 3 is mainly associated region with delayed flowering under high N levels.

To identify the candidate region on Chr. 3, we calculated pairwise linkage disequilibrium (LD) correlations around the peak and found an LD block from 30.30 Mb to 31.8 Mb (Fig. 1c). The candidate region of the peak included three polymorphisms with significant associations and amino acid exchanges or stop-loss variants (Supplementary Data 3). One of the polymorphisms was a stop-loss variant in LOC_Os03g55389 (Hd6), which encodes CK2α. Because Hd6 reportedly plays a major role in the photoperiodic control of flowering, we focused on this gene as a candidate gene of N-responsive flowering regulation. There were two haplotypes in the coding region of Hd6 in our population: haplotype (Hap) A and Hap B. Hap A reportedly is a loss-of-function type because of a premature stop codon, whereas Hap B is functional (Fig. 1d). We compared DTH at the different N levels between the Hd6 haplotypes (Fig. 1e). In all N conditions, DTH was shorter for Hap A than for Hap B. DTH was longer under high than under low or normal in Hap B, while differences among the various N conditions were not observed in Hap A, which was consistent with the GWAS results. These results suggested that the functional Hd6 allele is required for N-responsive flowering regulation.

Studies have demonstrated that Hd6 requires functional Hd1 and Hd2 to control flowering time in rice. Therefore, we investigated the epistatic interactions among Hd1, Hd2 and Hd6 in N-responsive flowering regulation by comparing DTH among functional (+/+) and non-functional (-/-) alleles in our population (Fig. 2a-c). Apparent differences were not observed among the three N levels, except for a combination of all functional alleles (Hd1 Hd2 Hd6). This result supported that an epistatic interaction among Hd1, Hd2 and Hd6 is essential for N-responsive flowering regulation. To confirm this genetic relationship, we compared the DTH of five near-isogenic lines (NILs; NIL-Hd1/Hd2/Hd6, NIL-Hd2/Hd6, NIL-Hd1/Hd6, NIL-Hd1 and NIL-Hd2) and Nipponbare (Hd1/Hd2) under normal and high conditions (Fig. 2d, e). NIL-Hd1/Hd2/Hd6 (Hd1 Hd2 Hd6) showed a delay in DTH under H compared to N, whereas the other NILs did not show such N response. These results demonstrated that Hd6 was a causal gene of the peak detected in GWAS for delayed flowering under high-N condition, and that the epistatic interaction among Hd1, Hd2 and Hd6 is a prerequisite for N-responsive flowering regulation.

To substantiate the genetic evidence at the molecular level, we first analysed the effect of N application on Hd6 expression. Consistent with previous observations, exogenous application of various concentrations of NHCl as an N source to rice seedlings in hydroponic culture induced the expression of known N metabolic genes, including GS1;2 and AS1 (Supplementary Fig. 6a, b). However, Hd6 expression did not respond to the N treatment, and mRNA levels of Hd1 and Ghd7 tended to be reduced after N supply (Fig. 3a and Supplementary Fig. 6c, d). To clarify the interaction between Hd6 protein levels and N, we generated an anti-Hd6 antibody to detect endogenous Hd6 protein expression (Supplementary Fig. 7). Hd6 protein levels in NIL (Hd6) seedlings harbouring the functional Hd6 allele clearly showed a N dose-dependent increase, which peaked around 1.6 mM NHCl in this experimental condition (Fig. 3b). These results indicated that N regulates the amount of Hd6 protein, not its expression.

Since, as mentioned above, a genetic interaction between Hd6 and Hd2 has been revealed and our results indicated that Hd2 is required for N-responsive flowering regulation along with Hd6, we hypothesised that N may modulate the Hd6-Hd2 module. Like Hd6 expression, endogenous Hd2 expression is not affected by different N concentrations (Fig. 3c). To investigate Hd2 protein function in vivo, we analysed transgenic rice plants harbouring either a Ubq:Ω:3FLAG-Hd2 or a Ubq:Ω:3FLAG-Hd2 construct in tandem with CaMV35S-Hd6 (35S:Hd6/Ubq:Ω:3FLAG-Hd2) in a NIL (hd6 hd16) with natural variation for defective alleles of both Hd6 and Hd16 in the Nipponbare background. The presence of Hd6 protein significantly enhanced the immunological detection of 3FLAG-Hd2 in transgenic plants using anti-FLAG antibody, indicating that Hd6 protein enhances 3FLAG-Hd2 protein accumulation in vivo (Fig. 3d). 3FLAG-Hd2 clearly accumulated in response to N dosage in the presence of Hd6 (Fig. 3d). Artificial overexpression of Hd6 did not increase 3FLAG-Hd2 transcript levels (Fig. 3e, f). These results indicate that the turnover of 3FLAG-Hd2 protein is attenuated by N in an Hd6-dependent manner.

Consistent with the roles of Hd6 and Hd2 in the rice photoperiodic flowering pathway, the magnitude of floral repression was significantly and additively enhanced in the presence of Hd2 and Hd6 (Fig. 4a and Supplementary Fig. 8). Among these lines, only 35S:Hd6/Ubq:Ω:3FLAG-Hd2 plants showed the effect of high-N on the delay of flowering (Fig. 4a and Supplementary Fig. 9), confirming that Hd1, Hd2 and Hd6 are essential requisites for N-responsive flowering regulation as suggested by the experimental results obtained with the Japanese varieties and NILs (Fig. 2c, e), because Nipponbare harbours functional Hd1 and Hd2 and non-functional Hd6 (Supplementary Data 1).

To further analyse the repressive action of the Hd6-Hd2 module in N-responsive flowering, we conducted RNA-sequencing analyses of the transgenic plants under low and high conditions at ZT0 when the transcripts of the floral inducer genes are most abundant (Fig. 4b and Supplementary Fig. 10). We identified 409 differentially expressed genes (DEGs) that were unique to 35S:Hd6/Ubq:Ω:3FLAG-Hd2 plants in response to the N condition. Among these, 306 DEGs were significantly upregulated, and 103 were significantly downregulated in the transgenic plants under the high condition. Gene Ontology (GO) analysis of the unique DEGs in each transgenic line revealed that the 6 genes among the 103 DEGs downregulated in 35S:Hd6/Ubq:Ω:3FLAG-Hd2 plants were related to flowering and photoperiodism, whereas these GO terms were not found for genes upregulated in 35S:Hd6/Ubq:Ω:3FLAG-Hd2 plants and for DEGs in NIL or Ubq:Ω:3FLAG-Hd2 plants (Fig. 4c and Supplementary Fig. 11). Transcriptome profiling and quantitative reverse transcription (RT-q)PCR assays revealed significant decreases in the transcript levels of genes related to flower development under high condition only in 35S:Hd6/Ubq:Ω:3FLAG-Hd2 plants (Fig. 4d and Supplementary Data 4). These genes included Ehd1, Hd3a, RFT1, OsMADS14, OsMADS18 and LHS1/OsMADS1, which are major floral inducer and floral development genes. Based on these results, we concluded that N signalling uses the Hd6-Hd2 module to repress floral induction in the rice photoperiodic flowering pathway.

To clarify the biochemical relationship between Hd6 and Hd2, we further analysed the biochemical properties of 3FLAG-Hd2 protein using 35S:Hd6/Ubq:Ω:3FLAG-Hd2 plants. 3FLAG-Hd2 accumulation was cancelled in the presence of the CK2 inhibitor benzimidazole but not the CK1 inhibitor PF60462 (Fig. 4e), indicating that Hd6-dependent phosphorylation is essential for the stabilisation of the Hd2 protein. 3FLAG-Hd2 protein tended to produce a broad band, suggesting the detection of the phosphorylated Hd2 during this growth stage. As expected, the upper portion of the 3FLAG-Hd2 signal in immunoblots disappeared after protein treatment with alkaline phosphatase, confirming the phosphorylation (Fig. 4f). Furthermore, immunoblotting analysis using Phos-tag, which specifically captures phosphate residues, demonstrated that phosphorylated 3FLAG-Hd2 protein was super-shifted under the high N condition (Fig. 4g), suggesting 3FLAG-Hd2 phosphorylation under this condition. Mass spectrometric analysis of a purified recombinant Hd2 (rHd2) pre-incubated with or without recombinant Hd6 (rHd6) revealed a phosphorylation site in rHd2 by rHd6 on the serine at position 369 (S369) (Supplementary Fig. 12a-c), which is consistent with a previous study showing that the middle region of Hd2 can be phosphorylated by Hd6. These results demonstrated that N-induced Hd6 protein accumulation enhances the phosphorylation of S369 of Hd2 protein.

Given that Hd1 function is epistatic to Hd2 function in the rice photoperiodic flowering pathway, floral repression conferred by Hd2 protein stabilised by Hd6 protein accumulated under high-N condition may be integrated into Hd1 action. The genetic relationships among Hd6, Hd2 and Hd1 have been revealed, but the molecular relationships remain unknown. Therefore, we analysed the physical relationships among these proteins. In rice mesophyll protoplasts, transiently expressed GFP-Hd6, GFP-Hd2 and GFP-Hd1 were observed in the nucleus (Fig. 5a). Hd6 was highly accumulated in the nucleolus, whereas Hd2 and Hd1 were not observed in the nucleolus (Fig. 5a and Supplementary Fig. 13a). A bimolecular fluorescence complementation (BiFC) assay revealed the presence of Hd6-Hd2 and Hd6-Hd1 complexes in the nucleus, and YFP signals were accumulated in the nucleolus (Fig. 5b, c and Supplementary Fig. 13b). The interaction of Hd2 and Hd1 was not observed in the presence or absence of Hd6. We confirmed the direct association of N-induced Hd6 protein produced in NIL plants expressing Hd6 with a recombinant Hd2 protein (rHd2) synthesised by E. coli, which indicated that the capacity of Hd6 proteins to form a complex with Hd2 in plants increases in response to exogenous N levels (Fig. 5d). We also performed domain analysis of Hd2 for interaction with Hd6 by BiFC assay. The Hd2 protein has two conserved domains: pseudo receiver (PR) and CCT (Supplementary Fig. 12a). Truncated Hd2 protein containing PR but not CCT (Nterm and ΔMCCT) interacted with Hd6 in the cytosol and truncated Hd2 containing the middle region (300-480) (Mreg and ΔCCT) interacted in the nucleus, whereas truncated Hd2 containing only the C-terminal region, including the CCT domain (Cterm), could no longer interact (Supplementary Fig. 13c). These results indicated that the PR domain and the middle region of Hd2 are essential for Hd6-Hd2 interaction and determine the subcellular localisation of the complex. In vitro pull-down assays using recombinant proteins produced in E. coli confirmed the Hd6-Hd2 interaction, whereas only a faint signal was observed for Hd1-Hd2 interaction and none for Hd1-Hd6 interaction (Supplementary Fig. 13d). Therefore, the physical interaction of the Hd6-Hd2 complex may be stable, whereas that of other combinations may be not or depend on the conditions.

As Hd1 protein can bind directly to CO response element 2 (CORE2) found in the Hd3a promoter proximal region, we developed an in vivo transient reporter assay system using rice protoplasts prepared from rice seedlings to investigate the potential roles of Hd6 and Hd2 protein in CORE2-dependent Hd3a transcription in the presence of Hd1 protein. We found that Hd1 significantly activated firefly luciferase (LUC) expression driven by an artificial promoter, i.e., the 35S promoter inserted four times of CORE2-containing oligonucleotide in front of the minimal promoter (MP) (Fig. 5e, f). The repressive effect of Hd2 on the CORE2-containing promoter was observed upon single transfection of Hd2 (Supplementary Fig. 14), and co-transfection of Hd1 and Hd2 significantly reduced the Hd1 solo activation of the CORE2-containing promoter (Fig. 5f). The counteractive effect of Hd2 on Hd1 was dose-dependent (Fig. 5g). Therefore, Hd2 protein antagonises the CORE2-dependent action of Hd1 protein. In other words, the regulation of Hd3a expression via Hd1 and Hd2 is mutually exclusive. Furthermore, Hd6 can enhance the repressive effect of Hd2 (Fig. 5f), which is consistent with the accumulation of Hd2 in the presence of Hd6 (Fig. 3d). These results demonstrated that N-induced Hd6 protein accumulation enhances the phosphorylation of Hd2 protein, thus stabilising it for the regulation of floral inducer gene expression together with Hd1 (Fig. 5h).