This experiment was conducted in a screenhouse equipped with two heat-treatment chambers at The University of Western Australia’s Shenton Park Field Station, Perth, Western Australia (31°56′54″ S, 115°47′44″ E). Four B. napus cultivars were used, including two early-flower cultivars (AV-Ruby and ZY821), and two late-flower cultivars (Alku and YM11) (Chen et al., 2023). To synchronize anthesis across all cultivars, the early-flowering cultivars were sown 4 weeks later than the late-flowering ones.

Plants were grown in the screenhouse until the second flower on the main stem was open. At 08:00 on that day, the 2nd to 5th flowers on the main stem were pollinated with pollen collected from freshly opened flowers of the same cultivar. Then the plants were covered with selfing bags and moved to two separate heat treatment chambers for 7 days, with a control treatment with transient daily maximum/minimum temperatures of 25 °C/15 °C, and a heat stress treatment with transient maximum/minimum temperatures of 32 °C/22 °C. The plant management and chamber settings followed the protocols described in Hu et al. (2025).

After temperature treatment, plants were returned to the screenhouse and grown to maturity. At maturity, pods formed from the 2nd to 5th flowers on the main stem were collected and dried at 30 °C for 14 days. Yield-related traits, including pod length, seed number per pod, seed yield per pod and 1,000-seed weight were measured in these pods.

For transcriptome analysis, the 2nd to 5th flowers and immature pods on the main stem were pooled into a single sample and collected at 14:00 on day zero after beginning of heat treatment (DAT0), DAT1, DAT3, and DAT6, following the protocol described in Hu et al. (2025). Samples were immediately frozen in liquid nitrogen and stored at −80 °C. This experiment included three biological replicates of four cultivars sampled at four time points under both heat and control conditions, resulting in a total of 96 samples.

Total RNA extraction was followed the same methods described in Hu et al. (2025) from each of the 96 samples. A total of 10 μg of RNA from each sample was aliquoted and sent to Biomarker Technologies (BMKGENE, Hongkong, China) for mRNA purification, cDNA library construction and sequencing by using Illumina technology.

RNA sequencing was performed using the Illumina Novaseq X platform (Illumina, San Diego, CA, United States), generating high-quality sequencing reads. Raw sequencing data in FASTQ format were initially processed using in-house Perl scripts to filter out adapter sequences, ploy-N-containing reads, and other low-quality reads. After this quality control step, clean reads were obtained, and their quality metrics, Q20, Q30, GC-content and sequence duplication levels were calculated.

The clean reads were then mapped to the B. napus reference genome (v4.1) (https://www.genoscope.cns.fr/brassicanapus/data/Brassica_napus_v4.1.chromosomes.fa.gz). Only reads with no more than one single mismatch were retained for further analysis and gene annotation.

The expression levels of genes were quantified using fragments per kilobase of transcript per million mapped reads (FPKM). Gene expression levels were compared between heat and control treatments. Differential expression analysis was conducted by DESeq2 (Love et al., 2014). The p-values were adjusted using the Benjamini and Hochberg’s approach to control the false discovery rate (FDR) (Benjamini and Hochberg, 1995). DEGs were considered significant if they met the threshold of fold change ≥1.5 and FDR <0.01. Venn diagrams were generated by using the Molbiotools (https://molbiotools.com).

Gene function annotation was conducted with multiple databases, including COG (Clusters of Orthologous Genes), GO (Gene Ontology), KEGG (Kyoto Encyclopedia of Genes and Genomes), KOG (Eukaryotic Orthologous Groups of Proteins), Pfam (Protein family), SWISS-PROT (a curated protein sequence database), eggNOG database (evolutionary genealogy of genes: Non-supervised Orthologous Groups), and NCBI-nr (NCBI non-redundant protein sequences).

Functional enrichment analysis of GO terms and KEGG pathways were performed at the Brassica napus multi-omics information resource (BnIR) (https://yanglab.hzau.edu.cn/BnIR) (Yang et al., 2023). Short time-series expression miner (STEM) analysis was performed using the OmicShare tools (www.omicshare.com/tools) (Mu et al., 2024). The STEM clustering method was applied to analyse the expression patterns of AV-Ruby specific DEGs, based on the FPKM at four time points under heat treatment. The number of model profiles was set to 10, and the threshold of significant profile was set as p < 0.05.

Quantitative real-time PCR (qRT-PCR) was performed to validate the RNA-Seq findings. One ug RNA of each sample was reverse-transcribed using iScript™ gDNA Clear cDNA Synthesis Kit (Bio-Rad, Hercules, CA, United States). Thirteen DEGs were randomly selected across four time points for the qRT-PCR experiment: BnaA10g20610D, which showed the highest average expression among the three consistently upregulated DEGs, BnaA01g23000D, BnaA02g28330D, BnaA02g36510D, BnaA05g09400D, BnaA09g52180D, BnaC04g13110D, BnaC04g22890D, BnaC07g39860D, BnaC09g16520D, BnaC09g45670D, BnaCnng12730D, and BnaCnng52120D. Primers were designed for each DEG using the NCBI online website Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primerblast) (Supplementary Table S1). BnActin7 served as an internal control for normalizing gene expression levels (Chen et al., 2010). qRT-PCR was performed with iTaq™ Universal SYBR^® Green Supermix (BioRad, Hercules, CA, United States). All experiments were carried out at least three times on three biological replicates. The relative expression levels were calculated by the method of fold change of 2^−ΔΔCT value (Schmittgen and Livak, 2008).

Analysis of variance (ANOVA) was performed using R software (version 4.3.3) with the ‘tidyverse’ R package (version 2.0.0) (Wickham et al., 2019; R Core Team, 2024). Two-way ANOVA was employed to assess interactions, and pairwise comparisons of means were conducted using the least significant difference (LSD) test. A difference between two means greater than the LSD value was considered significantly different at p < 0.05 (Butler, 2022).

All hand-pollinated flowers formed fertile pods at maturity, and each pod contained at least one seed. Highly significant main effects of cultivar (C) were observed for all four yield traits (p < 0.001), while significant main effects of heat treatment (T) were detected for pod length (p < 0.001), number of seeds per pod (p < 0.001) and seed yield per pod (p < 0.001), but not for 1,000-seed weight. Additionally, significant main effects of C × T interaction (p < 0.05) were observed in pod length and seed yield per pod (Supplementary Table S2).

Among the four cultivars, AV-Ruby, showed no significant reduction in any yield-related traits under heat compared to control treatment in flowers 2–5 (Figure 1; Supplementary Table S3). In contrast, YM11 and ZY821 showed a significant reduction in pod length after heat treatment, at 81.1% and 83.1% of the control, respectively. Moreover, Alku, YM11 and ZY821 experienced a notable decline in the number of seeds per pod and seed yield per pod after heat stress. ZY821 was the most heat-sensitive cultivar with seed yield per pod in the heat treatment was only 62.8% of control. However, 1,000-seed weight was not significantly affected by heat treatment in all four cultivars (Figure 1; Supplementary Table S3).

After filtering and trimming, an average of 46.7 million high-quality clean reads were retained for further analysis. The average Q20 and Q30 values were 97.7% and 93.6%, respectively (Supplementary Table S4). The clean reads were mapped to the B. napus reference genome with an average mapping rate of 87.78%.

Functional annotation was performed using multiple databases, with the following number of genes assigned: 26,566 (22.7%) to COG, 74,968 (63.9%) to GO, 63,879 (54.5%) to KEGG, 50,575 (43.1%) to KOG, 73,628 (62.8%) to Pfam, 69,008 (58.9%) to SWISS-PROT, 77,641 (66.2%) to eggNOG, and 109,161 (93.1%) to NCBI-nr. In total, 109,317 (93.2%) genes were successfully annotated across these databases (Supplementary Table S5).

Across the 4 sampling times and 4 cultivars, a total of 36,933 DEGs were identified in heat vs. control treatments (Supplementary Table S6). The total number of DEGs followed a regular temporal pattern across the four cultivars, with an increase from DAT0 to DAT3, then a decline from DAT3 to DAT6. This pattern occurred in both upregulated and downregulated DEGs which reached peak levels at DAT3 in all four cultivars. At DAT0, AV-Ruby exhibited the highest number of total DEGs, including both up- and downregulated DEGs, among the four cultivars. At DAT0, a total of 3,436 DEGs were identified in AV-Ruby, with 1,270 upregulated and 2,166 downregulated. On the contrary, at DAT6, AV-Ruby displayed the lowest DEG count (2,713 in total) among the four cultivars (Figure 2A).

Three DEGs - BnaC05g17990D, BnaC04g12620D, and BnaA10g20610D - encoded HSP and HSP-like proteins and were consistently upregulated across all four time points and in all four cultivars (Supplementary Table S6).

BnaC05g17990D encoded a HSP20-like chaperones superfamily protein. In AV-Ruby, BnaC05g17990D peaked earlier than in the other three cultivars at DAT1 (log₂[fold change] = 3.06), while its peak expression occurred at DAT3 in Alku (log₂[fold change] = 4.18), YM11 (log₂[fold change] = 3.55), and ZY821 (log₂[fold change] = 3.84) (Figure 2B; Supplementary Table S6).

BnaC04g12620D encoded a 70-kDa HSPA8, and was consistently upregulated in YM11 across all time points (average log₂[fold change] = 1.30) (Figure 2B; Supplementary Table S6).

BnaA10g20610D encoded a 17.6-kDa class II HSP, and exhibited the highest average expression levels among these three genes, but its expression varied considerably across all cultivars and time points. YM11 showed an increasing trend over time, while the highest expression was detected at DAT0 in Alku (log₂[fold change] = 3.74), at DAT1 in AV-Ruby (log₂[fold change] = 4.21), and at DAT3 in ZY821 (log₂[fold change] = 5.67) (Figure 2B; Supplementary Table S6).

At DAT0, DEGs were significantly enriched in several GO terms related to cell growth and development and to pollen tube growth. At DAT1, DEGs were notably enriched in GO terms related to unfolded protein binding and hormonal signalling. At DAT3, DEGs were enriched in biological process. Notably, response to high light intensity DEGs were consistently enriched across DAT0, DAT1, and DAT3 (Table 1; Figure 3A).

Among the top 20 GO terms at DAT6, 11 were related to secondary metabolites, including flavonoids, phenylpropanoids, and glucosinolates. Additionally, four GO terms at DAT6 were associated with DNA replication (Table 1; Figure 3A).

KEGG pathway enrichment analysis revealed several pathways consistently enriched across all four time points, including “glyoxylate and dicarboxylate metabolism,” “carbon fixation in photosynthetic organisms,” and “alpha-linolenic acid metabolism.” In the meantime, some pathways were uniquely enriched at specific time points. At DAT0, the pathway “lectins” was exclusively enriched, while DAT1 showed specific enrichment in “carotenoid biosynthesis.” At DAT3, the “autophagy–other” pathway was uniquely enriched, and DAT6 exhibited significant enrichment in “DNA replication” (Figure 3B).

Photosynthesis-related pathways, including “photosynthesis proteins,” “photosynthesis - antenna proteins,” and “photosynthesis,” were notably enriched at both DAT0 and DAT3. Additionally, several antioxidant-related pathways were highlighted, such as “glutathione metabolism” at DAT1 and “flavonoid biosynthesis” at both DAT0 and DAT6 (Figure 3B).

In this experiment, there were 4,741 DEGs specifically identified in AV-Ruby (Figure 4A) which were further classified by STEM analysis into ten profiles based on their expression patterns under heat stress. The most significantly enriched profiles (p < 0.05), accounting for over 10% of the grouped DEGs, were profile 5 (1,107 genes), profile 0 (888 genes), profile 9 (745 genes), and profile 6 (507 genes) (Figure 4B).

DEGs in profile 5 showed a late-response pattern, with no significant change at DAT0 but gradually increased from DAT1, peaked at DAT3, and then declined. In contrast, DEGs in profile 0 exhibited high expression level at DAT0, followed by a steady decline from DAT0 to DAT6, and DEGs in profile 6 exhibited an early-response pattern, with an initial spike in expression at DAT1, followed by a sharp downregulation. DEGs in profile 9 exhibited a delayed or gradual response, starting with low expression levels that increased over time.

To investigate the functional significance of gene expression changes, GO and KEGG enrichment analyses were performed on DEGs in profiles 5, 0, 9 and 6 (the most populated expression profiles).

In profile 5, 16 among the top 20 GO terms were related to RNA processing and DNA repair pathways, including “RNA methyltransferase activity” and “double-strand break repair via homologous recombination.” DEGs in profile 0, which were steadily downregulated over time, were predominantly enriched in cell growth and calcium-dependent kinase pathways, such as “regulation of pollen tube growth,” “regulation of cell development,” “calcium-dependent protein serine/threonine kinase activity” and “calcium-dependent protein kinase activity.” DEGs in profile 9, which demonstrated a gradual upregulation pattern, significantly enriched in seven GO terms, with the highest enrichment in “pyridoxal phosphate binding.” Meanwhile, DEGs in profile 6, which showed an early-response pattern with a sharp initial increase followed by a decline, were associated with vesicle trafficking, oxidative stress response, membrane fusion, and programmed cell death, suggesting their involvement in heat stress adaptation (Figure 5A).

KEGG pathway enrichment analysis further highlighted the functional differences among these profiles. DEGs in profile 5 were specifically enriched in pathways related to ribosome biogenesis and DNA repair, including “ribosome biogenesis in eukaryotes,” “mismatch repair,” “base excision repair,” “nucleotide excision repair,” and “RNA degradation.” DEGs in profile 0 were strongly associated with metabolic and antioxidant pathways, including “citrate cycle (TCA cycle),” “peroxisome,” and “ascorbate and aldarate metabolism“ Profile 9 was linked to various metabolic processes, such as “sulphur metabolism,” “glyoxylate and dicarboxylate metabolism,” and “amino sugar and nucleotide sugar metabolism.” DEGs in profile 6 were uniquely enriched in “SNARE (soluble N-ethyl maleimide sensitive factor attachment protein receptor) interactions in vesicular transport” (Figure 5B).

To further characterize the transcriptional responses underlying heat resilience, we examined the top 15 up- and downregulated AV-Ruby–specific DEGs across four time points (Figure 6; Supplementary Table S8). At DAT0, BnaA03g03760D (unknown protein) and BnaA05g09400D (light-harvesting chlorophyll-protein complex II subunit B1) (Supplementary Table S6) were markedly upregulated in AV-Ruby, with log₂(fold change) values of 2.20 and 1.94, respectively, whereas both DEGs were downregulated in heat-sensitive cultivars Alku, YM11, and ZY821 under heat stress relative to the control treatment.

Several of the top 15 AV-Ruby–specific DEGs displayed significant expression changes across multiple time points. For example, BnaA04g05090D, a cytochrome P450-related gene (Supplementary Table S6), was upregulated at both DAT0 (log₂[fold change] = 2.20) and DAT1 (log₂[fold change] = 2.69). Similarly, BnaC04g22890D, associated with GENE SILENCING 3, showed strong heat-induced expression at DAT3 (log₂[fold change] = 6.52) and DAT6 (log₂[fold change] = 4.16). In contrast, BnaC03g60950D, related to PsaB, was significantly upregulated at DAT0 (log₂[fold change] = 2.33) but notably downregulated at DAT6 (log₂[fold change] = 5.06). Likewise, BnaUnng00820D (unknown protein, Supplementary Table S6) showed increased expression at DAT3 (log₂[fold change] = −3.35) but downregulated at DAT6 (log₂[fold change] = −4.23).

The expression patterns of thirteen DEGs randomly selected across four time points to validate the RNA-seq expression data were generally consistent with the RNA-seq results (Figure 7). As expected, the DEG that was common to all cultivars (BnaA10g20610D) showed consistent upregulation across all samples, confirming the RNA-seq findings (Figure 7A). Among the 12 randomly selected DEGs, seven (BnaA01g23000D, BnaA02g28330D, BnaA02g36510D, BnaA09g52180D, BnaC04g22890D, BnaC09g45670D and BnaCnng12730D) exhibited similar expression trends between qRT-PCR and RNA-seq. However, five genes (BnaC09g16520D, BnaA05g09400D, BnaC07g39860D, BnaCnng52120D and BnaC04g13110D) displayed discrepancies in expression direction. For example, BnaC09g16520D was downregulated under heat stress in the RNA-seq data in YM11 and ZY821 but upregulated in the qRT-PCR analysis (Figure 7B).