Heat stress response in Chinese cabbage (Brassica rapa L.) revealed by transcriptome and physiological analysis

Plant material and treatment

Two Chinese cabbage highly inbred lines (heat-tolerant ‘268’ and the heat-sensitive ‘334’) were provided by the Chinese Academy of Agricultural Sciences, Beijing, China. Five hundred full and uniform seeds were selected for germination, sowing, and seedling raising. They were grown at 25 ± 2 °C with a photoperiod of 16 h light/8 h dark, a 150 µmol m⁻² s⁻¹ light intensity, and a humidity level of about 54% to avoid dehydration. After 30 days of growth, the ‘268’ and ‘334’ lines were separated into two groups: one was treated at 40 °C/30 °C (light/dark), and the other was treated at 25 °C (control). On the 10th day of treatment, the heat-stressed plants were transferred to the control environment for 4 d.

At 0 (CK), 4 (HT-4), 8 (HT-8), and 10 (HT-10) d of heat stress at 40 °C, and at 4 (RC-4) d after recovery treatment at 25 °C, three Chinese cabbage plants with the same level of growth (three biological replicates) were selected to sample the functional leaves, namely, the third leaf of each seedling center that was fully expanded from the inside to the outside. The samples were immediately pre-cooled in liquid nitrogen and then stored in an ultra-low temperature refrigerator at −80 °C for later RNA-seq. The sample used for measuring physiological indexes was selected from the functional leaves of the same part, cut into pieces, and mixed for later use.

Physiological measurements

The differences in the physiological responses of the ‘268’ and ‘334’ lines under control and heat stress treatments were measured during CK, HT-4, HT-8, HT-10, and RC-4 periods. In accordance with the manufacturer’s instructions, the corresponding kits (Solarbio Life Sciences, Beijing, China) were selected to determine various physiological indicators. The Catalase (CAT) Activity Assay Kit (Cat#BC0200; Solarbio, Beijing, China), Superoxide Dismutase (SOD) Activity Assay Kit (Cat#BC0170; Solarbio, Beijing, China), Proline (Pro) Content Assay Kit (Cat#BC0290; Solarbio, Beijing, China), Malondialdehyde (MDA) Content Assay Kit (Cat#BC0020; Solarbio, Beijing, China), and Plant Soluble Sugar (SS) Content Assay Kit (Cat#BC0030; Solarbio, Beijing, China) were used for subsequent spectrophotometric measurements. Three biological replicates of each sample (n = 3) were assessed.

Total RNA extraction, library preparation and assembly, and analysis

The heat-tolerant ‘268’ line and heat-sensitive ‘334’ line were used for testing. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and 30 cDNA libraries (two lines × five periods × three replicates) were constructed. The purity, concentration, and integrity of the RNA samples were tested using Nanodrop2000 (Shimadzu, Japan) and agarose gel electrophoresis. A total amount of 1 μg RNA per sample was used as input material for the RNA sample preparations. The Next Ultra TM RNA library preparation kit (Illumina, New England Biolabs, Ipswich, MA, USA) was used to construct the library. PCR was performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index Primer. After product purification (AM Pure XP system), the library quality was evaluated using the Agilent Bioanalyzer 2,100 system. Following testing, the library preparations were sequenced on an Illumina platform and paired-end reads were generated. The raw RNA-seq data are available in the NCBI SRA database under the accession number PRJNA663233.

Raw reads in fastq format were first processed using in-house perl scripts. In this step, clean reads were obtained by removing those containing adapters and ploy-N, as well as low-quality reads, from the raw data. At the same time, the Q30, GC-content, and sequence duplication level of the clean data were calculated. The filtered data were compared with the B. rapa V3.0 reference genome (http://brassicadb.cn/#/Download/) using HISAT2 (Kim, Langmead & Salzberg, 2015). String Tie (Pertea et al., 2015) was then used to assemble and quantify results. The Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) value of each sample was calculated using Cufflinks V2.2.1 (Trapnell et al., 2012). To evaluate the replicate reproducibility, Pearson correlation analyses were performed in R language.

Identification and functional annotation of differentially expressed genes

Power analysis calculations were performed on transcriptome samples using the RNASeqPower package (https://doi.org/doi:10.18129/B9.bioc.RNASeqPower) with the following parameter settings: Depth = 20, Effect = 2, and α = 0.05 (Conesa et al., 2016). The DESeq R package (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) was used to analyze the RNA-seq results (Love, Huber & Anders, 2014). Genes with FDR < 0.01 and |log₂(fold change)|> 2 were identified as differentially expressed in reference to previous methods used to adjust the false discovery rate (FDR) (Benjamini & Hochberg, 1995).

Based on the B. rapa V3.0 genome, all DEGs were annotated through the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to determine their biological functions and participation pathways. GOseq R was used to annotate GO (Young et al., 2010). Gene functions were annotated based on the Swiss-Prot, EuKaryotic Orthologous Groups (KOG)/Clusters of Orthologous Groups (COG), KEGG Orthology, Pfam, and Non-redundant (Nr) NCBI databases.

Gene temporal expression analysis and hub gene analysis

Through gene co-expression trend analysis, the expression patterns of Chinese cabbage genes under heat stress were studied. We compared DEGs in S268-CK_vs_S334-CK, S268-HT-4_vs_S334-HT-4, S268-HT-8_vs_S334-HT-8, S268-HT-10_vs_S334-HT-10, and S268-RC-4_vs_S334-RC-4 in order to cluster the change patterns of unigenes through the BMKCloud website (www.biocloud.net). DEGs were clustered into 15 expression profiles. Clustered profiles of DEGs with p < 0.01 were considered to be significantly different from the reference set for each sample.

Hub genes are genes with the most junctions in each cluster. The B. rapa reference genome was analyzed against each cluster’s DEGs to predict the interaction network of these DEGs. Protein interactions were obtained from the STRING database (http://string-db.org/). kME values were used to indicate height. Pearson correlation coefficients between expression levels and module eigengenes were calculated and genetic networks were mapped using Cytoscape software (Chin et al., 2014).

Fluorescent quantitative real-time PCR

RNA samples were reverse transcribed to cDNA using EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix Kit (TransGen, Beijing, China) for qPCR. The relative expression level of the DEGs were detected through quantitative real-time PCR using Taq Pro Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) and CFX-96 Real-time System (BIORAD, Hercules, CA, USA). Actin was used as an internal control, and the 2^−ΔΔCT method was used to calculate relative gene expression levels (Livak & Schmittgen, 2001). The primer sequences used in the qPCR analysis are shown in Table S13.

Statistical analyses

GraphPad Prism 8 (San Diego, CA, USA) and Microsoft Office Excel 2010 software (Redmond, WA, USA) packages were used to analyze the data collected in the physiological experiments and by qRT-PCR. IBM SPSS 25.0 (Armonk, NY, USA) tests were used to evaluate the significance of the differences. Tukey’s post-hoc test was used for mean comparisons. p < 0.05 was considered statistically significant.

Morphological and physiological biochemical changes

The ‘268’ line had a deeper leaf color, more compact phenotype, and thicker stem than the ‘334’ line (Figs. 1A and 1D). The moisture content (MC) of both lines showed a decreasing trend under heat stress, but the ‘268’ line had a more significant decline, from 91.8% to 82.7% (Fig. 1G). Under heat stress and recovery, the ‘268’ line showed continuous accumulation of Pro and SS content. However, the Pro content of the ‘334’ line changed little across different treatment stages, and the SS content increased first and then decreased (Figs. 2A and 2B). The trend of the MDA content in both lines was very complicated, as reflected by the accumulation level and the time taken to reach the extreme value (Fig. 2C). The SOD activity of the ‘268’ and ‘334’ lines displayed changes similar to that of the Pro content. The difference was that the SOD activity of the ‘268’ line decreased slightly during HT-8 (Fig. 2D). Under heat stress and recovery, the CAT activity of the ‘334’ line increased and was always greater than the activity of the ‘268’ line, but it declined rapidly during RC-4, and the activity level was lower than that of the ‘268’ line (Fig. 2E). Heat stress caused significant damage to the ‘334’ line phenotype, the petioles of all tested plants were too long, and the leaf width became narrow and brittle. After 10 days of heat stress at 40 °C, leaf edge scorching was obvious (Figs. 1E and 1F), but the change in the ‘268’ line was relatively small (Figs. 1B and 1C). This suggests that heat stress may be less damaging to heat-tolerant plants than heat-sensitive plants.

Transcriptome quality control and DEG analysis

The sequencing results for 30 samples were analyzed, producing 1,522,694,456 raw reads and 211.72 Gb of clean data. Each sample had a minimum of 5.79 Gb of clean data, a Q30 percentage of 94.18% or more, and a comparison efficiency ranging from 86.85% to 92.66%. After the removal of rRNA, an average of 90.52% reads were mapped to the reference database (http://brassicadb.cn) (Table 1). The mapped files were then processed by HISAT2 and String Tie to generate a consensus transcriptome assembly containing 44,161 known and 3,529 novel genes.

To ensure sufficient depth and biological replication of the RNA-seq data, we used RNASeqPower to calculate the statistical power, which was greater than 0.92 for each group of samples (Table S1). Pearson’s correlation coefficients were used to test for biologically repeated correlations between samples. The generated cluster heatmap was used to observe the overall transcriptome correlation of the ‘268’ and ‘334’ lines under different treatments (Fig. S1). The results showed that the three biological replicates from each time period and the transcriptome data both exhibited good correlation.

Gene expression was analyzed based on FPKM, the ‘334’ line under the same treatment was used as a control, and the transcriptional results of the ‘268’ line were analyzed to identify DEGs. In different treatment groups, S268-CK_vs_S334-CK and S268-HT-4_vs_S334-HT-4 had the lowest number of DEGs, while S268-HT-8_VS_S334-HT-8 had the most significant difference (Fig. 3A). Comparisons of these five datasets showed that 118 genes overlapped among HT-4, HT-8, HT-10, and RC-4. Among these, 73 genes were upregulated and 44 genes were downregulated. These genes may play crucial roles in the heat-tolerance of Chinese cabbage (Figs. 3B–3D; Tables S2 and S3).

GO and KEGG analyses of DEGs

We conducted GO and KEGG pathway analyses of the DEGs. A total of 2,815 GO terms were assigned to the 3,360 DEGs that responded to heat stress (Fig. 4A). Pathways with higher enrichment were screened in KEGG and formed five main categories and 39 subcategories (Fig. 4B). The most abundant category was ‘metabolism’, distributed over a large number of subcategories, followed by ‘genetic information processing’. In ‘metabolism’, the most abundant subcategories were ‘biosynthesis of amino acids’, ‘carbon metabolism’, and ‘starch and sucrose metabolism’ with enrichment levels of 50, 47, and 29 DEGs, respectively. The richest subcategory in ‘genetic information processing’ was ‘ribosome’ with 78 (10.96%) DEGs. ‘Plant hormone signal transduction’ was also an enriched subcategory with 52 (7.30%) DEGs.

Dynamic gene expression pattern analyses during heat stress

The K-means method was used to analyze the dynamic changes of gene expression in the ‘268’ line, and all DEGs were classified into 15 clusters with similar expression profiles (Fig. 5). In Clusters 1, 3, and 13, the expression level of enriched DEGs showed a gradual induction state under heat stress, and a decrease during RC-4. The expression trends of DEGs in the enriched Clusters 2, 4, and 11 continued to decrease under heat stress and increased during RC-4. GO and KEGG analyses indicated that ‘carbon metabolism’, ‘biosynthesis of amino acids’, ‘plant hormone signal transduction’, ‘purine metabolism’, ‘glutathione metabolism’, and ‘plant–pathogen interaction’ were jointly enriched in Clusters 1, 3, and 13 (Figs. 6A–6C; Table S4). Many genes in Clusters 2, 4 and 11 were classified together into ‘photosynthesis’, ‘photosynthesis-antenna proteins’, ‘ribosome’, ‘porphyrin and chlorophyll metabolism’, ‘carbon metabolism’, ‘plant hormone signal transduction’, and ‘biosynthesis of amino acids’. It should be noted that, in Cluster 4, the ‘photosynthesis-antenna proteins’ (nine genes) and ‘photosynthesis’ (seven genes) pathways were significantly enriched. (Figs. 6D–6F; Table S5). The results suggest that these genes may be closely related to heat stress.

Clusters 8 and 15 showed increased gene expression during heat stress and recovery compared with the S268-CK with the greatest differential expressions both occurring during HT-8 (Fig. 5). The KEGG enrichment pathway of Cluster 8 indicated extensive involvement of these genes in the synthesis and degradation of amino acids (glutamate, alanine, and aspartate), the degradation and metabolism of fatty acids, and the biosynthesis of various secondary metabolites (Table S6). Notably, Cluster 8 underwent a decrease in gene expression levels during RC-4 with little difference when compared with the S268-CK. Combined with KEGG annotations, Cluster 15 was enriched for a large number of genes in ‘starch and sucrose metabolism’ (ko00500), ‘glycolysis/gluconeogenesis’ (ko00010), ‘Galactose metabolism’ (ko00052), ‘fructose and mannose metabolism’ (ko00051), and ‘amino sugar and nucleotide sugar metabolism’ (ko00520) (Table S6), including the BraA06g003260.3C. gene (Hexokinase-Like 1, HKL1), BraA01g028060.3C. gene (TPS10), BraA03g011780.3C. gene (BGLU4), BraA05g040820.3C. gene (SS2), BraA01g001150.3C. gene (PCKA), and BraA10g021910.3C. gene (At5g18200). The HKL1, TPS10, BGLU4, and SS2 genes play important roles in starch and sucrose metabolism. Additionally, many sugar metabolism pathways were also enriched in Cluster 1. This indicates that the expression of genes related to sugar metabolism in the heat-tolerant ‘268’ line was generally upregulated under heat stress (Table S4).

The top 10 hub genes of Cluster 15 were screened with the DEG interaction network for analysis (Figs. 7A and 7B; Table S7). The results showed that five hub genes were associated with ABA levels: BraA03g054960.3C. gene (TZF3), BraA02g000470.3C. gene (RPT6), BraA03g052020.3C. gene (NBR1), BraA06g035430.3C. gene (STUBL1), and BraA03g008710.3C. gene (MIEL1). An interaction network analysis was also performed on the downregulated Cluster 4 (Fig. 7D). We found that all 10 hub genes clustered in the ‘photosynthesis’ (ko00195) and ‘photosynthesis-antennal protein’ (ko00196) pathways (Table 2). Six of the seven genes associated with the light-harvesting chlorophyll a/b-binding (LHC) gene family belong to the LHCB gene subfamily (LHCB3, LHCB5, LHCB1.4, LHCB1.3, and LHCB2.4) and one belongs to the LHCA gene subfamily (LHCA4). Under heat stress and recovery, these genes showed a downward trend in both lines, but the expression levels in the ‘268’ line were always higher than those in the ‘334’ line during the same period. This indicates that the photosystem may be destroyed under heat stress with a significant reduction in photosynthetic efficiency, which is consistent with the observed plant phenotype.

The interaction networks of DEGs were screened in Clusters 1 and 2, and the top 10 hub genes for both selected clusters were associated with ribosomal proteins (Fig. 7C; Tables S8 and S9). The difference was that genes related to 80S ribosomal proteins were significantly enriched in Cluster 1. In contrast, Cluster 2 was enriched for many genes related to 70S ribosomal protein, which is widely involved in the composition of chloroplast ribosomes (Ahmed, Shi & Bhushan, 2017). Under heat stress, the level of the chloroplast ribosome-encoded protein was weakened, and the expression of some light energy capture factors was suppressed, which may have led to a decrease in photosynthesis intensity.

Identification of ROS-scavenger and heat shock transcription factor genes

To elucidate the molecular mechanisms of the ‘268’ and ‘334’ lines under heat stress, 18 enzymes, five HSFs, and 14 HSP-related genes were identified in all DEGs (Fig. 8; Tables S10–S12). The BraA01g001360.3C. gene (Peroxidase 50, Prx50), BraA10g028260.3C. gene (Peroxidase 52, Prx52), BraA10g029420.3C. gene (Peroxidase 54, Prx54), and BraA10g006200.3C. gene (Phospholipid hydroperoxide glutathione peroxidase 1, GPX1) showed significantly higher expression levels in the ‘268’ line than in ‘334’ line. The expression levels of the BraA09g062610.3C. gene (Superoxide dismutase, SOD1) and BraA04g020160.3C. gene (Superoxide dismutase 2, SOD2) in the ‘268’ line increased continuously compared with the ‘334’ line, but the levels decreased during RC-4, indicating that the two genes were more responsive to heat stress (Fig. 8, Table S10). During treatment, the BraA04g028620.3C. gene (Heat stress transcription factor A-1, HSFA1) was expressed only in the ‘268’ line, and the expression level gradually increased with the treatment time and then returned to the CK level during RC-4. Activation of HSPs reduces the likelihood of protein denaturation and improves the stress capacity of an organism. Two new genes, Brassica_rapa_newGene_7058 and Brassica_rapa_newGene_7942, were identified. Before treatment, the expression levels of the two genes were similar in the ‘268’ and ‘334’ lines, but the expression levels of the ‘268’ line increased significantly under heat stress and suddenly decreased during RC-4. In contrast, the gene expression in the ‘334’ line didn’t change under different treatment. We speculated that these two genes likely play active roles in heat stress in the heat-tolerant ‘268’ line (Fig. 8, Tables S11 and S12).

The DEGs from the transcriptome analysis detected by qRT-PCR

To confirm the results of the RNA-seq analysis, eight DEGs were verified using qRTPCR in B. rapa, including three genes related to ABA levels (MIEL1, STUBL1, and TZF3), three genes related to sugar metabolism (HKL1, TPS10, and At5g18200), and two genes related to ROS-scavenger activity (Prx54 and SOD2) (Fig. 9). The RNA-seq analysis and the qRT-PCR results were consistent across time periods, indicating the reliability of transcriptome sequencing. The expression level of ROS scavenger-related genes in the ‘268’ line was higher than that of the ‘334’ line, which was the same as the RNA-seq data, suggesting that these genes may be involved in plant stress tolerance.