Integrated Physiological and Transcriptomic Analysis Reveals Key Transcriptional Responses to Prolonged Heat Stress in Chinese Cabbage

1 Introduction

Heat stress is a major environmental factor that severely affects crop growth, productivity, and quality worldwide, sometimes even threatening plant survival [1,2,3]. High temperatures disrupt plant development by inducing cellular and tissue damage, resulting in water loss, membrane destabilization, and impaired photosynthetic efficiency [2,4]. As a result, sessile plants have evolved specialized strategies to cope with such environmental challenges [3].

To respond to heat stress, plants employ diverse defense mechanisms, including metabolite accumulation and activation of signaling pathways [2,5]. At the molecular level, heat stress induces the production of reactive oxygen species (ROS) and activates antioxidant enzymes, thereby altering metabolic processes [6,7]. In particular, increased activities of antioxidant enzymes such as superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) help suppress ROS accumulation, thereby minimizing cellular damage [6,8].

Recent studies have further highlighted the complexity of plant responses to heat stress, including mechanisms related to stress perception, signal transduction, redox regulation, and adaptive responses, which collectively contribute to plant adaptation under prolonged heat stress conditions [8,9,10,11].

Chinese cabbage (Brassica rapa L.), a cool-season leafy vegetable, grows optimally at 18–22°C and is predominantly cultivated in autumn [12]. However, exposure to high temperatures delays growth, induces leaf yellowing and wilting, impairs head formation, and increases susceptibility to diseases, collectively leading to a marked decline in both yield and quality [13]. In contrast, heat-tolerant varieties exhibit a greater increase in antioxidant enzyme activity compared to heat-sensitive cultivars, which helps maintain head formation under high-temperature conditions [12,13]. Therefore, understanding the mechanisms underlying heat tolerance is important for the development of resilient cultivars.

Nevertheless, most existing studies have primarily focused on short-term heat stress models that are effective for capturing early responses but may not fully reflect the prolonged heat stress conditions frequently encountered in agricultural environments [14,15,16]. Even when the stress period has been extended to several days, studies have predominantly focused on comparisons between contrasting genotypes (e.g., heat-tolerant vs. heat-sensitive) [13].

In particular, previous transcriptomic studies in Brassica rapa have largely examined heat stress responses within short time frames, typically ranging from several hours to 24 h [14,15,16]. While these studies have provided valuable insights into early stress-responsive pathways, they do not fully capture the dynamic transcriptional changes associated with prolonged heat exposure over multiple days.

RNA sequencing (RNA-seq) is a powerful high-throughput approach for transcriptome analysis, enabling the identification and quantification of gene expression under various environmental conditions [17,18]. It allows the detection of differentially expressed genes (DEGs) and provides insights into regulatory pathways associated with stress responses.

Transcription factors (TFs) play an important role in regulating gene expression and are involved in various physiological processes and stress responses. TFs regulate plant growth, development, and environmental adaptation by binding to cis-acting elements located in the promoters of target stress-responsive genes, thereby activating or repressing their expression [19,20,21,22]. The HSF (heat shock transcription factors) family is a key group of transcription factors that respond to heat stress and play a significant role in heat resistance [19]. Additionally, TF families such as NAC, MYB, WRKY, bZIP, and AP2/ERF are known to play crucial regulatory roles under various abiotic stresses, including drought, salinity, and heat [20,21,22]. While the heat shock factor (HSF) family has been well-studied, especially under short-term stress conditions [19], the regulatory roles of these other TF families under prolonged heat stress remain poorly understood.

In this study, we aimed to fill this knowledge gap by analyzing the transcriptomic and transcription factor responses of Chinese cabbage under prolonged heat stress conditions. By integrating physiological assessments and transcriptome profiling, we sought to provide insights into the temporal regulation of heat stress responses. These findings may contribute to the development of heat-tolerant Chinese cabbage cultivars.

2 Materials and Methods

2.1 Plant Materials and Growth Conditions

The National Institute of Horticultural and Herbal Science in Korea grew the Chinese cabbage cultivar ‘Chunkwang’ (Brassica rapa subsp. pekinensis) in a greenhouse under natural light conditions. Seeds of the cultivar were purchased from SAKATA KOREA (Seoul, Korea). After sowing, the seedlings were grown in a greenhouse for 6 days. On the 6th day, they were transferred to a growth chamber with light intensity at 350 μmol m−2 s−1, temperature at 25°C/20°C (16 h Light/8 h Dark), and humidity at 75%/60% (16 h L/8 h D). The seedlings were allowed to adapt for 5 days under long-day conditions, reaching the 11-day-old seedling stage with fully expanded first and second true leaves. Subsequently, half of the plants underwent heat treatment. This treatment was conducted in a growth chamber set to 40°C/35°C (16 h L/8 h D) and 75%/60% (16 h L/8 h D) humidity.

Samples were collected at 0, 1, 2, 4, 5, 6, 7, and 9 days after heat treatment by harvesting the first and second true leaves. The samples were immediately stored at −80°C.

For short-term heat treatment experiments, the seeds were sown in plug trays filled with 110 mL of commercial potting soil (Baroker potting mix, Seoul Bio Co., Ltd., Korea) per cell. After 14 days, the seedlings were transplanted into 450 mL pots and provided with sufficient water. The plants were acclimated in the greenhouse for 2 days, and then transferred to a growth chamber (Controlled Environment Rooms, Environmental Growth Chamber, USA) set at 340 μmol m−2 s−1 light intensity, constant temperature of 21°C (day/night), and a 16/8-h light/dark photoperiod for 5 days of adaptation. For heat treatment, half of the potted plants at the seven-leaf stage (23 days old as the vegetative stage) were transferred to a growth chamber set at 40°C. Leaf samples were collected at 1, 3, 6, 9, 12, and 15 h, taking the 3rd and 4th true leaves from the bottom. The collected leaf samples were immediately flash-frozen in liquid nitrogen and stored at −80°C for RNA extraction.

2.2 RNA Preparation and Transcriptome Sequencing Analysis

Total RNA was extracted using the TRIzol Reagent method (Invitrogen, USA) according to the manufacturer’s instructions [23]. Approximately 100 mg of leaf samples were pre-processed with 1 mL of TRIzol Reagent, followed by the addition of 100 μL of chloroform to disrupt the cell membranes. The subsequent steps were carried out using the Qiagen RNeasy Mini Kit, starting from the RNeasy Mini spin column purification.

To obtain pure RNA, the extracted RNA was treated with DNase I (Qiagen) to remove any DNA contamination. The quality and purity of the RNA were verified using a NanoDrop 2000 spectrophotometer and 1% agarose gel electrophoresis.

For transcriptome sequencing, four biological replicates were processed. Library preparation was conducted using the TruSeq Stranded Total RNA Library Prep Plant Kit (Illumina), ensuring depletion of ribosomal RNA and inclusion of coding and non-coding transcripts. The prepared libraries were sequenced on the Illumina platform, generating paired-end reads with a read length of 101 bp.

The resulting raw sequence data were quality-checked using FastQC (v0.11.9) [24] and subjected to adapter trimming and low-quality read removal using Trimmomatic (v0.39) [25]. High-quality reads were aligned to the Brassica rapa reference genome (GCF_000309985.2_CAAS_Brap_v3.01) using HISAT2 (v2.1.0) [26]. The resulting alignments were processed using SAMtools [27]. Gene/transcript abundance was calculated with StringTie (v2.1.3b) [28], and normalized expression values were represented as Transcripts Per Kilobase Million (TPM).

2.3 GO and KEGG Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to identify functional categories and pathways associated with DEGs. GO enrichment analysis was conducted for the three main categories: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), using the DAVID bioinformatics resources (https://david.ncifcrf.gov/).

KEGG pathway analysis was carried out using the KEGG Mapper tools (http://www.kegg.jp/kegg/pathway.html). DEGs were mapped to KEGG pathways, and enrichment was analyzed based on pathway-specific gene mapping and pathway completeness using the species-specific KEGG database (Brassica rapa reference genome: GCF_000309985.2). Statistical significance of enriched pathways was determined using Fisher’s exact test, with results adjusted for multiple testing using Bonferroni correction and false discovery rate (FDR) methods.

2.4 Protein Quantification and Analysis of Malondialdehyde (MDA) and H₂O₂ Contents

The same samples used for RNA extraction were processed for analysis. Each sample (100 mg) was homogenized with 1 mL of PBS (potassium phosphate buffer saline, pH 7.0) and centrifuged at 13,000 rpm for 30 min at 4°C. The supernatant was collected for further analysis. Protein quantification was conducted using the Bradford method [29]. MDA and H2O2 contents were measured using the EZ-Lipid Peroxidation (TBARS) Assay Kit and the EZ-Hydrogen Peroxide/Peroxidase Assay Kit (DoGenBio Co., Korea), respectively, following the manufacturer’s instructions.

2.5 Antioxidant Enzyme Activity Analysis

The activities of antioxidant enzymes, including ascorbate peroxidase (APX), catalase (CAT), peroxidase (POD), and superoxide dismutase (SOD), were analyzed using the following assay kits: Ascorbate Peroxidase Activity Assay Kit (Elabscience Biotechnology Inc., USA), EZ-Catalase Assay Kit, EZ-Hydrogen Peroxide/Peroxidase Assay Kit, and EZ-SOD Assay Kit (DoGenBio Co., Korea), respectively. The assays were performed according to the manufacturer’s protocols.

2.6 Quantitative RT-PCR Analysis

The synthesized cDNA was used for qPCR analysis using EvaGreen qPCR Master Mix (BIOFACT, Daejeon, South Korea) and a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Relative mRNA levels were determined by normalizing the PCR threshold cycle number of each target gene with that of the reference gene, Actin [30]. In the qPCR analysis, three technical repeats were measured for each biological replicate analyzed. The primers used for qPCR analyses are presented in Table S1.

3 Results

3.1 Growth of Chinese Cabbage under Heat Stress

In this study, we used the Chinese cabbage cultivar ‘Chunkwang’ (Brassica rapa subsp. pekinensis) to compare gene expression profiles of young seedlings under heat stress using RNA-seq analysis.

First, to investigate the effects of heat stress on the growth and physiological changes in 11-day-old Chinese cabbage seedlings, we compared the length of the first true leaf in the control group (O_group, 24°C) and the heat-treated group (H_group, 40°C) over a period of 9 days. Plants exposed to heat stress exhibited smaller leaf sizes compared to those grown under optimal conditions in the control group.

Under optimal temperature conditions (24°C), leaf growth exhibited a consistent progression, characterized by a gradual increase in the length of the first true leaf throughout the experimental period. In contrast, the group subjected to heat treatment at 40°C exhibited inhibited growth, as evidenced by a slower rate of leaf elongation compared to the control group, which resulted in relatively smaller leaf sizes from the initial stages of the treatment (Fig. 1a). Notably, at 9 days after heat treatment (DAH), the average leaf length of the heat-treated group was approximately 39% shorter than that of the control group (Fig. 1b).

Following a period of 5 days of heat treatment, the cotyledons of the treated plants displayed significant curling and initiated a yellowing process (Fig. 1a) [12]. These observations indicate stress-induced damage, and further analyses were conducted at 2 DAH (short-term heat stress) and 5 DAH (long-term heat stress) to investigate physiological changes under heat stress. These time points were selected to capture distinct stages of prolonged heat stress, with 2 DAH representing an early transitional phase beyond short-term responses (≤24 h) before visible phenotypic damage, and 5 DAH corresponding to the onset of visible stress symptoms (Fig. 1a).

Figure 1: Phenotypic changes and leaf length analysis in Chinese cabbage under heat stress: (a) Phenotypic changes in Chinese cabbage (Brassica rapa subsp. pekinensis) under heat stress were observed. 11-day-old plants in the control group (24°C, upper row) and the heat-treated group (40°C, lower row) were monitored daily over a 9-day period. On day 5 after heat treatment (DAH), noticeable yellowing of cotyledons was observed in the heat-treated group, highlighted with red arrows to emphasize the phenotypic differences caused by heat stress. Scale bar = 10 cm; (b) Changes in the length of the first true leaf were analyzed in the control and heat-treated groups over 9 days period after heat treatment. The graph illustrates changes in leaf length over time. Error bars represent SD (n ≥ 4). Asterisks indicate statistically significant differences between groups (***p < 0.001).

3.2 Oxidative Stress and Antioxidant Responses in Chinese Cabbage

ROS scavenging system and the activity of antioxidant enzymes were compared between the control and heat-treated groups at different time points.

SOD activity increased by approximately 18% in the heat-treated group compared to the control group at 2d and by 41% at 5d (Fig. 2a), while POD activity decreased by approximately 32% in the heat-treated group at 2d (Fig. 2b). In contrast, CAT activity showed no significant differences between the control and heat-treated groups at both 2d and 5d (Fig. 2c). Ascorbate peroxidase (APX) activity increased by 28% in the heat-treated group at 2d and by 73% at 5d (Fig. 2d). H2O2 concentration in the heat-treated group decreased by approximately 65% at 2d and by 61% at 5d, showing lower levels than in the control group (Fig. 2f). Malondialdehyde (MDA) content showed no significant differences at 2d, but increased by approximately 33% in the heat-treated group at 5d (Fig. 2e). These results indicate that 2-day and 5-day heat stress have distinct effects on plants.

Figure 2: Antioxidant enzyme activities and oxidative stress markers in the control and heat-treated groups at 2 and 5 days after heat treatment (DAH): (a) Superoxide dismutase (SOD) activity; (b) Peroxidase (POD) activity; (c) Catalase (CAT) activity; (d) Ascorbate peroxidase (APX) activity; (e) Malondialdehyde (MDA) content; (f) Hydrogen peroxide (H₂O₂) content. Bars represent the mean values ± standard error (SE) of three biological replicates. Asterisks indicate statistically significant differences between groups (*p < 0.05; **p < 0.01; ***p < 0.001).

3.3 Transcriptome Profiling of Chinese Cabbage under Heat Stress

To investigate the transcriptomic response of Chinese cabbage to heat stress, RNA-seq was performed on samples collected at 2 DAH and 5 DAH as well as from control plants grown under optimal conditions. The results showed that more than 55 million high-quality reads were generated per sample, with over 94% of the reads successfully mapped to the reference genome (Table 1).

The raw sequencing reads were processed and subjected to quality control using FastQC [24], and trimming was performed with Trimmomatic to remove low-quality bases and adapter sequences [25]. The trimmed reads were then aligned to the Brassica rapa reference genome (GCF_000309985.2_CAAS_Brap_v3.01) using HISAT2 [26]. The high mapping rates (94.74%–97.41%) and minimal levels of unmapped reads confirmed the reliability and quality of the data, ensuring suitability for downstream analyses. Transcript abundance was calculated using StringTie [28], and differential expression analysis was conducted using edgeR [31], identifying significant DEGs with |log2 fold change| ≥ 2 and p-value < 0.05. The high-quality reads produced in this study have been deposited in the NCBI BioProject database (GenBank accession number: PRJNA1184374).

To evaluate the overall transcriptomic changes induced by heat stress, multidimensional scaling (MDS) analysis was conducted based on the gene expression profiles of all samples. The MDS plot (Fig. 3) reveals clear clustering between the heat-treated (H_group) and control (O_group) samples. Control samples (O_group) exhibited a close clustering, reflecting a consistent pattern of gene expression under normal conditions. In contrast, the heat-treated samples (H_group) were distinctly separated from the control samples along the primary component axis. The x-axis (Component 1) accounted for 71.4% of the total variance, while the y-axis (Component 2) explained 20.7% of the variance, indicating transcriptomic changes induced by heat stress. Notably, the 2-day heat-treated sample (2-H) exhibited a clear separation from the 2-day control sample (2-O) in clustering. Furthermore, the separation between the 5-day heat-treated sample (5-H) and the 5-day control sample (5-O) was even more pronounced than that between the 2-H and 2-O samples, indicating that short-term and long-term heat stress have distinct effects on gene expression.

Figure 3: MDS (Multidimensional Scaling) Plot analysis: The MDS plot visualizes the transcriptomic differences between heat-treated (H_group) and control (O_group) samples. The x-axis (Component 1) explains 71.4% of the total variance, while the y-axis (Component 2) accounts for 20.7% of the variance. Control samples are blue dots (2 and 5 days), and heat-treated samples are red dots.

3.4 Analysis of Differentially Expressed Genes (DEGs) in Chinese Cabbage under Heat Stress

RNA-seq analysis identified a substantial number of DEGs at both 2 and 5 days after heat treatment, consistent with prior transcriptomic studies in Brassica species under abiotic stress [32,33]. At the 2-day time point (2-H vs. 2-O), a total of 4137 genes were up-regulated, while 3928 genes were down-regulated. At the 5-day time point (5-H vs. 5-O), 5738 genes were up-regulated, and 6048 genes were down-regulated (Fig. 4a). The increased number of DEGs at day 5 indicates an intensified transcriptional response under prolonged heat stress.

Subsequently, a Venn diagram illustrated the overlap and distinct sets of DEGs across the two time points. Among these, 2507 genes were consistently up-regulated, and 2133 were consistently down-regulated across both time points (Fig. 4b), suggesting the presence of shared heat-responsive gene sets that may contribute to sustained thermotolerance in Chinese cabbage.

Figure 4: Analysis of Differentially Expressed Genes (DEGs); (a) The number of differentially expressed genes (DEGs) between the heat-treated group (H_group) and the control group (O_group) is presented based on Fold Change (|FC| ≥ 2) and statistical significance (raw p < 0.05). The graph on the left represents data from day 2, while the graph on the right represents data from day 5. The yellow bars represent the number of up-regulated genes, while the blue bars indicate the number of down-regulated genes. (b) Venn diagrams showing the overlap of differentially expressed genes (DEGs) between 2-H vs. 2-O and 5-H vs. 5-O comparisons. The left panel represents up-regulated genes, and the right panel shows down-regulated genes. The overlapping regions indicate unigenes consistently expressed across both time points.

3.5 Functional Annotation of DEGs Responsive to Heat Stress

The functional categorization of up-regulated DEGs under heat stress exhibited distinct patterns at day 2 (2-H vs. 2-O) and day 5 (5-H vs. 5-O). Gene Ontology (GO) analysis of up-regulated genes at day 2 revealed that in the BP (Biological Process) category, RNA modification (GO:0009451) was the most enriched term, with 96 genes involved. In the CC (Cellular Component) category, mitochondrion (GO:0005739) showed the highest level of enrichment, encompassing 142 genes. In the MF (Molecular Function) category, protein binding (GO:0005515) was the most significantly enriched term, with 568 genes involved (Fig. 5a). At day 5, the enrichment pattern had shifted, with translation (GO:0006412) becoming the most enriched BP term, involving 163 genes. Mitochondrion (GO:0005739) remained the dominant CC term, increasing to 196 genes, and protein binding (GO:0005515) continued to be the most enriched MF category, now with 676 genes (Fig. 5a). These results suggest a transition from early regulatory processes to enhanced protein synthesis and energy metabolism during prolonged heat exposure.

For down-regulated genes, GO analysis indicated a consistent repression of signaling-related processes. On day 2, protein phosphorylation (GO:0006468) was the most suppressed BP term (182 genes), while the extracellular region (GO:0005576) and protein kinase activity (GO:0004672) were the leading CC and MF terms, respectively (Fig. 5b). By day 5, repression intensified, with protein phosphorylation (GO:0006468) expanding to 324 genes, membrane components (GO:0016020) emerging as the most affected CC category (1224 genes), and DNA-binding transcription factor activity (GO:0003700) dominating the MF category (267 genes) (Fig. 5b).

KEGG pathway analysis further highlighted these temporal shifts. Among up-regulated genes on day 2, the top enriched pathways were “Protein Processing in the Endoplasmic Reticulum” (81 genes), “Plant Hormone Signal Transduction” (71 genes), and “Spliceosome” (64 genes). On day 5, “Ribosome” (186 genes) became the most significantly enriched pathway, alongside sustained enrichment of “Protein Processing in the Endoplasmic Reticulum” (118 genes) and “Spliceosome” (117 genes) (Fig. 5c).

In contrast, down-regulated DEGs were predominantly enriched in pathways associated with energy metabolism and hormonal signaling. On day 2, the top suppressed pathways included “Plant Hormone Signal Transduction” (92 genes), “Ribosome” (81 genes), and “Oxidative Phosphorylation” (56 genes). By day 5, “Plant Hormone Signal Transduction” remained significantly repressed (146 genes), followed by “Photosynthesis” (72 genes) and “Plant-Pathogen Interaction” (67 genes) (Fig. 5d). Collectively, these findings indicate a coordinated shift from early stress perception and regulatory adjustments to metabolic reprogramming and energy conservation strategies under prolonged heat stress.

Figure 5: GO and KEGG analysis: (a) Gene Ontology (GO) analysis of up-regulated DEGs at 2 days (left panels) and 5 days (right panels) after heat stress treatment; (b) Gene Ontology (GO) analysis of down-regulated DEGs at 2 days (left panels) and 5 days (right panels) after heat stress treatment. The GO categories are divided into Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Enrichment significance is indicated by *p-value thresholds (***p < 0.001); (c) KEGG pathway enrichment analysis of up-regulated DEGs at 2 days (2-H vs. 2-O, left) and 5 days (5-H vs. 5-O, right) after heat stress; (d) KEGG pathway enrichment analysis of down-regulated DEGs at 2 days (2-H vs. 2-O, left) and 5 days (5-H vs. 5-O, right) after heat stress. The top 20 significantly enriched pathways are displayed, ranked by GeneRatio (x-axis). The size of each dot represents the number of significantly enriched genes (Number of SigGenes), while the color gradient indicates the significance level (p-value), with red denoting higher significance. Arrows show key pathways of heat stress.

3.6 Transcription Factors with Altered Expression Levels by Heat Stress in Chinese Cabbage

Transcription factors were classified based on the DEG analysis to evaluate their differential expression under heat stress conditions. A total of 24 TF families were identified as being upregulated by more than twofold on both day 2 and day 5, whereas 16 families were found to be downregulated (Fig. 6).

Figure 6: Up- or down-regulated expression of transcription factors identified by RNA-seq: The bar graph illustrates the distribution of transcription factors (TFs) with a fold change of 2 or greater (|FC| ≥ 2) under heat stress conditions at both day 2 and day 5. The graph shows the numbers of upregulated TFs (red bars) and downregulated TFs (blue bars) across TF families. The red bars represent the number of TFs whose expression is upregulated under heat stress, whereas the blue bars indicate the number of TFs that are downregulated during the stress period.

This is broad representation of TF families highlights the extensive transcriptional reprogramming in Chinese cabbage in response to prolonged heat stress. Previous studies have demonstrated that TF families such as HSF, NAC, bZIP, HD-ZIP, MYB, and bHLH are critically involved in plant responses to various abiotic stresses, including heat stress [19,20,21,22].

Next, we selected TFs with |FC| > 10 on both day 2 and day 5 (Table 2) to identify those that are highly sensitive to heat stress. Among these, four TFs belonging to the HSF family were identified, and the TF with the highest FC value on day 5 (LOC103862791) was designated as BrHSF B-1. For the NAC family, two TFs (LOC103829311 and LOC103875105) were identified and named BrNAC13 and BrNAC59, respectively. The bZIP family included one TF (LOC103853909), which was named BrbZIP43. For the HD-ZIP family, one TF (LOC103842010) was identified and designated as BrATHB-12. The C2H2 family included one TF (LOC103859816), which was named BrAZF2. For the COL family, two TFs (LOC103874939 and LOC103868662) were identified and named BrCOL10 and BrCOL8, respectively. The bHLH family included one TF (LOC103874966), which was designated as BrbHLH78. Lastly, the MYB family included one TF (LOC103855087), which was designated as BrMYB34. The selected TFs span diverse functional families, indicating that heat stress modulates a broad array of transcriptional regulators. In total, six up-regulated TFs and four down-regulated TFs were selected for further expression pattern analysis.

3.7 Validation and Expression Pattern Analysis of Transcription Factors through qRT-PCR

To validate the RNA-seq results and investigate the expression patterns of 10 selected transcription factors with significant fold change values on day 2 (2d) and day 5 (5d) in Chinese cabbage, time-course qPCR was conducted under both prolonged heat stress and short-term conditions. The analysis was performed using RNA extracted from the leaves of 11-day-old (seedling stage) and 21-day-old (vegetative stage) plants.

During prolonged exposure, spanning 9 days, transcription factors previously identified as up-regulated in the RNA-seq analysis exhibited substantial increases in transcript levels. These up-regulated genes, including BrHSF B-1, BrNAC13, BrbZIP43, BrATHB-12, BrAZF2, and BrCOL10, showed marked induction in the heat-treated group (H_group), with pronounced increases observed at 2d and 5d following heat treatment. In the H_group, expression levels increased by 10- to 135-fold at 2d compared to the control group (O_group), and further escalated to 45- to 775-fold at 5d relative to the O_group (Fig. 7a–f). Among these, BrHSF B-1 displayed the most pronounced increase, being markedly upregulated, reaching approximately 120-fold at 2d, 1013-fold at 5d, and peaking at 1273-fold at 7d relative to 0d in the H_group. Meanwhile, the O_group exhibited no noticeable changes in BrHSF B-1 expression at either 2d or 5d compared to 0d (Fig. 7a).

BrNAC13 expression in the O_group remained largely unchanged at 2d and increased moderately by 5-fold at 5d. In contrast, expression in the H_group was markedly upregulated, reaching approximately 45-fold at 2d and 346-fold at 5d, before peaking at 887-fold at 9d (Fig. 7b). BrbZIP43 followed a similar trend, showing little change in the O_group at 2d but increasing by 5-fold at 5d. In the H_group, expression was markedly upregulated, recording an increase of about 124-fold at 2d and 243-fold at 5d, and reaching a peak of 283-fold at 7d (Fig. 7c). BrATHB-12 expression in the O_group showed minimal variation at 2d and a modest 5-fold increase at 5d. Conversely, expression in the H_group was markedly upregulated, reaching approximately 32-fold at 2d and 607-fold at 5d (Fig. 7d). BrAZF2 in the O_group exhibited a slight increase of 1.3-fold at 2d and 3.7-fold at 5d, whereas in the H_group, expression was markedly upregulated, with levels surging by 14-fold at 2d and dramatically climbing to 743-fold at 5d (Fig. 7e). BrCOL10 expression in the O_group remained unchanged, whereas in the H_group, expression was markedly upregulated, increasing by about 72-fold at 2d and peaking at 236-fold at 5d (Fig. 7f).

In contrast, the down-regulated genes identified by RNA-seq, including BrNAC59, BrCOL8, BrbHLH78, and BrMYB34, generally exhibited lower expression levels in the H_group compared to the O_group. BrNAC59 expression increased modestly in both groups during the early stages, by approximately 1.7-fold at 2d and 7.2-fold at 5d in the H_group relative to 0d, with no clear difference between the two groups. However, divergence became apparent at 9d, when expression increased 73-fold in the O_group but only 43-fold in the H_group, indicating a 1.7-fold lower level in heat-treated plants (Fig. 7g). BrCOL8 expression declined by approximately 14 percent in the O_group and 55 percent in the H_group at 2d relative to 0d, resulting in a 1.9-fold greater reduction in the H_group. At 5d, expression in the O_group increased by about 2.3-fold, whereas in the H_group it decreased by 13-fold, leading to a 31-fold difference between the groups (Fig. 7h). BrbHLH78 expression was reduced by 31 percent in the O_group and by 78 percent in the H_group at 2d, corresponding to a 3.1-fold greater decrease in the H_group. Although the difference between the groups narrowed to approximately 1.17-fold at 5d, a substantial gap reappeared at 6d, with expression in the H_group approximately fourfold lower than in the O_group (Fig. 7i). BrMYB34 exhibited the most pronounced suppression, with expression in the H_group decreasing 37-fold relative to the O_group as early as 1d. Transcript levels remained low until 6d, after which a gradual increase was observed at 7d and 9d (Fig. 7j).

Figure 7: The expression level of 10 selected transcription factors (TFs) in Chinese cabbage (Brassica rapa subsp. pekinensis) under prolonged heat stress, validated by qRT-PCR. Seedlings were grown for 11 days under controlled conditions (25°C) following greenhouse acclimation and then subjected to high-temperature treatment (40°C). The first and second true leaves were collected from both control and heat-stressed plants at 0, 1, 2, 4, 5, 6, 7, and 9 days following the initiation of treatment. Among the TFs analyzed, six genes were upregulated in response to heat stress (Fig. 7a–f), whereas four genes were downregulated (Fig. 7g–j). Data are presented as means ± standard deviation of three biological replicates.

Additionally, to examine the expression profiles of transcription factors during the vegetative stage under short-term heat stress, qPCR analysis was performed using leaves from 21-day-old plants collected at 0, 1, 3, 6, and 12 h after heat treatment. The analysis revealed that all six transcription factors previously identified as upregulated under prolonged stress reached their peak expression levels as early as 1 h following heat exposure (Fig. S1a–f). In contrast, the four transcription factors classified as downregulated during prolonged heat stress exhibited variable expression patterns under short-term conditions. BrNAC59 and BrbHLH78 consistently showed lower expression in heat-treated plants compared to controls throughout the treatment period, with BrNAC59 declining by approximately 70% and BrbHLH78 by 93% at 12 h (Fig. S1g,i). Meanwhile, BrCOL8 and BrMYB34, which were suppressed under prolonged stress, exhibited transient upregulation at early time points following heat exposure (Fig. S1h,j). However, as the prolonged and short-term heat treatments were conducted at different developmental stages, direct comparisons of expression patterns between these conditions should be made with caution.

4 Discussion

Prolonged heat stress has become an increasingly significant constraint on global crop productivity under ongoing climate change [3,14]. Chinese cabbage (Brassica rapa subsp. pekinensis), a cool-season vegetable widely cultivated in East Asia, is particularly sensitive to sustained high-temperature conditions [12]. While previous studies have largely focused on short-term or acute heat stress models [14,15,16], such conditions may not fully represent the prolonged heat exposure commonly encountered in agricultural environments. Previous transcriptomic studies in Brassica rapa have predominantly focused on short-term heat stress responses, typically within several hours to 24 h [14,15,16]. While these studies have provided important insights into early stress-responsive mechanisms, they may not fully reflect the progressive transcriptional reprogramming that occurs under prolonged heat stress over multiple days. In this context, our study aimed to examine plant responses under extended heat stress by integrating physiological, biochemical, and transcriptomic analyses.

Physiological assessments revealed that prolonged heat exposure suppressed growth, as evidenced by reduced leaf elongation and visible stress symptoms such as chlorosis and cotyledon curling (Fig. 1). These responses were accompanied by progressive accumulation of reactive oxygen species (ROS) [7,34,35], as reflected by increased MDA content and the temporal upregulation of antioxidant enzymes, particularly APX and SOD (Fig. 2a,d,e). In contrast, catalase (CAT) activity remained relatively stable (Fig. 2c), suggesting a more constitutive role in ROS detoxification, whereas APX appeared to be more responsive under stress conditions (Fig. 2f). These findings are generally consistent with previous reports describing differential regulation of ROS-scavenging enzymes under prolonged abiotic stress conditions [34,35].

These transcriptional changes, particularly the upregulation of NAC and bZIP transcription factors, may be associated with the observed activation of antioxidant enzymes such as APX and SOD, indicating a possible relationship between transcriptional regulation and redox responses under prolonged heat stress. Overall, these biochemical responses may reflect coordinated changes in stress-responsive genes and transcription factors, suggesting a potential link between redox homeostasis and gene expression under prolonged heat stress [8,9,10,36].

Expanding upon these physiological and biochemical observations, our transcriptomic analysis revealed distinct shifts in gene expression dynamics. Multidimensional scaling (MDS) analysis demonstrated a clear separation between control and heat-treated samples, with divergence becoming more pronounced at day 5 (Fig. 3). This trend mirrored the cumulative physiological impairments and supported the idea that transcriptional reprogramming appears to become more pronounced as heat stress persists [32,33]. Rather than focusing on the magnitude of differential expression, these results suggest temporal changes in gene regulation between early and later stages of heat stress. Differential expression analysis further highlighted temporal shifts in gene regulation between the early and later stages of heat stress, suggesting early activation of stress-responsive pathways followed by broader transcriptional reprogramming as heat stress persisted [37].

Functional annotation of DEGs via GO and KEGG analyses provided insights into the biological processes associated with heat stress responses. In contrast, downregulated genes were primarily linked to signaling and regulatory functions. The suppression of “protein kinase activity” (GO:0004672) at day 2 and “DNA-binding transcription factor activity” (GO:0003700) at day 5 may reflect a progressive reduction in growth- and regulatory-related processes, potentially reallocating resources toward stress adaptation (Fig. 5b), thereby indicating a shift in resource allocation toward stress survival rather than growth. KEGG pathway analysis corroborated these findings, highlighting early activation of pathways such as “protein processing in the endoplasmic reticulum” and “hormone signaling,” followed by sustained engagement of “ribosome” and “protein folding” pathways at day 5 (Fig. 5c). These results suggest a shift from early stress signaling to sustained metabolic and protein homeostasis processes under prolonged heat stress. Meanwhile, the downregulation of “photosynthesis” and “hormone signaling” pathways further emphasized a shift towards energy conservation and stress prioritization (Fig. 5d) [38,39,40].

Complementing these findings, transcription factor analysis uncovered dynamic expression patterns of key regulators (Fig. 6). Heat stress prominently induced TF families such as HSF, NAC, bZIP, and HD-ZIP, known for their roles in abiotic stress tolerance [21,22,41,42,43,44], while repressing several TFs involved in growth and development [45]. Notably, these changes were more pronounced at day 5 than at day 2, indicating an intensification of transcriptional reprogramming as heat stress persisted. To validate the transcriptomic findings, we performed qRT-PCR analysis of six key transcription factors that were markedly upregulated under heat stress: BrHSF B-1, BrNAC13, BrbZIP43, BrATHB-12, BrAZF2, and BrCOL10. All six genes showed strong induction, with expression levels peaking at day 5, and notable upregulation was evident as early as 1 h after heat exposure (Fig. 7a–f and Fig. S1a–f). These patterns were consistent with RNA-seq data, supporting their potential roles in heat stress responses (Table 2). Among them, HSFs are known to regulate heat shock proteins (HSPs) and play central roles in heat stress signaling by activating HSPs and other protective genes [19,41,44]. This pattern is consistent with the activation of protein folding and stress response pathways under prolonged heat stress. Similarly, the expression patterns of BrNAC13 and BrbZIP43 suggest their involvement in transcriptional regulation under prolonged heat stress. NAC and bZIP transcription factors have been widely reported to play important roles in abiotic stress responses, including the regulation of ROS detoxification and hormone-mediated signaling pathways such as ABA signaling [45,46]. These findings are consistent with previous studies, supporting their potential roles in stress responses under prolonged heat stress. These functions suggest that the observed transcriptional changes may be associated with the modulation of antioxidant defense systems under prolonged heat stress conditions [8,9,10,36].

Conversely, four transcription factors—BrNAC59, BrCOL8, BrbHLH78, and BrMYB34—exhibited modest downregulation under prolonged heat stress (Fig. 7g–j). Although RNA-seq data indicated a decline in their expression, qRT-PCR revealed less pronounced changes, likely due to methodological differences: RNA-seq offers broad coverage but may underrepresent low-abundance transcripts, whereas qRT-PCR provides high sensitivity for individual targets [17]. Interestingly, BrCOL8 and BrMYB34 showed transient upregulation at early time points, suggesting potential roles in the initiation of heat response signaling, followed by transcriptional suppression as part of an energy-conserving strategy under prolonged stress [22,34]. It is important to note that these trends were observed under distinct developmental contexts for short-term and long-term treatments, warranting cautious interpretation when comparing dynamic expression profiles. Nonetheless, these complementary datasets offer valuable insight into the temporal coordination of stress-responsive gene networks in Chinese cabbage.

Moreover, promoter analysis of the heat-responsive TFs revealed the presence of multiple stress-related cis-elements, including motifs for HSF, NAC, bZIP, and other TF families (Fig. S2). This observation suggests that these genes may be involved in complex transcriptional regulatory networks under heat stress conditions [47]. While further functional validation is required, the identification of these regulatory motifs provides valuable insights into the regulatory landscape governing heat stress responses in Chinese cabbage.

Taken together, our results suggest that prolonged heat stress is associated with changes in antioxidant activity and gene expression in Chinese cabbage. These changes reflect a transition from early stress responses to sustained adaptive processes. The data indicate dynamic shifts in stress-responsive pathways, including activation of antioxidant systems and transcriptional regulators, followed by broader metabolic adjustments as heat exposure persists. These findings contribute to improving our understanding of heat stress responses in Brassica crops. Future studies focusing on functional validation of candidate genes will further clarify their roles in stress adaptation.

5 Conclusions

This study provides a comprehensive overview of the physiological, biochemical, and transcriptional responses of Chinese cabbage (Brassica rapa subsp. pekinensis) to prolonged heat stress. Sustained high-temperature exposure resulted in significant morphological damage, ROS accumulation, and activation of antioxidant defense systems. Transcriptome-wide analysis revealed dynamic transcriptional reprogramming over time, with notable shifts from RNA processing to protein synthesis pathways. We identified six heat-inducible transcription factors (BrHSF B-1, BrNAC13, BrbZIP43, BrATHB-12, BrAZF2, and BrCOL10) and four repressed TFs with potential regulatory roles in long-term heat adaptation. These findings offer novel insights into the molecular mechanisms underlying thermotolerance and provide promising targets for the development of heat-resilient Brassica cultivars. Future studies involving functional characterization of these transcription factors may help clarify their roles in heat stress response and crop improvement.