Genome-wide association study and gene network analyses reveal potential candidate genes for high night temperature tolerance in rice

High night temperatures (HNT) are shown to significantly reduce rice (Oryza sativa L.) yield and quality. A better understanding of the genetic architecture of HNT tolerance will help rice breeders to develop varieties adapted to future warmer climates. In this study, a diverse indica rice panel displayed a wide range of phenotypic variability in yield and quality traits under control night (24 °C) and higher night (29 °C) temperatures. Genome-wide association analysis revealed 38 genetic loci associated across treatments (18 for control and 20 for HNT). Nineteen loci were detected with the relative changes in the traits between control and HNT. Positive phenotypic correlations and co-located genetic loci with previously cloned grain size genes revealed common genetic regulation between control and HNT, particularly grain size. Network-based predictive models prioritized 20 causal genes at the genetic loci based on known gene/s expression under HNT in rice. Our study provides important insights for future candidate gene validation and molecular marker development to enhance HNT tolerance in rice. Integrated physiological, genomic, and gene network-informed approaches indicate that the candidate genes for stay-green trait may be relevant to minimizing HNT-induced yield and quality losses during grain filling in rice by optimizing source-sink relationships.

Rice is a major source of calories for more than 50% of the world's population and plays a key role in sustaining food security in Asia 1 , Africa, and the Americas 2 . Current global climate models predict a mean temperature increase by 1.0-3.7 °C by 2100 3,4 , which would potentially increase the occurrence of heat stress events in ricegrowing tropical and subtropical countries. Disaggregating the predicted mean temperature increase has allowed unraveling a rapid rise in night temperature compared to the maximum day temperature at the farm 5,6 , region 7,8 , and global 9 scale. A rapid rise in night temperature is fast becoming a significant challenge for sustaining crop yields, including rice 5,10-13 , wheat (Triticum aestivum) 14,15 and barley (Hordeum vulgare) 16,17 . Previous studies have shown that an average season-long increase in minimum night temperature by 1 °C reduced tropical irrigated rice (Oryza sativa L.) grain yield by 10% 5 . Additionally, field-based studies document the susceptibility of rice inbred cultivars or commercial hybrids exposed to HNT (28 or 29 °C), leading to an average grain yield reduction of ~ 15% 13,18 to 90% when exposed to severe HNT (32 °C) under controlled environments 19 . An increase in night temperatures have resulted in a negative impact on grain yield and quality in the major rice-growing regions of South and Southeast Asia 5,7 and the United States 20 . Therefore, breeding programs need to focus on enhancing high night temperature (HNT) stress tolerance in rice to minimize yield and grain quality losses.
Unlike the pre-flowering phase, the post-flowering phase is most affected by HNT, wherein the rate of senescence is advanced leading to lower grain-filling duration, resulting in lower grain yield and quality 13,14 . In response to high day temperature (HDT), rice can benefit from early morning flowering mechanism to escape or effectively employ transpiration cooling to avoid heat stress during the day 21,22 . However, rice is not equipped with similar physiological processes that can minimize the negative impact of HNT [22][23][24] . Interestingly, the physiological routes through which HDT and HNT stress induce yield losses are shown to differ, although both the stresses finally result in significant yield and quality losses 22 . Experiments under controlled environments [25][26][27][28][29] and field conditions 30,31 have convincingly demonstrated that HDT stress negatively impacts reproductive processs in rice, leading to spikelet sterility. However, induction of HNT from panicle initiation until physiological maturity under field conditions did not induce a significant reduction in spikelet sterility and had no loss in seed number 11 . The same response was captured even with different levels of nitrogen application 11,32 and across a range of inbreds and hybrids 13,18,33 , providing strong support for the differential route through which HNT induced yield and quality losses in rice. These results make a strong case that progress achieved through genetic mapping for HDT stress in rice may not be directly translated to induce HNT tolerance in rice 22 .
Rice exposed to HNT stress during the grain-filling stage significantly reduces whole plant and panicle carbohydrate content due to increased night respiration (carbon loss) resulting in lower biomass, lower grain yield (due to lower 100-seed weight), and enhanced chalkiness 11,13,15,34 . Grain yield and quality, specifically chalkiness, are key traits that drive market value and are considered significant aspects for assessing HNT impact in rice 10,13 . Despite the increasing knowledge of the metabolomic 35 and transcriptomic 36,37 responses of rice to HNT, our understanding of the genetic basis of rice response to HNT tolerance is mostly unexplored. Hence, to the best of our knowledge, this is the first study to identify genetic loci and putative candidate genes controlling HNT response in rice that will help in developing HNT tolerant varieties and hybrids for current and future warmer climates.
To fill the above knowledge gaps, an indica panel was phenotyped for grain yield, 100-seed weight, harvest index, grain size and chalkiness, under HNT using controlled environment walk-in chambers and field-based heat tents. The diversity panel consisting of improved and traditional genotypes with sub/tropical adaptation was phenotyped to (1) identify genomic regions associated with phenotypic traits in control and HNT in controlled environments and field conditions; (2) compare the extent of co-localization of the genetic loci across treatments under controlled environment and field conditions, and (3) identify putative candidate genes associated with enhanced HNT tolerance in rice using gene network analysis.

Materials and methods
Plant materials. The indica rice diversity panel (PRAY-Phenomics of Rice Adaptation and Yield) involving cultivars from tropical and subtropical rice-growing regions were assembled as a part of the Global Rice Phenotyping Network (GRPN, http:// ricep henon etwork. irri. org). A sub-set of two-hundred and nine genotypes was chosen based on the uniformity in flowering phenology, using data generated under the International Rice Research Institute (IRRI) farm conditions 38 . Two independent experiments were conducted during 2015 (Experiment 1-controlled environment facility) and 2016 (Experiment 2-field-based heat tents) at IRRI, Los Baños (14° 11′ N, 121° 15′ E, 21 masl), Philippines. In Exp. 1, the 209 genotypes with 3 replicates and two treatments resulted in 1254 pots. In Exp. 2, the same 209 genotypes with each grown in a row (each row having 15 plants) and two treatments resulted in about 6270 plants ( Supplementary Fig. S1). Plants in both experiments were grown under fully flooded conditions to study the effect of HNT stress starting from panicle initiation until physiological maturity. Unlike short episodes of high day heat stress, HNT is acknowledged to persist for longer growth periods encompassing different developmental stages, leading to yield and quality losses 11,13 . Previous studies using multiple genotypes with HNT stress imposed in the same field-based tents concluded that the impact of HNT was significant when coincided with panicle initiation and lasting until physiological maturity, but not during the early vegetative stage 33 .

Experiment 1 (walk-in chamber facility)
Crop management. Seed dormancy of rice genotypes was broken by exposing them to 50 °C for three days and were then hand sown in seeding trays. Fourteen-day-old seedlings were transplanted into plastic pots (25 cm in height; 26 and 20 cm diameter at the top and bottom, respectively), filled with 6 kg of clay loam farm soil and 2.0 g ammonium sulfate [(NH 4 ) 2 SO 4 ], 1.0 g single superphosphate (SSP) and 1.0 g muriate of potash (KCl). An additional 2.5 g (NH 4 ) 2 SO 4 was top-dressed at 30 days after transplanting. One plant per pot was grown under fully flooded conditions and was grown in the greenhouse conditions before imposing stress at the targeted stages (panicle initiation until physiological maturity). The 209 genotypes were grouped into five batches based on flowering time information from Kadam et al. 38 . The five batches were stagger sown with a 5-day interval between each batch (Fig. 1a) and transplanted with the same intervals to ensure best possible synchrony in panicle initiation and flowering. Representative plants for each batch were monitored and used to determine the panicle initiation stage by destructive inspection 39 .  Fig. S1a). Each chamber was fitted with six independent units of 1 kW high-intensity discharge lamps, providing a photosynthetic photon flux density (PPFD) of ~ 650 µmol m −2 s −1 at 1.2 m from the base of the chamber. A photoperiod of 11 h (13 h of dark) and a constant relative humidity of ~ 70% was maintained throughout the stress period. The HNT temperature in the chambers was maintained at 29  . All the fertilizers were incorporated into the plots as a basal dose (1 day before transplanting) except for urea, which was applied in four splits of 3:2:3:2 at basal, mid-tillering, panicle initiation, and heading stage, respectively. All the agronomic and crop management practices were followed, similar to Shi et al. 11 and Bahuguna et al. 13 .

Walk-in chambers and
HNT stress using field-based tents. All genotypes were exposed to two-night temperature treatments i.e., control night (CNT; 24 °C) and high night temperatures (HNT; 29 °C) with a common day-time temperature across both treatments. Plants were exposed to 24 °C (CNT) or 29 °C (HNT) night temperature by using field-based tents between (18:00 to 06:00 h) and during the day-time (06:00 to18:00 h) the tents were open, exposing the plants to natural conditions [ Supplementary Fig. S1]. To start the HNT stress, the sides of the tents were manually rolled down at 18:00 h every evening and re-opened at 6:00 h in the morning, exposing the plants to12 h of photoperiod. Each tent was outfitted with an air conditioner, which was programmed to automatically regulate control (24 °C) and stress (29 °C)  Non-stress period High night temperature Treatment A Figure 1. Synchronization of the reproductive stage by stagger sowing (A) for phenotyping indica panel under HNT stress (B) in walk-in chambers (Exp. 1) and field-based heat tents (Exp. 2). DFF days to 50% flowering, S seeding, T transplanting, PI panicle initiation, FL flowering, GF grain filling, PM physiological maturity. Amplitude is calculated as a difference between average day temperature and average night temperature for each treatment separately. VPD vapor pressure deficit, RH relative humidity and SD standard deviation. Schematic representation of materials and methods followed to impose HNT stress is visualized in Supplementary Fig. S1 www.nature.com/scientificreports/ night (18:00-06:00 h) temperature during stress treatment (from panicle initiation until maturity) was 24.3 °C (SD = 0.7) for the control and 28.3 °C (SD = 0.53) with HNT in Exp. 2 (Field experiment; Fig. 1b). Amplitude (the difference between day and night) temperatures during the stress treatment was 5.8 °C (Fig. 1b).
Phenotypes and statistical analysis. At physiological maturity, about three (Exp. 1) to twelve plants (Exp. 2) were harvested separately to determine the impact of HNT on the grain yield per plant (YPP) and 100grain weight (HWT), harvest index (HI), and chalkiness (CHALK). Panicles were hand-threshed to estimate grain yield per plant and 100-grain weight. The harvest index was calculated as the ratio of grain yield to the plant's total dry matter. In Exp. 2, representative samples from each genotype were analyzed for grain size and grain chalkiness at the Grain Quality and Nutrition Center, IRRI, Philippines, following the ISO 17025 certified protocols of the IRRI GQNSL (http:// www. knowl edgeb ank. irri. org/ riceb reedi ngcou rse/ bodyd efault. htm# Grain_ quali ty. htm). All phenotypic traits were subjected to statistical analysis to detect the effects of genotype and treatment using the library ("agricolae") RStudio 3.6.1 (https:// rstud io. com/). Least square difference (LSD) was used to compare the mean values differences between control and HNT. Pearson's correlation coefficient was carried out between all traits using the library ("corrplot") in RStudio 3.6.1. The relative performance of a genotype under stress was determined using the following equation.

Genome-wide association analysis (GWAS).
One hundred and seventy-three genotypes in Exp. 1 (panicles did not emerge in 36 genotypes under HNT), and 206 genotypes in Exp. 2 had complete phenotypic data. However, the genotypic data for twenty genotypes were missing, and hence only those genotypes that had both genotypic and phenotypic data (158 genotypes in Exp. 1 and 188 in Exp. 2) were used for GWAS analysis. GWAS was performed to estimate the association between marker-trait using the Genomic Association and Prediction Integrated Tool 40 . A mixed-model approach was implemented in the Compressed Mixed-Model Association (CMLM), which controls the familial relatedness and population structure and kinship (as a random effect to control population structure) to minimize false-positive associations 41 . A total of 45, 200 single nucleotide polymorphisms (SNPs; minor allele frequency > 0.05) that were generated under the Global Rice Phenotyping Network (http:// ricep henon etwork. irri. org) were used for performing GWAS 38,42,43 . Estimated absolute phenotypic or relative values along with PC and Kinship were used to discover marker-trait associations (MTAs). The Bonferroni correction had a more stringent threshold; therefore it did not result in many significant associations, so we have used a significant (-log(p) = 3.5) threshold to discover MTAs similar to several previous studies 42, 44 . Manhattan plots were constructed using the library ("qqman") in RStudio 3.6.1. (https:// cran.r-proje ct. org/ web/ packa ges/ qqman/ vigne ttes/ qqman. html). The SNP with the lowest p-value within 100 kb is considered as a candidate marker-trait association (cMTA) to represent that locus. Thus, the SNP with minimum p-value for each locus was defined as cMTA. Then we compared cMTAs/peak SNP to detect shared loci and unique loci across treatments and experiments. Candidate genes were searched in ± 50 kb from the cMTAs in rap-db on IRGSP-1.0 (http:// rapdb. dna. affrc. go. jp/).

Prioritizing gene candidates.
Part of this study was aimed to reduce the inherent bias in selecting candidate genes based on existing literature by incorporating alternate ways to identify the top genes of interest 45 . To the 638 candidate genes identified in the 100-kb window surrounding the unique 51 cMTAs, three potential gene prioritization approaches were applied: (1) Comparison to genes previously identified as responsive to HNT via differential expression analysis, (2) identification of genes highly connected to transcriptional regulators through a gene regulatory network, and (3) genes orthologous to genes previously identified as associated with traits through QTG-Finder2 46 . These three comparisons integrate data from different experimental approaches to narrow down a list of potential candidate genes.

Prioritization based on differential expression analysis
This approach used altered expression in response to HNT as a filter to prioritize candidate genes. A prior experiment performed RNA-Sequencing analysis on the O. sativa genotype IR64, grown under similar experimental conditions at IRRI, Philippines (data accessible at NCBI GEO database 37 , accession GSE159073). RNA-Seq data from this prior time course study from flag-leaf tissue was normalized, mapped, and expression count per gene obtained as described in Desai et al. 37 . Only 354 genes, 223 in the genes identified from the absolute analysis and 131 from the relative analysis, were in the candidate gene list from trait loci and had detectible expression levels in the initial transcriptomic analysis. The expression data were analyzed using the R package edgeR 47 to identify genes differentially expressed at any time point. Candidate genes were prioritized when identified as significantly differentially expressed (FDR adjusted p-value < 0.05) for any of the comparisons between HNT and control.
Prioritization based on connectedness in a gene regulatory network. We used a previously published gene regulatory network created using expression data from different rice developmental stages 37 . The connectedness of transcriptional regulators and their target genes in the network was visualized and examined using the R package networkD3 48 . In this network, the confidence that a transcription factor interacts with a target gene is based on the strength of the association. Interaction scores ranged from 0 (the transcription factor is unlikely to interact with the target) to 871 (a high-confidence interaction www.nature.com/scientificreports/ with strong associations. Only 0.3% of the interactions passed this stringent filter, leaving 45,879 interactions between transcription factors and target genes. The 638 annotated genes identified in the GWAS were compared to these networks. The gene candidates were evaluated based on the strength and number of interactions it had with transcriptional regulators. Orthologous genes. The last approach used to prioritize candidate genes used the relationship to known causal genes to determine the best candidate genes. Lin et al. 46 developed a machine learning approach to prioritizing association-identified loci called QTG-Finder. QTG-Finder takes a list of QTL loci and their associated potential candidate genes from a locus (or loci) and ranks each gene in order from the most to least likely to be a causal gene for each respective locus. This approach has been tested and works well in O. sativa and A. thaliana 46 . This has recently been extended in QTG-Finder2 to use orthologous genes to identify causal genes in Setaria and Sorghum based on rice and Arabidopsis studies 46 . The list of candidate genes from the GWAS study was submitted to QTG-Finder2 to rank the genes per locus. The genes were translated from RAP_ID to MSU_ID format to be compatible with the software. This reduced the total to 523 genes, 333 candidate genes from 34 absolute loci and 190 genes from the 20 relative loci. QTG-Finder2 could not be run on two loci due to an error we could not resolve. Of the remaining loci, the top-ranked genes for each locus were considered to be a prioritized candidate gene. Top-ranked genes were the top two genes for each locus ranked by QTG-Finder2.

Genome-wide association study
Experiment and treatment specific MTAs. We identified genomic regions associated with different traits under control (representing controlled non-stress night temperature), HNT (representing high night temperature) and relative (representing genotype response as a ratio of HNT to CNT), by CMLM (Supplementary Figs. S2 and S3). A summary of GWAS results is given in Table 1, and the list of peak SNPs or candidate MTAs (cMTAs) for each locus are presented in Table 2. Significant MTAs varied between traits (absolute and relative values), from a minimum of five to a maximum of 25 ( Table 1). The phenotypic variance (R 2 ) explained by each marker ranged from 0.07 (relative yield in Exp. 2) to 0.37 for chalkiness (in Exp. 2) under HNT (  (Table 2). Twenty candidates MTAs were identified for the yield and quality-related traits i.e., yield (5), 100-seed weight (3), harvest index (6), chalkiness (4) and one each for grain width and grain length under HNT across experiments (Table 2). We also noticed that the identified genetic loci for grain width (Chr5pos5366489.1) and grain length (Chr3pos16794259.1) co-localized on the well-characterized, previously identified genes across treatments (Fig. 4). Three co-localizing cMTAs were detected between relative and absolute values for grain yield (Chr1pos22255654.1, Chr2pos34468097.1) and harvest index (Chr11pos3631830.1).
Harvest index has been used as a potential proxy for partitioning capacity between source and sink both under non-stress and stress conditions. Six loci (Q16-Q21) were detected for harvest index under control. Minor alleles     (Table 2). In Exp. 2 the minor allele of Q26 had a negative impact on the harvest index, while the minor allele at Q27 affected the same trait positively under HNT (Table 2). Genomic regions (Q22-Q27) associated with harvest index under HNT were in close proximity to the region of genes controlling plant nutrient status (phloem sucrose transport, nitrogen metabolism, photoperiodic control of flowering, osmotic stress response, and lipid homeostasis) (Supplementary Tables S2 and S3). In Exp. 1, GWAS identified a key region on chromosome 2 (34.46 Mb) that was commonly associated with harvest index (Q23) and yield (Q6) under HNT ( Table 2). Out of four cMTAs detected across experiments for relative harvest index, one candidate SNP (Chr11pos3631830.1) was shared between absolute and relative harvest index in Exp. 1 ( Table 2). Chalkiness is one of the primary factors determining rice milling and appearance qualities, both under control and stress conditions 11 . Eleven cMTAs (3 in control, 4 in HNT and 4 for relative) were identified for chalkiness in Exp. 2 (Table 2, Supplementary Fig. S3). Under control conditions, GWAS identified locus (Q30, Chr11pos22646338.1) on chromosome 11 (22.6 Mb) that was highly associated with chalkiness ( Supplementary  Fig. S3b), but not for HNT ( Supplementary Fig. S3c). The minor allele of the lead SNP (Chr3pos1662112.1) underlying the peak Q33 had a negative effect (− 17.8%) on chalkiness ( Table 2). Another cMTA (Q31) for chalkiness had the highest positive allele effect (39.6%) under HNT, detected on chromosome 1. A minor locus (Q34) on chromosome 4 (Chr4pos20331507.1) induced chalkiness (effect = 24.3%) under HNT. Also, one genomic region on chromosome 4 (19.9-20.3 Mb) for chalkiness was found to be co-localized both under control and HNT conditions (Supplementary Fig. S3b, c). Among the identified loci for relative chalkiness in Exp. 2, a SNP (Chr2pos18751285.1) on chromosome 2 was strongly (p = 6E−07) associated with relative chalkiness ( Table 2). The Chr2pos18751285.1 SNP marker's minor allele had a positive effect on relative chalkiness ( Table 2).
Grain shape/size (GS) is usually evaluated by grain width (GW)and grain length (GL), which are the key grain quality traits that influence the rice consumer preferences. Interestingly, a highly significant association (p = 1E−06 and p = 7E−05) was detected for grain width (5.3 Mb) on chromosome 5 across treatments (CNT, and HNT, respectively) (Fig. 4c,d), which is in the region of a well-characterized gene, qSW5/GW5 (encoding a calmodulin-binding motif family protein) that regulates grain width in rice 49,50 . No significant common loci were detected for grain width between absolute and relative values (Fig. 4e). Similar to absolute grain width, the SNP (Chr3pos16794259.1) associated with grain size (GS3) 51,52 , was detected as most significantly associated locus with absolute grain length across treatments (Fig. 4f,g). No loci were shared between the absolute and relative grain length. However, one SNP (Chr3pos33471970.1) on chromosome 3 was significantly associated with relative grain length but not co-located with the GS3 locus (Fig. 4h).
Overall, we found two prominent GWAS loci peaks on chromosome 2 (region 1: 31.4-31.9 Mb and region 2: 34.6-34.7) influencing grain yield and harvest index in rice ( Supplementary Fig. S2). Region 1 (31.4-31.9 Mb) is associated with Q2, Q3, Q8, and Q12 that control yield and 100-seed weight either under control or HNT ( Table 1, Supplementary Fig. S2). Similarly, region 2 (34.5 on chromosome 2; Q6 and Q23) commonly associated  Exp. 2, C). In experiment 1 (walk-in chambers), temperature treatments included control night temperature (CNT) and high night temperature (HNT), and in experiment 2, two different night temperature treatments i.e., control night (CNT) and high night temperature (HNT) were imposed with common day-time ambient temperature across treatments. YPP-grain yield (g per plant); HWT-100-seed weight (g); HI Harvest index and CHALK chalkiness (%). *p < 0.05, **p < 0.01, ***p < 0.001 indicate significant correlations among traits. Prioritizing gene candidates. An in-silico search for all annotated genes in a 100-kb window surrounding the unique 51 candidate SNPs yielded 638 (396 for absolute, 209 for relative and 33 shared between absolute and relative) annotated genes, ranging between 2 (Q20) and 23 (Q19) genes ( Table 2, Supplementary Table S2). We found that these genes encode several classes of protein with known or predicted functions in abiotic and grain development pathways, including helix-loop-helix DNA-binding domain-containing proteins, similar to starch synthase, GRAS transcription factor domain-containing protein, a calcium-binding protein, mitogen-activated protein kinase 4, ROS scavenging, sugar signaling and metabolism, and auxin-responsive protein (Supplemental Table S2). Further, a search for characterized genes within candidate loci yielded 127 functional candidate genes (http:// funri cegen es. ncpgr. cn/) with a range of 1-17 candidate genes per locus (Supplementary Table S3). www.nature.com/scientificreports/ A common approach is prioritizing candidates for each locus based on functional roles from prior literature and coding variants. However, the negative effects of HNT manifest through different mechanisms than HDT, and little is known about the underlying genes driving the genetic architecture of these traits under HNT. Therefore, previous literature may not cover the genes involved in HNT response regulation. Moreover, since non-coding variants can play critical roles in trait expression, prioritizing coding variants ignores significant non-coding regulatory variants. Therefore, three independent approaches were applied to prioritize candidate causal genes driving the traits' genetic architecture under both control and HNT conditions. The annotated genes in the 100-kb window surrounding each locus were examined for overlap in three ways: (1) with genes differentially expressed in response to HNT, (2) strength of association in a gene regulatory network, and (3) genes orthologous to prior causal genes. This approach reduces the bias by prioritizing candidate genes based on other data sets.
The first approach, which examined if candidate genes were differentially expressed in HNT, found 34 significant candidate genes, 22 of the genes from loci identified in the absolute analysis and 12 from the relative loci. The gene regulatory network approach prioritized 76 genes based on their interaction with transcription regulators, 50 from the absolute candidates and 26 from the relative candidates (Fig. 5a,b). The final method, which prioritized candidate genes based on their orthogonality to known causal genes, ranked 97 genes as top priority candidates, 66 absolute and 34 relative. To further prioritize and select the most promising candidates, each approach's selected genes were crosschecked for their performance in the other two techniques. Candidates were considered highly promising if they had two methods supporting them. Additionally, since many of the genes could not be tested in one or more of the approaches, candidates could still be considered highly promising if supported by one method. However, they could not be tested in both methods. For example, Os02g0635600 ranked highest in Q22 (harvest index in absolute HNT) based on ortholog to known causal genes using QTG-Finder2, but this gene was not present in the full gene regulatory network and was not detected in the transcription data (Table 3). Therefore, we only penalized against genes that did not perform well where they could be tested. A total of 56 genes are identified as the "high-priority" gene candidates based on this method, 35 genes from the loci identified in the absolute HNT or control and 21 genes from the relative analysis (Supplementary Table S4). Seventeen candidate genes were identified as strong candidates in two of the three evaluations, 12 in the absolute and 5 in the relative analysis. An additional 23 candidate genes were selected as priority candidates in absolute HNT or control based on one criterion but could not be tested in either of the other two methods. Similarly, 16 additional priority candidates identified by one method, but not testable in the other two were identified from Table 1. Summary of significant marker-trait associations (MTAs) with yield, 100-seed weight, harvest index, and chalkiness traits in experiment 1 (walk-in chamber study) and experiment 2 (field-based heat tent study) in rice. The list of significant MTAs and Manhattan plots are presented in Supplementary Table S1

Discussion
Increasing tolerance to HNT while simultaneously minimizing yield and quality losses is an emerging goal in rice breeding. A night temperature of > 27 ℃ has been used to determine differential physiological, biochemical, and agronomic responses both under controlled 19,53 and field experiments 11,13,32 , using a limited number of genotypes [two 11 to 36 33 ]. Here, we carried out stage-targeted phenotyping of indica diversity panel both under controlled-environmental and field conditions ( Supplementary Fig. S1), which allowed us to identify genetic loci associated with key phenotypic traits that determine HNT responses (Table 2).

Phenotypic variation explains underlying trait correlations and physiological responses under HNT.
Our study involving rice diversity panel (entries from different rice-growing regions) recorded a significant decline in yield, 100-seed weight, and increased chalkiness under prolonged HNT stress (Fig. 2), in both controlled environment and field studies, supporting previous findings 5,11,13 . A more significant reduction in grain yield and 100-seed weight, at least in part, is known to be related to reduced carbohydrate availability due to higher night-time respiration or reduced assimilates supply to sink or a combination of both, under HNT in rice 13 . HNT stress during grain development terminates or reduces the activity of a series of enzymes involved in starch synthesis in the endosperm, negatively affecting grain size under HNT in rice 13,35 and other crops 14 . For example, a significant correlation between changes in the fructose-6-phosphate content of the panicles and the yield reduction in rice emphasized the significant metabolic changes in the sink during grain development in response to HNT 35 . Thus, minimizing yield losses and quality deterioration caused by HNT during grain filling is becoming extremely important to sustain global rice production and end use quality. Unlike HDT stress at anthesis, post-flowering HNT reduces grain yield and quality by lowering grain weight caused by heat-induced premature senescence 14 and quicker break-down of stored assimilates during the night 13,32,35 . In particular, heat stress-accelerated chlorophyll degradation, inhibition of chlorophyll biosynthesis 54 ,  Table 2. List of candidate MTAs identified under control night (CNT) and high night temperature (HNT) in walk-in chambers (Experiment 1-E1) and field-based heat tents (Experiment 2-E2). Allelic effect with respect to minor allele (AE) = (average phenotypic value of accessions with the minor allele-average phenotypic value of accessions carrying the major allele). Annotated genes housed in ± 50 kb of the cMTAs/locus identified in the study. Rel-relative (high night temperature/control night temperature).
Scientific Reports | (2021) 11:6747 | https://doi.org/10.1038/s41598-021-85921-z www.nature.com/scientificreports/ and loss of carbon already fixed by daily photosynthesis 55 are key aspects that limit the source capacity during grain filling. The depletion of carbohydrates or imbalance with the source (decrease in photosynthesis and increase in respiration) during grain filling exhausts the assimilate reserves that are particularly needed for grain development under heat stress 14,34,56 . Similar to the impacts of HNT stress, low solar radiation also interferes with starch accumulation during grain filling, which significantly affects the seed weight and market value (quality) of rice 57 . An increase in average night-time VPD seems to be coupled with warmer night and lower relative humidity (Fig. 1b). Such an increase in night VPD has been shown to impact the plant productivity 58 . Hence, further studies are needed to disentangle the HNT vs. low solar radiation, differential VPD levels under HNT scenarios and their combined effect on yield and grain quality losses in rice. Interestingly, canopy functional stay-greenness is known to have the potential to enhance radiation use efficiency and photo-assimilates supply by improving the longevity of photosynthetic machinery during postanthesis. We hypothesize that the canopy's extended ability to generate additional assimilates under HNT could potentially fill the carbon gap created by increased carbon losses due to higher night respiration in response to HNT. Therefore, exploring the relevance of functional stay green in combination with an efficient transfer of photo-assimilates (water-soluble carbohydrates) from source to sink is proposed to minimize the HNT-induced yield and quality losses during grain filling in rice.

Comparison of GWAS revealed common loci for phenotypic traits across treatments. Grain
size/shape and chalkiness are two crucial components that determine rice grain quality 59 , with increased chalkiness having a significant negative impact on the economic returns 20 . Both grain size and chalkiness are polygenic quantitative traits 60 . The high genetic and phenotypic correlations among measured traits across treatments, particularly for seed weight and grain size, indicate strong potential surrogate traits for breeding HNT tolerance in rice. A weak correlation was observed between controlled environment and field for some traits, which might be due to one or more factors (thermal amplitude [Fig. 1b], low light condition, relative humidity) differing between these environments 23,30 . In spite of these differences, several GWAS loci were in the genic region or close to previously reported genes controlling grain shape traits, e.g., GS3 51,52 and GW5/qGSW5 49,50,61 across treatments. This indicates that (1) GWAS is an effective way to identify putative genes/target regions for HNT tolerance in rice and (2) exhibited the presence of strong and common genetic regulation across control and HNT in rice (Fig. 4). For example, a strong association for grain width was identified on chromosome 5 under both control and HNT conditions, which is in the region of a well-characterized gene, GW5 (Os05g0187500, encoding a calmodulin-binding motif family protein) that regulates grain width in rice 50 (Fig. 4c,d). Another prominent peak on chromosome 3 (Chr3pos16794259.1) associated with the gene, GS3 (Os03g0407400, a regulator of grain size and organ size), was detected for grain length across treatments (Fig. 4f,g). Co-localization of GWAS loci for grain size (length or width) variations under control and warmer nights suggest that tolerance to HNT stress for grain size is likely to have common regulators in rice (Fig. 4). Another cMTA was identified on chromosome 3 (Q28) for grain chalkiness, which is also near the GS3, indicating that the molecular mechanisms of grain chalkiness in rice might be rather complex than the grain size 59 . Therefore, further functional studies are needed to explore the association of grain size and chalkiness genes in regulating rice grain weight and quality under heat stress conditions. For instance, one genomic region on chromosome 2 (31.4-31.9 Mb) associated with Q2, Q3, Q8, and Q12 was shown to control yield and 100-seed weight either under control or HNT (Supplementary Fig. S2). Loci (Q3 and Q8) identified under control and HNT grain yield were close (25 kb) to the trehalose-6-phosphate phosphatase gene (TPP7, Os02g0753000), which is shown to display significant improvement in abiotic stress tolerance in   Table S3), which is found to enhance heat stress tolerance by modulating the reactive oxygen species (ROS) production, particularly H 2 O 2 and reproductive organ development under different abiotic stresses 65,66 . A similar family gene (OsANN10) has been mapped for the heat stress tolerance index in rice 67 . A defined role for these genes in combating heat stress at reproductive and grain-filling has not been fully explored.
Despite the variable number of genotypes and different experimental conditions, our phenotypic correlations and co-localized GWAS loci exhibited the presence of strong and common genetic regulation across treatments in rice. A genomic region associated with the relative grain length on chromosome 3 (33.47) co-localized with locus influencing absolute spikelet fertility under high temperature 68 and relative spikelet fertility 69 . Two cMTAs for absolute grain yield under HNT on chromosome 1 (Q5, Chr1pos22255654.1) and chromosome 2 (Q6, Chr-2pos34468097.1) were related to genomic regions for relative grain yield in experiment 1 (Table 2). Similarly, a common MTA was observed between harvest index under HNT (Q25, Chr11pos3631830.1) and relative harvest index in experiment 1 (Table 2). Therefore, further characterization of loci (31.46 to, 34.46 Mb on chromosome 2 and 33.47 on chromosome 3) would be interesting to decipher the role of these pleiotropy regions for enhancing heat stress resilience (both day-time and night-time) in rice. Relatively little is known about the interaction between high day and night temperature stress on grain yield and quality 70 . Identifying genomic regions associated with combined stress tolerance would set a stage to breed rice with improved heat stress tolerance.
A limited number of shared loci between relative and absolute phenotype were observed across treatments suggesting that mapping genetic loci based on stress index (stress effects relative to control) may not always reflect the true stress tolerance 68 . While comparing the relative change in phenotypes between control and HNT allows the evaluation of sensitivity to HNT stress across all phenotypic ranges. In addition, evaluating the candidate genes in these loci may provide firsthand targets for pathways that alter adaptation to HNT stress. Nevertheless, quantifying stress tolerance by combining both absolute and relative traits would help develop rice accessions with greater productivity and adaptability under stress. The identified loci can be used as a starting point for future candidate gene validation and breeder-friendly markers development to enhance HNT tolerance in rice.

Refinement of candidate genes by available molecular data.
Although the possible genetic determinants of phenotypic trait variations in response to stress have been effectively determined based on GWAS in crops, detecting their functional relevance or prioritizing candidate genes is a severe bottleneck 46 . Using available prior data, we refined the possible candidate genes associated with the traits for each genomic region through three independent methods: evaluating the transcriptional response of the genes in each locus to HNT, determining the connectedness of each gene to transcriptional regulators in a gene regulatory network, and determining the orthology to genes previously identified as causative in other studies. Validation of the causal gene for each GWAS-identified locus is a key component to move forward. However, it is not practical to experimentally test all the genes that overlay the genomic regions. Therefore, there is a need to prioritize the target genes. Previous genetic mapping or GWAS studies have routinely used currently available literature to decide on the potential candidates to validate based on annotations and known functions. This approach, although widely used, introduces bias at different levels 45,46 . Hence, developing new approaches that reduce this bias would be beneficial. The three methods used to find top candidates from this study demonstrate integrating multiple ways to prioritize candidates for further validation. One limitation of this approach is that not all candidate genes could be tested for many of these datasets. For example, genes not on the Agilent rice 4 × 44 K microarray were not included in the gene regulatory network. Likewise, genes not expressed in the tissue evaluated for HNT response were not included in the transcriptional analysis. However, as data resources extend and more data becomes available in accessible formats, these approaches will become even more reliable for narrowing down the candidate genes to test for causality.  Table 3. List of genes closely related to absolute traits under HNT and relative across experiments. Genes located within the vicinity of significant GWAS sites were used to identify closely connected genes using the Network analysis. A list of annotated genes around the cMTAs is given in the Supplementary Table S2. E1-walk-in chambers (Experiment 1) and E-2-field-based heat tents (Experiment 2). The In_DE_Set, Rank, and In_Network columns correspond to the different gene candidate prioritization methods. These three columns indicate the gene's presence in the DE the gene's rank in its respective loci and network datasets, respectively. YPP grain yield, HWT 100-seed weight, HI Harvest index, CHALK chalkiness, GL grain length, GW grain width.

Conclusion
The current study is the first report investigating the natural variation of rice genotypes to HNT using walk-in chambers and field-based heat tents. The observed significant positive correlations (Fig. 3) and common GWAS loci across treatments (Fig. 4) for grain size indicate that selection for higher phenotypic values from the control conditions could be a proxy for high yields under HNT. Advanced gene sorting methods combined with GWAS allowed us to narrow down the candidate region to potential genes controlling yield and related traits across control and HNT. Exploring the genetic regulation of functional stay green and source-sink carbohydrate relationships under night temperature is an exciting area for further research. Our findings provide insights towards integrating multi-disciplinary approaches such as physiological, genomic, and network algorithms to identify possible causal loci/genes. The identified favorable alleles for HNT tolerance would help develop molecular markers to support ongoing heat stress breeding programs in rice.