A panel of 620 rice accessions from the 3K Rice Genome Project (3K RG) (Wang et al., 2018) was used to evaluate heat tolerance at the seedling stage. The accessions contained 173 Geng, 411 Xian, 19 admix, 7 Aus, 9 Basmati and 1 unknown accessions ( Supplementary Table S1 ). Twenty-four uniformly germinated seeds per replicate of each accession were sown in 96-well plates with holes at the bottom of each well. Then, the seeds were soaked in a container with tap water by placing the plates on scaffolds and were cultured in a phytotron at 28°C/25°C, 70% relative humidity and a 13-h light/11-h dark photoperiod. After 7 d, the seeds were transferred to Yoshida solution (pH 5.8-6.0), which was replaced every 3 d (Li et al., 2015). 13-day-old seedlings were exposed to 45°C for 3 d in a phytotron and then returned to normal conditions (28°C) for 7 d of recovery. The phytotron was set at 60% relative humidity and low light intensity (50–80 µM m⁻² s⁻¹) to minimize the influence of high light and hydrophobic stress (Hasanuzzaman et al., 2013). Then, the survival rate (SR) was calculated as the proportion of surviving seedlings. Based on the evaluation system ( Table 1 and Supplementary Figure S1 ), the leaf score of heat tolerance (SHT) was determined by visual inspection. At least three biological replicates were performed.

A total of 2,802,578 SNPs with a missing rate < 0.1 and minor allele frequency (MAF) ≥ 0.05 in the GWAS panel were filtered from the 3K-RG 4.8M SNP dataset (Alexandrov et al., 2014) by PLINK (Purcell et al., 2007). The GWAS based on a mixed linear model was performed with EMMAX (Kang et al., 2010) to identify the associations between SNPs and heat tolerance. The kinship matrix was calculated with an identical-by-state matrix using the pruned SNP subset (with the parameter “indep-pairwise 50 10 0.1” in PLINK) as a measure of relatedness between accessions. The eigenvectors of the kinship matrix were calculated using GCTA (with the parameter “-make-grm”) (Yang et al., 2011) and the first three principal components were used as covariates to control population structure. The effective number of SNPs (N) was calculated by the GEC software (Li et al., 2012), and a suggestive significance threshold of association (P = 2.29E-06) was determined by the Bonferroni correction method (1/N) for claiming significant SNPs. Manhattan plots of the GWAS results were plotted by the R package “qqman” (Turner, 2014). The significant SNPs within the 300-kb region were considered as a locus based on the previously reported linkage disequilibrium (LD) decay in 3K RG (Wang et al., 2018). The leading SNP within a locus was defined as the SNP with the lowest P value. Local LD block analysis was performed within 150 kb upstream and downstream of the leading SNP using the LDBlockShow (Dong et al., 2021).

The haplotype analysis was performed on each annotated gene in qHT7 to identify candidate genes and unearth favorable haplotypes. The gene haplotypes were constructed with all SNPs in the coding sequence (CDS) and 1-kb promoter regions, respectively. The synonymous SNPs were merged into one haplotype following the method by Zhang et al. (Zhang et al., 2021). Duncan’s multiple range post-hoc tests were used to compare phenotypic differences between haplotypes (n ≥ 40 rice accessions). The module of Custom Genotyping and Comment (Rice) in MBKbase database (http://www.mbkbase.org/rice/customGT ) (Peng et al., 2019) was used to construct the candidate gene’s variety groups based on the SNP genotype of the published wild rice accessions and 3K RG with parameters “Sample Num: ≥ 40, ALT ≥ 5%, Missing ≤ 20%”. The haplotype network of a candidate gene was drawn by the minimum-spanning tree in Popart (Leigh and Bryant, 2015).

One representative heat-tolerant Xian accession, FACAGRO 64 (F64), and a representative heat-sensitive Xian accession, PUILLIPINA KATARI (PK), were selected for transcriptome analysis. Shoot samples before and after 24 h of heat-stress treatment were collected and stored in liquid nitrogen, each with three biological replicates. Total RNA was extracted from shoot samples using the TRIzol reagent (Invitrogen) and then treated with RNase-free DNase I (Takara) to remove genomic DNA. Sequencing libraries were constructed according to the standard protocols provided by Illumina. The libraries were sequenced using Illumina NovaSeq 6000platform (150-bp paired ends) in Novogene (China). The raw sequence data reported have been deposited in the Genome Sequence Archive (Chen et al., 2021b) in National Genomics Data Center (CNCB-NGDC Members and Partners, 2022) with accession number CRA008760 that are publicly accessible at https://ngdc.cncb.ac.cn/gsa.

After removing adaptor and low-quality reads, clean reads were aligned to the Nipponbare reference genome (MSU v7.0) using HISAT2 (Kim et al., 2015). The gene expression levels based on fragments per kilobase of exon per million mapped fragments (FPKM) were calculated by StringTie (Pertea et al., 2015). Differentially expressed genes (DEGs) between two samples were identified with the DESeq2 package (Love et al., 2014) in R. The threshold for claiming DEGs was set as adjusted P-value (FDR) ≤ 0.05 and log fold change (FC) absolute value ≥ 1. Functional enrichment analysis of Gene Ontology (GO) and KEGG pathway was performed by clusterProfiler software (Yu et al., 2012). The threshold of adjusted P-value (FDR) < 0.05 was used to identify significantly enriched GO terms and KEGG pathways.

Total RNA (1 ug) was reverse-transcribed into cDNA using FastKing gDNA Dispelling RT SuperMix kit (Tiangen; KR118-02). qRT-PCR analyses were performed with SuperReal PreMix Plus (SYBR Green) kit (Tiangen; FP205-2), including three biological replicates. UBQ was used as the internal control, and the relative expression levels of the target genes were quantified using the comparative cycle threshold (2^-ΔΔCT) method (Livak and Schmittgen, 2001). Primers used for qRT-PCR are listed in Supplementary Table S2 .

The shoots of F64 and PK under heat stress (45°C) for 72 h and control conditions were collected, respectively. The levels of superoxide dismutase (SOD), malondialdehyde (MDA), peroxidase (POD) and hydrogen peroxide (H₂O₂) in shoot tissue were determined using SOD, MDA, POD, H₂O₂ commercial kits following the manufacturer’s instructions (Suzhou Grace Biotechnology Co., Ltd.). Three biological replicates were included.

The phenotypic measurements of SR and SHT showed a considerable variation in heat tolerance at the seedling stage among the 620 rice accessions ( Figure 1 and Supplementary Table S1 ). The mean SR in the whole population was 46.9%, ranging from 0 to 100.0%. Similarly, the mean SHT was 6.82, with a range of 1.29 to 9.00. Heat tolerance at the seedling stage did not vary significantly between Xian and Geng subpopulations ( Figures 1A, C ). Among the four Geng subgroups, most accessions with high heat tolerance belonged to GJ-trp ( Figures 1B, D ). Within Xian subpopulation, the average heat tolerance of XI-3 accessions was higher than those of other Xian subgroups. The heat tolerance of the GJ-trp accessions was similar to that of the XI-3 subgroups. Moreover, most GJ-trp and XI-3 accessions are both from Southeast Asia islands.

Thirty-one SNPs significantly associated with heat tolerance were identified in the 620 accessions, including 29 and 2 SNPs associated with SR and SHT, respectively ( Figures 1E, F and Supplementary Table S3 ). Among the significant SNPs, 6, 3 and 22 were located in the promoter, CDS and intergenic regions of 17 annotated genes, respectively. Based on the local LD block analysis, we combined these significant SNPs into six loci distributed on rice chromosomes 1, 3, 6, and 7 ( Figure 1G ). By comparing the 58 previously reported QTLs for heat tolerance (Xu et al., 2021) and 23 known genes involved in heat tolerance (Huang et al., 2022), three genes/QTLs, qHTB1 (Zhu et al., 2017), qHTB3-3 (Jagadish et al., 2010; Zhu et al., 2017; Kilasi et al., 2018), and TT2 (Kan et al., 2022) were also found in the region of qHT1, qHT3.2, and qHT3.3, respectively. Out of the six loci, qHT7 (Chr7: 29067638-29223510 bp) was determined as the major locus since it contained the most and strongest association signals ( Figure 1G ).

F64 was highly tolerant to heat stress, with an average 95.8% ± 7.2% SR, which only slightly dried at the tips after 7 d of recovery from heat stress ( Figures 2A, B ). PK was extremely sensitive to heat stress, with an average SR of 0% ± 0%, and all of the seedlings were apparently dead. To examine the cell membrane damage and redox homeostasis caused by heat stress in the two rice accessions, we compared the physiological traits between F64 and PK under heat stress for 72 h ( Figures 2C–F ). Although there was no statistically significant difference in the relative MDA content between the two accessions ( Figure 2F ), the relative H₂O₂ content of F64 after 72 h of heat stress was significantly lower than that of PK ( Figure 2E ). Moreover, the relative activity of the antioxidant enzymes SOD and POD were significantly higher in F64 than in PK ( Figures 2C, D ). These results suggested that F64 suffered less damage to cell membranes under heat stress than PK, possibly due to more effective active detoxification by ROS scavenging regulation in F64.

In order to reveal the differences in transcriptome response to heat stress at the seedling stage between rice accessions with different levels of heat tolerance, we compared F64 (a representative heat-tolerant accession) with PK (a representative heat-sensitive accession) using RNA-seq analysis ( Supplementary Table S4 ). A total of 2056, 8303, 4070 and 8717 DEGs were identified for G1 (F64 vs PK under control conditions), G2 (heat stress vs control in F64), G3 (F64 vs PK under heat stress) and G4 (heat stress vs control in PK), respectively. Among them, 1202, 4287, 2311 and 4365 DEGs were upregulated, and 854, 4016, 1759, and 4352 DEGs were down-regulated in G1, G2, G3 and G4, respectively ( Figures 3A, B ). A series of biological processes and pathways involved in response to heat stress were commonly identified in both heat-tolerant and heat-sensitive accessions. Specifically, gene ontology (GO) analysis for the G2 and G4 DEGs, which were significantly regulated by high temperature in heat-tolerant and heat-sensitive accessions, respectively, showed that the common biological processes were mainly upregulated in protein folding (GO: 0006457) and RNA processing (GO: 0006396) ( Figure 3C ), and were primarily downregulated in carbohydrate metabolic process (GO: 0005975), biosynthetic process (GO: 0009058), metal ion transport (GO: 0030001), and glycolytic process (GO: 0006096) ( Figure 3E ). Similarly, six KEGG pathways, including spliceosome (map03040), protein processing in endoplasmic reticulum (map04141), RNA degradation (map03018), RNA transport (map03013), ribosome biogenesis in eukaryotes (map03008), and valine, leucine and isoleucine degradation (map00280), were significantly enriched both in the G2 and G4 upregulated DEGs ( Figure 3D ). For the G2 and G4 downregulated DEGs, 12 common KEGG pathways were significantly enriched, such as carbon metabolism (map01200), biosynthesis of amino acids (map01230), etc ( Figure 3F ). The results suggest that the aforementioned biological processes and pathways mentioned above should be the components of regulatory mechanisms underlying heat tolerance in rice.

Moreover, several unique GO terms and KEGG pathways were identified in G2 DEGs compared to the G4 DEGs. Four specific GO terms and one KEGG pathway were significantly enriched in G2 upregulated DEGs compared to G4 upregulated DEGs, including cell redox homeostasis (GO: 0045454), ribosome biogenesis (GO: 0042254), carbohydrate metabolic process (GO: 0005975), protein metabolic process (GO: 0019538), and RNA polymerase (map03020). In contrast, the divergence between G2 and G4 downregulated DEGs was much greater according to the number of unique GO terms and KEGG pathways ( Figures 3E, F ). Interestingly, cell redox homeostasis (GO: 0045454) was specifically enriched in G4 downregulated DEGs compared to G2 downregulated DEGs.

Genes involved in cell redox homeostasis are usually triggered in plants tolerant to abiotic stresses (Awasthi et al., 2015). Given that the genes related to cell redox homeostasis exhibited contrasting responses to heat stress in the two rice accessions with different heat tolerance, we further compared the expression profiles of 62 DEGs related to cell redox homeostasis between F64 and PK under control and heat stress conditions ( Supplementary Figure S2 ). There were 34 (55%) common, 12 (19%) F64-specific and 10 (16%) PK-specific DEGs regulated by heat stress in the two accessions. The results suggest that cell redox homeostasis should play an important role in rice heat tolerance.

Based on the Nipponbare reference genome IRGSP 1.0, 28 genes were annotated in qHT7 ( Figure 4A and Supplementary Table S5 ). Candidate genes for heat tolerance were selected based on the following criteria: (1) functionally related to abiotic stresses based on the annotation of Nipponbare reference genome, GO annotation, and literature search; (2) significant differences in heat tolerance among gene haplotypes. Consequently, 14 candidate genes were identified ( Figure 4B and Supplementary Table S5 ). To further screen the promising candidate genes, we examined the expression profiles of the 14 candidate genes using the transcriptomic datasets of F64 and PK. As a result, five DEGs (LOC_Os07g48830, LOC_Os07g48630, LOC_Os07g48710, LOC_Os07g48570, and LOC_Os07g48760) were selected ( Figure 4B ). We also verified the expression of the five genes by qRT-PCR ( Figure 4E and Supplementary Figure S3 ). Among the five genes, only the expression level of LOC_Os07g48710 was significantly higher in the heat-tolerant accession F64 than that in the heat-sensitive accession PK under heat stress ( Figure 4E ), which was consistent between the RNA-seq and qRT-PCR results. Thus, LOC_Os07g48710, encoding a VQ domain-containing protein, was determined as a promising candidate gene for the further analysis.

Mining heat-tolerant allele is helpful in improving the heat tolerance of rice through molecular breeding. To examine the favorable haplotype of the promising candidate gene of qHT7, LOC_Os07g48710, we performed the haplotype analysis using CDS and 1-kb promoter SNPs in 3K RG. Due to no SNPs detected in the CDS of LOC_Os07g48710, we identified four major haplotypes (n ≥ 40 accessions) using 20 SNPs (MAF ≥ 0.05 and heterozygous rate < 0.05) in its 1-kb promoter region in the GWAS panel ( Figure 4C ). Among the four haplotypes, Hap4 with significantly higher SR was determined as the favorable haplotype ( Figure 4D ), which was significantly enriched (P = 2.63E-54) in GJ-trp accessions of the GWAS panel ( Figure 4F ). For the 3K RG, the heat-tolerant haplotype LOC_Os07g48710 ^Hap4 was also highly enriched in the GJ-trp accessions (P = 1.60E-287) ( Figure 4G ). In contrast, LOC_Os07g48710 ^Hap4 was virtually absent in GJ-tmp subpopulation and Xian subpopulation ( Figures 4F, G ). To explore the origin and spread of Hap4, the haplotype network of LOC_Os07g48710 was analyzed, showing that Hap4 possibly evolved from Hap1 ( Figure 4I ). Furthermore, the proportion of Geng accessions with LOC_Os07g48710 ^Hap4 dropped dramatically from 60.3% in landrace to 33.2% in modern variety ( Figure 4H ).