Seeds of spring oilseed rape (Brassica napus L., cv. Westar) were surface-sterilised and sown into seed trays filled with a seed potting mix (Clover Peat, Dungannon, Northern Ireland). Trays were then placed in polyethylene tunnels for 4 weeks at the Crops and Environment Laboratory, University of Reading, to germinate. Twenty-eight days after sowing (DAS), the plants were transplanted into 3-L pots containing a peat based potting mix (Clover Peat).

At green bud stage (38 DAS), the plants were moved to two controlled environment chambers [Fitotron, Weiss Technik (UK) Ltd], for a period of 1 week to adjust to the new environment before the start of the heat experiment. Each cabinet contained 16 plants and was maintained under control conditions of 20°C/14°C day/night and a photoperiod of 16/8-h light/dark.

Since the flowering stage of winter-sown B. napus occurs during May in the UK, an analysis of daily weather temperatures during this developmentally crucial period was carried out using the UK Meteorological (Met) Office weather data obtained from the Met office (n.d.) and from WeatherOnline (n.d.) platforms. The data showed an increase in the frequency of high temperature fluctuations during May, ranging between 26°C and 28°C over several days ( Supplementary Figure S1 ). Given the continuous rise in global temperature, it is expected that future heatwaves would increase in severity and duration (Seneviratne et al., 2012). Based on these data, the current study temperature simulates a heat stress event through a temperature increase to 30°C/24°C day/night for a period of 6 days and 20°C/14°C day/night for control. Forty-five DAS, the temperature in one cabinet was kept under control conditions (20°C/14°C) and denoted as control cabinet, while the second cabinet was set to a gradual increase in temperature and denoted as heat treatment (HT) cabinet. To mimic a heatwave, the temperature in the HT cabinet was increased gradually from 20°C to 30°C between 9:00 and 12:00 in three steps: 20°C at 9:00 h, 24°C at 10:00 h, 28°C at 11:00 h, and 30°C at 12:00 h. The temperature was held at 30°C until 22:00 h ( Supplementary Figure S2B ), before it was dropped gradually and maintained at 24°C until 8:00 of the next day ( Supplementary Figure S2B ). This cycle of gradual increase (day) and decrease (night) of temperature was held for 5 days ( Supplementary Figure S2C ). On the sixth day of treatment, the temperature was gradually decreased to 20°C/14°C day/night cycle ( Supplementary Figure S2D ) and held for a recovery period of 7 days ( Supplementary Figure S2E ). During heat treatment, the plants were frequently irrigated to avoid drought.

Five biological replicates from each treatment condition were collected at days 0, 1, 2, 6, and 12 of the experiment. Each replicate was made up of three individual plants. Three flower buds were collected per plant (nine buds per replicate). All samples were immediately snap frozen in liquid nitrogen before they were stored at −80°C. Bud samples were ground into fine powder using pestle and mortar, with the frequent addition of liquid nitrogen to stop enzymatic reactions. The ground material was then stored at −80°C until RNA extraction.

Total RNA from bud samples was extracted using the Spectrum™ Plant Total RNA Kit (Sigma-Aldrich Dorest, UK) in accordance with the manufacturer’s protocol and treated with genomic DNAse using DNASE 70-On Column DNase I Digestion set (Sigma-Aldrich) to eliminate DNA contamination. RNA quantification was estimated on both NanoDrop 2000 (Thermo Scientific) and Qubit 2.0 fluorometer (Invitrogen, USA), and its quality was evaluated on 1% (w/v) denaturing formaldehyde agarose gel (MOPS). Samples were shipped over dry ice to Novogene Europe, where sequencing was performed on the Illumina Novaseq™ 6,000 (PE150) platform.

Raw sequence reads were assessed for quality using FastQC tools (v0.11.5) (Andrews, 2012). The mean sequence lengths were 150 bp, and the mean sequence GC content was 43%. The mean quality scores in each base position were higher than 36, and the mean quality scores per sequence were 36. As a result, sequence trimming was not necessary.

STAR software (v2.7.10a) was used to map the clean reads to the B. napus cv. Westar v0 reference genome (Song et al., 2020), and the percentage of aligned reads was calculated by using the flagstat command from samtools on each alignment file generated by STAR aligner ( Supplementary Figure S3 ). To estimate transcript abundance, RNA-Seq by Expectation-Maximisation (RSEM) software (v1.2.18) was used in three steps to prepare the reference, calculate expression, and generate the count matrices. First, rsem-prepare-reference script was used with “–gff3” option to extract reference sequences from the genome and “—star” option to build STAR indices. Next, rsem-calculate-expression script was used with “–paired-end” option to align input reads against the reference transcriptome with STAR and calculated expression values using the alignments. Finally, rsem-generate-data-matrix script was used to generate the count matrices from expression results files.

For differential expression analysis, DESeq2 (v1.41.10) was used in R environment (Love et al., 2014) with default parameters. A pre-filtering step was performed to keep transcripts that have a minimum count of 10 reads in a minimum of four samples. To identify transcripts that were significantly differentially expressed, the two conditions (heat stress versus control) were contrasted in pairwise comparisons. Thresholds of |log2(foldchange)| ≥ 0.5 and adjusted p-value < 0.05 (using the Benjamini and Hochberg method) were considered to identify significantly differentially expressed transcripts between the two treatment conditions.

Significant DEGs in each contrast were further analysed using K-mean clustering based on the kmeans() function in R (distance: Euclidean), with log2FC as the input.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted on the DEGs using OmicsBox (BioBam Bioinformatics, 2019). The significantly enriched GO terms and biological pathways in HT samples as compared to control at different timepoints were identified to determine heat-stress-related functions and pathways. The analysis was carried out using Fisher exact statistical test with FDR adjustment cutoff <0.05. The background dataset consisted of all B. napus cv. Westar identifiers present in the assembly’s annotation file (Song et al., 2020).

Sucrose Synthase 5 (SS5) and Heat Shock Protein 20 (HSP20) were chosen from the pathways that were involved in the heat stress response to further examine their expression in real time using real-time quantitative PCR (RT-qPCR). Five timepoints were chosen for each gene to represent the time course of the experiment. SS5 and HSP20 were not analysed at 1 day after treatment (DAT), as they exhibited similar expression profile at 0 and 1 DAT; therefore, one timepoint was selected ( Supplementary Table S2 ).

Total RNA used for transcriptomic sequencing was used for cDNA synthesis. First, cDNA was synthesised from 1 µg total RNA by Invi trogen SuperScript IV VILO Master Mix kit according to the manufacturer’s instructions (Thermo Fisher Scientific). The RT-qPCR reaction volume was 20 µL consisting of 3 µL cDNA, 10 µL SYBR Premix Ex Taq II (Takara, China), 1 µL (200 nM final concentration) of 4 µM of each primer (forward and reverse, Supplementary Table S3 ), and 5 µL RNAse-free water. Three technical and three biological replicates for each sample were measured with the following protocol: 50°C for 2 min, 95°C for 2 min, followed by 40 cycles of 95°C for 3 s, 60°C for 30 s, followed by dissociation curve analysis ramping from 60°C to 95°C with a ramp rate of 0.3°C s⁻¹ on an ABI StepOne Plus RT-qPCR platform (Applied Biosystems, USA). BnaActin was used as housekeeping gene to normalise the data. The relative expression level of all selected genes at each timepoint was calculated using the 2^−ΔΔCT method of the StepOne software (v2.2.2).

Five biological replicates from each treatment condition were collected at days 0, 1, 2, 5, 6, 8, and 12 of the heat treatment experiment. Each replicate was made up of three youngest fully expanded leaves collected from three plants. All leaf samples were immediately snap frozen in liquid nitrogen before being stored at −80°C until further analysis. Frozen leaf samples were freeze-dried for 5 days, then ground into a fine powder using a Precellys 24 lysis and homogenisation (Stretton Scientific, Alfreton, UK) and stored at −40°C in sealed sample bags until extraction and analysis.

GLS extraction was carried out as per the protocol of Bell et al. (2015) as follows: 40 mg of ground leaf powder was heated in a dry block at 75°C for 2 min. This step was done as a precautionary measure to inactivate as much myrosinase enzyme as possible before extraction (Pasini et al., 2012). Afterwards, 1 mL of preheated (70% v/v) methanol (70°C) was added to each sample and placed in a water bath for 20 min at 70°C. Samples were then centrifuged for 5 min (6,000 rpm at 18°C) to collect the loose material into a pellet. The supernatant was then transferred into fresh labelled Eppendorf tube and stored at −80°C until further analysis. Before LC-MS analysis, the samples were filtered using 0.25-µm filter discs, diluted in 50% methanol at a 1:4 ratio (dilution factor =5) and spiked with 100 µL of the internal standard Sinigrin (100 ng mL⁻¹). The GLS content was analysed by SCIEX QTRAP 6500+ LC-MS.

For LC separation, a Waters Acquity BEH C18 column (particle size, 1.7 µ, 2.1 × 50 mm) with a Security VanGuard system from Waters (UK) was used. The mobile phase consisted of water with 0.1% formic acid (A) and methanol with 0.1% formic acid (B). The gradient started at 5% B and was raised to 90% B in 3.5 min, held at 90% B for 0.5 min and re-equilibrated at 5% B for 1 min. The total time of the gradient program was 5 min. The flowrate was 0.4 mL min⁻¹ with a column temperature of 60°C. A 1-µL aliquot of sample was injected for analysis.

The LC system was interfaced with a SCIEX QTRAP 6500+ mass spectrometer equipped with an IonDrive Turbo V ion source. Multiquant software was used to control sample acquisition and data analysis. The QTRAP 6500+ mass spectrometer was tuned and calibrated according to the manufacturer’s recommendations. For quantification, an external standard calibration curve was prepared using a series of standard samples containing the following GLS: Gluconapin (GNA), Progoitrin (PRO) and Glucobrassicin (GBS), with concentrations of 1,000, 500, 250, 100, 50, 10, and 1 ng mL⁻¹. Two-way mixed ANOVA in R was performed at 0.05 significance level to calculate the significance variation in concentration under heat treatment as compared to control.

The extraction of soluble sugars (glucose, fructose, and sucrose) was performed as follows: 50 mg of ground leaf powder was extracted with 1 mL of 62.5% methanol/37.5% HPLC-grade water (v/v) at 55°C over a period of 15 min in a shaking water bath and vortexed every 5 min. Afterwards, the samples were allowed to cool for 2 min and then centrifuged at 13,000 rpm for 10 min and filtered using 0.25-µm filter discs. The supernatant was then transferred to a clean labelled Eppendorf tube and stored at −80°C until further analysis. Before HPLC analysis, the samples were diluted in HPLC-grade water in a 1:2 ratio (dilution factor = 3). Chromatographic separation of sugars content was performed using HPLC (Agilent 1260, Infinity Series) equipped with an Evaporative Light Scattering Detector (ELSD) system and an Asahipak NH2P-50 4E column (250 × 4.6 mm; Shodex, Tokyo, Japan) in an isocratic elution mode. The mobile phase consisted of 75% acetonitrile/25% water at a flowrate of 1 mL min⁻¹, a column temperature of 40°C, and an injection volume of 20 µL. The analysis had a total run time of 25 min. For quantification, an external standard calibration curve was prepared using a series of standard samples containing glucose, fructose, and sucrose with the concentrations: 0.025, 0.05, 0.1, 0.5, and 1 mg mL⁻¹.

Following the treatments, the plants were allowed to complete their lifecycle, and total aboveground biomass and seed yields were collected from individual plants. Seed samples were scanned with a DA 7250 NIR analyser (PerkinElmer, Beaconsfield, UK). The seed samples were scanned for moisture, oil, and protein content and for fatty acid composition of the oil.

Two-way mixed ANOVA was performed at 0.05 significance level to calculate the statistically significant differences between the means of sugars and glucosinolate concentration under heat treatment as compared to control. For a significant two-way interaction between treatment and timepoint, the following analysis were conducted: the simple main effects of treatment on sugars concentration were examined at each timepoint, using Bonferroni-adjusted p-values for multiple comparisons. This was followed by pairwise comparisons between treatment levels (heat stress vs. control) at each timepoint to determine whether the mean sugars concentrations differed significantly between heat stress and control. In addition, near-infrared spectroscopy (NIRS) data were analysed for the investigated components using a t-test. All statistical analysis was performed using R.

Heat treatment significantly altered plant growth and appearance, with visible physiological changes on the flowers and buds ( Figure 1 ). At 5 DAT, flowers displayed smaller and paler petals compared to plants growing under control conditions. At 3 days of recovery (DOR), many buds suffered from abscission, appeared brown, and shrunken in size on plants subjected to the heat treatment. At maturity, heat-treated plants exhibited less total biomass and total seed weight relative to control plants ( Supplementary Figure S4 ). In terms of oil quality, heat treatment significantly decreased linoleic acid, linolenic acid, and palmitic acid concentration but increased oleic acid and stearic acid concentration in seeds compared to control plants ( Supplementary Figure S5 ).

Libraries were constructed and sequenced on the Illumina NovaSeq platform. An average of 47 million reads were generated per sample after low-quality reads were filtered out ( Supplementary Table S1 ). Following this cleaning step, approximately 90%–95.49% total reads per sample were mapped to the B. napus cv. Westar v0 reference genome ( Supplementary Figure S3 ) (Song et al., 2020). Principle component analysis (PCA) was used to analyse the samples that drove the group separation ( Supplementary Figure S7 ). PC1 and PC2 captured 48.9% and 8.2% of the total variance of the samples, respectively. The analysis showed that both control and treatment samples clustered near each other at 7 DOR, which demonstrated that heat-treated plants might have started to return to normal after 7 days of recovery.

To understand the effect of heat stress on gene expression in B. napus, differentially expressed genes (DEGs) were identified in a pairwise comparison between heat treatment and control conditions at five timepoints during and after heat treatment. A total of 1,892, 3,253, 4,553, 4,165, and 1,713 genes were differentially expressed at days 0, 1, and 2 DAT and 1 and 7 DOR, respectively ( Supplementary Figure S8 ). Analysis of the top differentially expressed genes after 24 h of heat treatment revealed that eight transcripts with logFC ranging between 4.4 and 10.8 mapped to small heat shock proteins (sHSP20), in addition to transcripts mapped to other heat shock factors such as Elongation Factor 1-beta 1-like and Chaperone Protein ClpB1 ( Supplementary Table S7 ). On the other hand, top downregulated genes featured transcripts mapped to sugars transporters, nitrogen transport and storage, cell wall modification, and methylation ( Supplementary Table S8 ). Similarly, after 48 h of heat treatment, many of the identified transcripts remained among the top DE, both up- and downregulated, highlighting the significant effect of heat stress on the functions performed by these enzymes and proteins.

To identify trends in the expression of genes across the time course of the experiment, K-means clustering was applied to 9,933 DE genes with logFC ≥ 0.5. The algorithm randomly assigned each gene into one of the clusters based on the Euclidean distance between the gene and the cluster mean. Six distinctive clusters based on expression changes across the five timepoints were identified ( Figure 2 ). Interestingly, clusters 1 and 3 (475 and 607 genes), showed opposite expression patterns towards the end of the heat treatment and the start of the recovery phase. In cluster 1, the expression profile of genes decreased at 1 DOR (HS6vsC6), before returning to normal levels at 7 DOR (HS12vsC12), while in cluster 3, a sharp increase in expression was observed at 1 DOR (HS6vsC6), which returned to normal levels at 7 DOR (HS12vsC12).

In clusters 2 and 5 (3,470 and 1,333 genes), genes were downregulated along the course of the experiment but more pronounced during the heat treatment phase especially in cluster 5, which exhibited a greater decrease in expression. Both clusters 4 and 6 (398 and 3,518 genes) showed increased expression during the treatment before it started to decline towards the end, reaching normal levels at 7 DOR (HS12vsC12) ( Figure 2 ).

To verify the reliability of RNA sequencing results, two genes with diverse expression profiles, including upregulated or downregulated at different timepoints of the experiment, were selected for real-time qPCR to measure expression levels. As a result, Heat Shock Protein 20 (HSP20) and Sucrose Synthase 5 (SS5), exhibited similar expression profiles between RT-qPCR and RNA-seq data ( Supplementary Figure S9 , Supplementary Table S4 ), with correlation coefficients of r=0.86 and r=0.68, respectively.

Gene Ontology (GO) enrichment analysis of the DEGs was performed to identify enriched GO terms. Terms such as “binding,” “catalytic activity,” “metabolic process,” “hydrolase activity,” and “ion binding” were among the predominant enriched terms at all timepoints with more transcripts mapped to these terms during heat treatment than during recovery ( Supplementary Figures S10 - S14 ). Oxidoreductase, hydrolase, and transferase were among the top identified enzymes ( Supplementary Figure S15 ). These results suggest that high temperature activates several metabolic processes through the expression of numerous enzymes involved in alleviating the impact of heat stress not only during the stress period but also during recovery.

Based on Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis (Kanehisa and Goto, 2000), pathways in the respiratory metabolism, namely, glycolysis and pentose phosphate pathway were enriched in all heat treatment comparisons at all timepoints, while citrate cycle and oxidative phosphorylation were enriched at 2 DAT and 1 DOR, respectively ( Supplementary Figures S16 - S19 ). Other related pathways were also enriched along the course of the treatment and during recovery period, and these include one carbon pool by folate (enriched at 1 and 2 DAT and 1 DOR), carbon fixation in photosynthetic organisms, cysteine and methionine metabolism (enriched at 1 and 2 DAT), glycerolipid and tryptophan metabolism (enriched at 2 DAT), fatty acid biosynthesis and pyruvate metabolism (enriched at 1 DOR), and fructose and mannose metabolism (enriched at 7 DOR). In the present study, most of the genes involved in these pathways were responsive to high temperature stress, with their expression declining with the onset of heat treatment, especially after 48 h ( Supplementary Figures S20 - S27 ).

In this study, comparative transcriptome analysis between heat treatment and control identified 9,933 significantly differentially expressed transcripts (|log2(foldchange)|≥0.5, FDR ≤ 0.05) ( Supplementary Figure S8A ) with a balanced proportion of up- and downregulated genes in contrast HS0vsC0 (0 DAT), and more downregulated than upregulated genes in the rest of the contrasts were identified. This suggests that the impact of heat treatment on the plant’s transcriptome after 24 h of heat stress, resulted in downregulation of many genes involved in multiple functional and metabolic processes such as glycolysis, PPP, citrate cycle, and oxidative phosphorylation. In contrast, HS12vsC12 (7 DOR), although the proportion of downregulated genes was higher than the upregulated ones, the total number of DEG decreased dramatically as the plant started to restore its metabolic functions after 7 days of recovery.

Differential expression analysis showed that some of the genes encoding sugar transporters, nitrogen transport and storage, cell wall modification, and methylation were downregulated ( Supplementary Tables S8 , S10 ). In contrast, several genes mapping to small heat shock proteins (sHSP), other heat shock transcription factors, methylation, and redox balance-related genes were upregulated ( Supplementary Tables S7 , S9 ). Among the top DE transcripts were several transcripts that mapped to HSPs and other HSFs. Notably, members of the sHSP family, also known as the HSP20 family proteins, such as HSP17.6, HSP22, and HSP23.6, showed high expression levels throughout the heat treatment, with expression reaching a log2FC of 11 at 2 DAT. This protein family is the most abundantly produced in plants in response to heat stress, and their role in plant responses to various abiotic stresses has been established (Sun et al., 2020). Other studies also reported the involvement of HSFs in different biotic and abiotic stress (Mishra et al., 2002, 2019; Chowdhary et al., 2023). In A. thaliana, thermoprotection tests demonstrated the role of HSFA6b in thermotolerance acquisition (Huang et al., 2016). Moreover, Giorno et al. (2009) showed that HsfA2 and Hsp17-CII are induced under high temperature in tomato, which are responsible for activating protection mechanisms during heat stress.

The analysis of KEGG-enriched pathways revealed the main underlying biochemical pathways altered in B. napus in response to heat stress ( Supplementary Figures S16 - S19 ). In accordance with other reports (Estravis-Barcala et al., 2021; Li et al., 2022), DEGs showed an enrichment of pathways related to respiratory metabolism, namely, glycolysis, citrate cycle, oxidative phosphorylation, and pentose phosphate pathway (PPP). At 2 DAT, the impact of heat stress on glycolysis was manifested by the downregulation of the key enzymes phosphoglucomutase, aldose 1-epimerase, hexokinase, glucose-6-phosphate isomerase, and 6-phosphofructokinase ( Supplementary Figure S21 ). Moreover, most of the citrate cycle genes were also downregulated at 2 DAT ( Supplementary Figure S22 ). In the electron transport chain (ETC), gene expression analysis showed that most of the involved genes were downregulated under heat treatment, except for succinate dehydrogenase (SDH), which was upregulated ( Supplementary Figure S23 ). In mitochondrial metabolism, SDH also known as complex II, plays a central role in the citrate cycle and ETC (Huang and Millar, 2013). Additionally, SDH is thought to impact different components of the plant stress response, including stomatal conductance and ROS scavenging (Huang et al., 2019).

Although respiratory pathways have been found to have potential roles in adaptive response to abiotic stress, mainly through reducing the production of ROS (Van Dongen et al., 2011), under extreme conditions such as severe or prolonged heat stress, respiratory enzymes are deactivated, and proteins are denatured leading to culminating of ROS and a total breakdown of mitochondrial respiration (Scafaro et al., 2021). This negatively impacts oxygen and carbon fluxes and eventually leads to more severe yield penalty (Rasmusson et al., 2020). In line with this, results from heat stress experiment on seagrass showed that high midday temperature stress of 40°C was associated with significant decrease in biomass (George et al., 2018). Likewise, Impa et al. (2019) found that under high night temperature, wheat tolerant cultivar exhibited smaller reductions in biomass and lower rates of both net photosynthesis and respiration compared to other wheat cultivars (Impa et al., 2019).

In the PPP, the oxidative phase (OPPP) was mainly affected at 1 DAT, and this was indicated by the differential expression of glucose-6-phosphate dehydrogenase (G6PDH) and 6-phosphogluconate dehydrogenase (6PGDH) enzymes ( Supplementary Figure S24 ). G6PDH and 6PGDH are the key enzymes in the OPPP. They catalyse the first and third steps of the pathway, respectively (Long et al., 2016), and have been reported to be associated with the response to various abiotic stresses in plants (Gong et al., 2012). Esposito (2016) reported the involvement of the OPPP in the early response to abiotic stress factors, forming a true metabolic sensor to oxidative stress. At 2 DAT, most of the oxidative and non-oxidative phase enzymes were downregulated, including transketolase (TK), transaldolase, and phosphofructokinase enzymes ( Supplementary Figure S25 ). TK has a central role in primary metabolism, where the products of the involved reactions produce the precursors for nucleic acids biosynthesis, aromatic amino acids, and vitamins (Bi et al., 2013). Bi et al. (2013) showed that the cucumber TK gene (CsTK) was sensitive to temperature and light (Bi et al., 2013, 2019). Likewise, the transaldolase gene was found to be involved in the regulation of expression of genes involved in the ABA signalling pathway and the enzymes responsive to ROS (Rong et al., 2021). The results of this study add further evidence to the involvement of the PPP enzymes in the heat stress response in B. napus.

At the transcriptomic level, heat treatment differentially affected the expression of GLS synthesis and regulatory genes ( Figures 3 , 4 ). Several genes encoding different UDP glycosyltransferases and glutathione S-transferases appeared downregulated, while others encoding GLS-related transcription factors were upregulated, suggesting that these genes might have participated at different stages of the stress response.

The influence of heat stress on individual GLS was further evaluated through metabolomic analysis of two aliphatic GLS [progoitrin (PRO) and gluconapin (GNA)] and one indolic GLS [glucobrassicin (GBS)]. The role of these GLS in response to different abiotic stress factors has been previously reported (Jasper et al., 2020; Ljubej et al., 2021). In the current study, the concentration of PRO and GNA significantly increased in response to heat stress before it started to decline at 5 DAT ( Figures 6 , 7 ). This decline could have resulted from membrane damage caused by prolonged exposure to heat, leading to GLS degradation. In addition, the removal of heat stress after 5 DAT could also be another factor for the decline in GLS content. In addition to its primary role in response to plant pathogen interaction (Hopkins et al., 2009; Touw et al., 2020), GLS accumulation in response to heat stress has also been previously reported (Valente Pereira et al., 2002; Martínez-Ballesta et al., 2013; Jasper et al., 2020). When plants are stressed, growth is reduced, and carbon utilisation is predominantly diverted towards the production of secondary metabolites. As part of the plant defence mechanism, GLS play important role as osmoprotective compounds (Martínez-Ballesta et al., 2013), where their concentration increases in response to stress, providing protection against oxidative damage. At 30°C/15°C (day/night) temperature regime, Brassica oleracea seedlings had significantly higher GLS concentration than plants cultivated at lower temperatures (22°C/15°C and 18°C/12°C) (Valente Pereira et al., 2002). Similarly, an Arabidopsis GLS mutant experienced a remarkable decline in growth and development and was found to be less heat tolerant than wild-type plants upon exposure to high temperature (Ludwig-Müller et al., 2000).

In contrast, indolic glucosinolate GBS showed a downward trend where concentration decreased gradually in response to heat ( Figure 8 ). Similar results have been found where GBS concentration was reduced in response to heat treatment (Bohinc and Trdan, 2012). In addition, during short-term high-temperature stress, the transcript level of the indolic glucosinolates synthetic genes was predominantly downregulated (Rao et al., 2021). Several studies reported that indole GLS are much more sensitive to heat treatment (Bones and Rossiter, 2006; Bohinc and Trdan, 2012) and can demonstrate more thermal degradation than aliphatic GLS at lower temperatures (Orlemans et al., 2006).

The results of our study indicate that the concentration of GLS in the leaves of B. napus is influenced by high-temperature stress, which was previously reported in other Brassica species (Rask et al., 2000; Bohinc and Trdan, 2012). It also shows that different groups of GLS (aliphatic versus indolic) differ in their response to heat stress. As reported by Ljubej et al. (2021), the different accumulation trends of these groups could correspond to different protective mechanisms driven by individual GLS with diverse chemical structures and bioactivities (Ljubej et al., 2021). These mechanisms may prioritise energy towards aliphatic GLS production, which are known to play crucial roles in ROS mitigation (Martínez-Ballesta et al., 2013). Additionally, the reduced availability of GLS precursors, such as nitrogen and sulphur, under stressful conditions may influence the allocation of GLS across plant tissues. In response, plants might break down their GLS and utilise the released sulphur to support primary metabolic processes, such as protein synthesis, within specific organs (Brown et al., 2003; Martínez-Ballesta et al., 2013). This plasticity highlights GLS players in balancing stress resilience and adaptation. By selectively modifying GLS content, plants can fine-tune their responses to stress, meeting immediate defence demands while ensuring long-term survival and reproductive success.

In the current study, the transcript levels of sulphur assimilation and transport genes were significantly downregulated especially at 2 DAT ( Figure 5 ). This included several adenylyl sulphate kinases and reductases, ATP sulphurylase (ATPS), and sulphate transporter genes. These genes are important for the activation, catalysis, and transport of sulphur compounds (Capaldi et al., 2015). Moreover, enrichment pathway analysis showed that cysteine and methionine metabolism pathway was among the top significantly enriched pathways at 1 and 2 DAT ( Supplementary Figures S16 , S17 ). In this pathway, 11 transcripts mapping to S-adenosylmethionine synthetase and 11 transcripts mapping to homocysteine S-methyltransferases were downregulated ( Supplementary Figure S20 ). These enzymes are among the key enzymes that control the methionine pool (Capaldi et al., 2015). Given that GLS compounds are rich in sulphur, assimilation of sulphur is closely associated with GLS biosynthesis. Cysteine, which is the terminal metabolite in sulphur assimilation, acts as a sulphur donor for methionine, a precursor for aliphatic GLS. This renders cysteine an important intersection point between sulphur assimilation and GLS biosynthesis (Rao et al., 2021). Such proximity between sulphur and these amino acids makes cysteine and methionine metabolism pathway directly related to GLS biosynthesis.

In the present study, heat treatment altered carbohydrate concentration as the amounts of glucose, fructose, and sucrose increased with heat treatment ( Figures 11 - 13 ). At the transcriptional level, the starch hydrolysing enzyme alpha-amylase was induced. Conversely, the starch-synthesizing enzymes such as phosphoglucomutase (PGM), AGPase, and hexokinase were inhibited ( Supplementary Figure S21 , S27 ). Likewise, heat induced the upregulation of sucrose-synthesizing enzymes and the downregulation of sucrose-hydrolysing enzymes such as cell wall invertase and sucrose synthase ( Supplementary Tables S5 , S6 ). These results suggest that heat-induced starch hydrolysis and promoted sucrose synthesis. Furthermore, heat stress downregulated many sugars transporters such as SWEET and ERD6-like ( Figures 9 , 10 ) and led to sugar accumulation in the leaves. Several studies have shown that various environmental stresses significantly impact sugar metabolism and transport in plants (Xalxo et al., 2020; Mathan et al., 2020). For instance, drought and salinity stress have been reported to increase sucrose content in both leaf and root tissues of rice (Mathan et al., 2020). Similarly, heat-stressed moth bean seedlings exhibited a notable accumulation of total sugars and proline across different genotypes (Harsh et al., 2016).

Heat stress not only impairs carbohydrate assimilation but also reduces sugar export from leaves (Julius et al., 2017). Marias et al. (2017) observed increased glucose and fructose levels in the leaves of Coffea arabica under heat stress. In maize, high temperatures decreased the 13C export rate from ear leaves, leading to enhanced growth of vegetative parts but reduced grain yield (Suwa et al., 2010). It is proposed that SWEET transporters efflux sucrose into the phloem apoplast, while SUT transporters import sucrose into the sieve element-companion cell complex (Osorio et al., 2014). During heat stress, leaves are the first tissues to sense and suffer from heat damage, and they utilise sugars to scavenge reactive oxygen species (ROS) and alleviate oxidative stress. Additionally, sugars serve as an energy source in high-respiration environments, helping maintain cell–water balance and membrane integrity through osmotic adjustments (Burke, 2007). Given the critical role of sugars in stress responses, the findings of this study, combined with previous observations, suggest that heat stress inhibits sugars export by SWEETs and/or its import by SUTs in heat-stressed leaves, resulting in sugars accumulation. Thereby, maintaining leaf sugar content through reducing source to sink export is considered a tolerance strategy (Kaushal et al., 2013). It is also worth noting that the decrease in sink demand due to growth limitation could also contribute to the accumulation of sugars in the source leaves under stressful conditions (Hummel et al., 2010; Lemoine et al., 2013).