The rice seeds of ZH11 and its OsGRF4^AA near-isogenic line, GS2, were provided by Dr. Jiang Hu from the State Key Laboratory of Rice Biology, China National Rice Research Institute (Hu et al., 2015). The seeds were first grown in the China National Rice Research Institute (Hangzhou, China, 119°57′E, 30°03′N) under natural conditions until the pollen mother cell meiotic stage (Feng et al., 2018). Thereafter, the rice plants were divided into two groups and moved into separate transparent growth chambers with a relative humidity of around 70%–80% and natural sunlight conditions. Group 1 was subjected to a HS condition (i.e., 28°C from 16:01 to 8:59, and 38°C from 9:00 to 16:00), whereas the other group was grown under a normal temperature condition (i.e., 23°C from 16:01 to 8:59 and 28°C from 9:00 to 16:00). The treatment time was 6 days. Anthers were sampled at 16:00 on the last day of the treatment. Sampled anthers were immediately submerged in liquid nitrogen for 1 min and stored at –80°C until RNA extraction and sequencing library preparation.

Pollen viability was determined using the method described by Feng et al. (2018). Briefly, mature pollen grains were removed from the anthers of the florets, and placed on a glass slide with a drop of potassium iodide/iodine (KI/I₂) solution. The slide was observed using a light microscope (DM4000B; Leica, Wetzlar, Germany). Eight replicates were investigated and photographed for the measurement of pollen viability.

Mature rice plants and panicles were photographed for phenotyping. Then the number of filled grains (FG) and abortive grains (AG) per panicle were measured. The seed-setting rate was calculated as FG/(FG + AG) × 100%.

RNA sequencing (RNA-Seq) and primary data analysis were performed at Shanghai Majorbio Bio-pharm Biotechnology Co., Ltd. (Shanghai, China). Total RNA was isolated from rice anthers using a plant RNA purification reagent (Invitrogen, California, USA). The quality and integrity of the RNA were measured by an Agilent Bioanalyzer 2100 system (Agilent Technologies, USA) and only high-quality RNA samples [OD_260/280 = 1.8–2.2, OD_260/230 ≥ 2.0, RNA integrity number (RIN) ≥ = 8.0, 28S : 18S ≥ 1.0, > 1 μg] were used to construct the sequencing library. RNA-Seq libraries were sequenced on an Illumina NovaSeq6000 platform to generate high-quality paired-end reads following the Illumina manufacturer’s recommendations. Next, data processing for differentially expressed gene (DEG) analysis was carried out using Kallisto (Bray et al., 2016) and 3D RNA-Seq App (Guo et al., 2021), with default settings (adjusted p-value < 0.01 and absolute log₂-fold change > 1) unless specifically indicated otherwise. Pre-mRNA splicing analysis was performed using rMATS (replicate multivariate analysis of transcript splicing) (Shen et al., 2014). RNA read coverages were visualized using the program Integrative Genomics Viewer (IGV) (http://www.igv.org) by Rice Annotation Project (RAP) International Rice Genome Sequencing Project (IRGSP) 1.0 genome (https://rapdb.dna.affrc.go.jp/) as the reference genome.

Reverse transcription-PCR (RT-PCR) and reverse transcription-quantitative PCR (qRT-PCR) were slightly modified compared with the methods previously described in (Chen et al., 2022). DNA digestion of total RNA samples was performed after RNA extraction by using an RNase-Free DNase Set (Vazyme, Catalog No. R223-01). Total RNA was then reverse transcribed using a SuperScript First-Strand Synthesis System to generate the first-strand cDNA libraries. qPCR was performed using Power SYBR Green PCR Master Mix (Yeasen Biotech; Catalog No. 11201ES08) under the following conditions: 95°C for 5 min, then cycles of 95°C for 15 s and 60°C for 1 min. The validation was assessed with three biological replicates for each sample. An actin gene (Os03g0718150) was used as the internal control for normalization of qPCR. The fold change was evaluated using the standard 2^(−ΔΔCT) method. Statistical analysis of gene expression in qRT-PCR was performed using Microsoft Excel 2013. In RT-PCR validation of AS events, primers flanking different exons were used. The PCR program was as follows: initial denaturation for 3 min at 95°C, followed by 40 cycles at 95°C for 30 s, 56°C for 30 s, and 72°C for 2 min. The final extension at 72°C lasted 5 min.

Pollen grains from ZH11 and GS2 with and without HS treatment were observed and measured using a microscope after KI/I₂ staining ( Figure 1 ). The viability of pollen grains from both genotypes grown under control conditions was high. No significant difference was found between ZH11 and GS2. Dead pollen grains accounted for only 5.4% and 3.8% of pollen grains in ZH11 and GS2, respectively. The viability of pollen grains was decreased significantly by HS treatment in both genotypes. However, more severe heat-induced damage of pollen grains was found in GS2 anthers. The percentage of dead pollen grains was 34.6% in ZH11 after HS treatment, whereas this proportion increased significantly, to 62.6%, in GS2.

Next, we calculated the seed-setting rates of both genotypes. The seed-setting rates of ZH11 and GS2 under control conditions were 82.3% and 83.3%, respectively. The seed-setting rate was decreased significantly by HS in both genotypes. In addition, the rate of decrease was greater in GS2 than in ZH11, as we found that fertilized grains accounted for 35.4% of all grains in ZH11 after HS treatment, whereas only 21.1% of grains were fertilized in GS2 after HS treatment.

By running a RNA-Seq experiment with the collected anthers, we obtained 40–45 million raw reads from each of the 12 cDNA libraries, and at least 95.6% of raw reads from each library could be assembled to the reference sequence of the rice genome ( Supplementary Table 1 ). In the 3D RNA-Seq analysis pipeline, low-count reads, that is, reads lower than one per million, and reads found only in one library, were removed from further analysis ( Figures 2A, B ). In principal component analysis, different samples separated from each other, while all three replicates from the same sample grouped together, suggesting high similarity between replicates ( Figure 2C ). The large distances between the ZH11 control group and other samples indicate that their transcriptomic profiles are very different. In contrast, the ZH11 and GS2 HS groups were close to each other, suggesting less transcriptomic difference between these two samples. Correspondingly, the analysis of DEG numbers showed that the numbers of down-regulated and up-regulated DEGs were, respectively, 2,112 and 816 in the comparison of the ZH11 and GS2 control groups; 3,256 and 1,976 in the comparison of the ZH11 control group and the ZH11 HS group; 1,476 and 1,429 in the comparison of the GS2 control group and the GS2 HS group; and 451 and 507 in the comparison of the ZH11 and GS2 HS groups ( Figure 2D ). Furthermore, the Venn diagram demonstrated that there were considerable numbers of overlapping and specific DEGs in different comparisons ( Figure 2E ).

Further analysis revealed a group of 961 genes that were up-regulated by HS in both ZH11 and GS2, and another group, of 820 genes, which were down-regulated in both varieties by HS ( Supplementary Tables 2, 3 ). Gene ontology (GO) analysis revealed that HS up-regulated genes were enriched in the biological processes “response to heat”, “protein folding”, “response to hydrogen peroxide”, “response to high light intensity”, and so on ( Figure 3A ). On the other hand, HS-inhibited genes tend to be enriched in biological processes involved in cell differentiation and cell structure control, flower development, and gene transcription ( Figure 3B ). Inhibition of such genes explained why we observed a significant reduction in the viability of pollen grains from both genotypes.

Next, we employed two extra sets of RNA-Seq data,which have also been used to explore HS-induced transcriptomic changes in rice anthers (Zhang X. et al., 2012; Liu et al., 2020). Our analysis demonstrated that 103 and 84 genes from the commonly up-regulated DEGs in ZH11 and GS2, respectively, overlapped with one of the two RNA-Seq data sets employed, whereas another group of 55 DEGs overlapped with both data sets. In the case of HS-inhibited genes, 100, 138, and 32 genes from commonly down-regulated DEGs in ZH11 and GS2 were also found to be repressed in at least one of the employed data sets ( Figure 3C ). Further study showed that 40% of 55 commonly up-regulated genes in the three studies were enriched in the GO term of “response to heat”. This group of genes comprised, among others, HSFA2a, Heat Shock Protein 101 (HSP101), Heat Stress-associated 32kd protein (HSA32), and a group of sHSPs. In addition, a gene encoding a splicing factor, RNA-binding Protein 1 (RNP1) (Os01g031660), was found to be commonly up-regulated by HS in anthers from different rice genotypes. In contrast, another group of 32 genes were commonly down-regulated by HS treatment. These commonly suppressed genes included genes involved in cell division and flower development, such as an aquaporin protein (Os01g0232000), two GDSL-like lipase/acylhydrolases (Os02g0816200 and Os06g0156600), a lipid transfer protein (LTPL65, Os01g0814100), a transcription factor (Os02g0791300), a hormone-responsive gene (Os11g0128300), a core histone H2A (Os05g0113900), and a chloroplast precursor (Os04g0678700) (Figure 3D). These results suggest that these 87 genes may play important roles in the HS response of anthers from different rice genotypes.

As it has been found that the number of dead pollen grains in anthers after high-temperature treatment is greater in GS2 than in ZH11, we next aimed to uncover the transcriptomic difference between anthers from the two genotypes after HS. A heatmap comparing the patterns of expression of 958 DEGs in GS2 and ZH11 after HS treatment is shown in Figure 4A . The heatmap shows that the heat-mediated regulatory patterns of many genes were similar in the two genotypes. However, the same genes were expressed at different levels after HS treatment, resulting from either originally different levels under normal temperature conditions or different intensities of regulation between the two genotypes by HS treatment, or both. Interestingly, HS-mediated gene expression changes seemed to be more marked in GS2. After HS treatment, 80 HS-induced DEGs had higher levels of expression in GS2 than in ZH1, whereas 121 heat-suppressed DEGs were, in contrast, expressed at lower levels in GS2 than in ZH11 after HS treatment; one gene from this category was expressed at lower level in ZH11 than in GS2 after HS treatment ( Figure 4B ).

The change trends of genes detected from RNA-Seq data could be validated through qRT-PCR. Genes whose regulation under HS showed the same trend in ZH11 and GS2 included a B‐box containing protein BBX19, HSFA2a, MALE STERILITY 2 (MS2), and ABA 8′-Hydroxylase (ABA8ox2) ( Figure 4C ). Regulation of some other genes under HS showed opposing trends in ZH11 and GS2. For example, Dehydration Responsive Element Binding Protein 1A (DREB1A) and Peroxidase 1 (Pox1, Os05g0499300) were induced in GS2, but inhibited in ZH11, after HS treatment.

Under normal conditions, a group of genes with higher expression levels in ZH11 than GS2 was enriched in GO terms including “carbohydrate metabolic process”, “lipid metabolic process”, “lipid transport”, and “photosynthesis”. In contrast, genes expressed at higher levels in GS2 were highly enriched in biological processes involved in pollen tube germination and growth, pollination, cell growth, transmembrane transport, and negative regulation of catalytic activity ( Supplementary Figure 1 ). However, we did not find a significant difference in pollen viability under control conditions.

Further analysis revealed that DEGs with higher levels of expression in ZH11 after HS treatment are involved in tens of biological processes ( Figure 5A ). About 60 DEGs are thought to play roles in response to abiotic stress, and 40 genes are involved in the response to endogenous stimuli. Interestingly, a group of DEGs with significantly higher levels of expression in ZH11 than in GS2 after HS treatment play important roles in chloroplast development and photosynthesis. In contrast, further GO analysis demonstrated that DEGs with lower levels of expression in ZH11 than in GS2 after HS treatment were greatly enriched, and thus with relatively low numbers, in five biological processes ( Figure 5B ), implying that the biological processes these genes are involved in have not yet been identified.

In fact, a great number of DEGs involved in chloroplast and photosynthesis are also responsive to abiotic stress ( Figures 5C, D ). Only a few abiotic stress-responsive genes, such as HSFA2a and OsBiP5, were up-regulated by HS, whereas most of the remaining genes were down-regulated in both genotypes by HS treatment ( Figure 5C ). Similarly, many genes related to photosynthesis were also down-regulated in both genotypes. However, these photosynthesis genes were down-regulated to different levels in ZH11 and GS2 by the same HS treatment. All genes, except CYP19-4, were expressed at higher levels in ZH11 than in GS2 after HS treatment ( Figure 5D ). Given that the abundance of dead pollen grains was greater in GS2 than in ZH11 after HS, we supposed that sustaining the expression of these genes at relatively high levels might contribute to reduced heat-induced damage to rice anthers.

Next, we investigated a group of genes that have been identified as male sterile genes (MSGs) in rice, the mutation of which would cause total or partial male sterilization in rice [as reviewed by Abbas et al. (2021)]. Interestingly, 8 out of 77 MSGs were expressed at higher levels in ZH11 than in GS2 after HS treatment (Figure 5E). In contrast, no MSGs were expressed at significantly higher levels in GS2 than in ZH11 after the same HS treatment. The eight MSGs with higher levels of expression in ZH11 than in GS2 after HS were five genes involved in lipid metabolism [i.e., OsC6 (Os11g0582500), OsABCG26 (Os10g0494300), OsDPW (Os03g0167600), OsSTRL2 (Os03g0263600), and OsABCG3 (Os01g0836600], a bHLH transcription factor [TIP2/TDR/bHLH142 (Os01g0293100)], an R2R3 MYB transcription factor Carbon Starved Anther (CSA, encoded by Os01g0274800), and a triterpene synthase [OsC12/PTS1 (Os08g0223900)] (Figure 5F). OsC6 was also one of the top 50 DEGs with much higher expression levels in ZH11 than in GS2 after HS treatment ( Supplementary Figure 2 ).

As the difference between the two rice genotypes used here was thought to have resulted from a couple of nucleotide variations (1,187TC→AA) in a GRF, OsGRF4, we were curious to know the patterns of expression of this gene and of its family members in the anthers of the two genotypes. As shown in Figure 6A , among these 12 GRF genes, OsGRF11 showed the highest level of expression in anthers, followed by OsGRF8, under either normal or high temperature. Almost all GRFs were suppressed by HS in both genotypes, except for GRF10 and GRF11: the former was suppressed by HS in GS2 but induced in ZH11, whereas the latter was induced by HS in GS2 but suppressed in ZH11. The expression level of OsGRF4 in anthers was lower than that of many other GRF genes. Furthermore, the transcripts of this gene seemed to accumulate less in ZH11 under HS conditions when tested by qPCR ( Figure 6B ).

When observing read coverages of OsGRF4 using the IGV program, we found another sequence variation, in addition to the previously reported variation (1,187TC→AA). This novel variation is located at 3′-UTR of the last exon, where an ‘A’ was replaced by a ‘T’ in GS2 (as this variation was found in all GS2 genotype samples but not in any ZH11 genotypes, we suggest that it is not a random error caused by sequencing). More interestingly, we found that the pre-mRNA splicing patterns of OsGRF4 were different in ZH11 and GS2 genotypes after HS, when read coverages of the target genes were visualized using the IGV program ( Figure 6C ). The last intron of this gene underwent intron retention (IR) in ZH11 after HS. This AS regulation was not found in any of the GS2 samples, either with or without HS treatment. In other words, the splicing patterns of the gene in the GS2 background remained the same as under the normal temperature condition. We further tested this AS event using RT-PCR. We found that an extra DNA fragment (252 bp) was amplified in ZH11 HS group samples ( Figure 6D ). This novel amplified fragment was predicted to be the intron-retained isoform owing to its size. Moreover, the AS of OsGRF4 seemed to be carried out only in ZH11 anthers after HS treatment, as we could not detect the same AS event in leaves from either genotype under the same temperature conditions ( Figure 6E ). Bioinformatic analysis revealed that the resulting mRNA isoform of OsGRF4 in ZH11 anthers induced a stop codon in the retained intron, which would induce RNA degradation if it triggered the nonsense-mediated degradation (NMD) system. Otherwise, the resulting mRNA isoform was predicted to be translated into a sequence-varied protein, with the last amino acid, D (Asp), at the C-terminal replaced by three contiguous T (Thr) residues ( Figure 6F ). Thus, the mechanism for the appearance of AS in OsGRF4 is yet not clear, and determination of its biological significance for the heat tolerance of anthers and pollen grains in rice requires further observation.

Next, we employed a widely used pre-mRNA splicing analysis program, rMATS, to explore AS events in the RNA-Seq samples. As shown in Figure 7A , AS events were defined as a skip exon (SE), a retained intron (RI), mutually exclusive exons (MXE), an alternative 3′ splice site (A3’SS), and an alternative 5′ splice site (A5’SS) in our analysis. There were 2,838, 2,442, 1,673, and 1,205 significant differential AS (DAS) events in comparisons of ZH11 control group vs. ZH11 HS group, GS2 control group vs. GS2 HS group, ZH11 control group vs. GS2 control group, and ZH11 HS group vs. GS2 HS group, respectively. SE and A3′SS events accounted for more than two-thirds of the total AS events in each comparison, with IR and A5′SS events the next most frequent. Interestingly, only 23 genes were considered to be DEGs and also DAS genes across ZH11 and GS2 after HS, that is, less than 3% of these 1,960 disturbed genes were regulated together at transcriptional and splicing levels ( Figure 7B ). GO analysis revealed that DAS genes tend to be highly enriched in nucleic acid-binding processes, such as DNA metabolism, RNA processing, gene transcription, etc. A group of genes involved in pollen grain development could also be differentially spliced. As we found that DAS genes were also enriched in GO terms, such as “starch biosynthetic process”, “mitotic cell cycle”, and “floral organ formation” ( Figure 7C ).

We randomly picked up several genes that was considered to be differently spliced in ZH11 and GS2 under heat stress conditions, and validated them through RT-PCR. It was demonstrated that encoding genes of important regulators, such as CCA1, OsSPO11, HSFB2c, and OsDR11, could amplify more than one target band ( Figures 7D–J ), suggesting that different mRNA isoforms could be generated after the pre-mRNA splicing process. Moreover, the intensities of amplified fragments of the same genes varied between ZH11 and GS2 after HS, suggesting different abundances of alternatively spliced mRNA isoforms in these two samples, which corresponds to snapshots from the IGV program (shown below each RT-PCR gel panel). Notably, not all AS events detected through RT-PCR and the IGV program were considered significant DAS events in the rMATS program. Reasons for this gap could be the relatively low read coverage of AS events in RNA-Seq data or the technical preference of the data analysis program (i.e., rMATS focusing on variations of junction counts between different exons, but not IR (Shen et al., 2014).