Transcriptome analysis and differential gene expression profiling of wucai (Brassica campestris L.) in response to cold stress
BMC Genomics volume 23, Article number: 137 (2022)
Wucai suffers from low temperature during the growth period, resulting in a decline in yield and poor quality. But the molecular mechanisms of cold tolerance in wucai are still unclear.
According to the phenotypes and physiological indexes, we screened out the cold-tolerant genotype “W18” (named CT) and cold-sensitive genotype “Sw-1” (named CS) in six wucai genotypes. We performed transcriptomic analysis using seedling leaves after 24 h of cold treatment. A total of 3536 and 3887 differentially expressed genes (DEGs) were identified between the low temperature (LT) and control (NT) comparative transcriptome in CT and CS, respectively, with 1690 DEGs specific to CT. The gene ontology (GO) analysis showed that the response to cadmium ion (GO:0,046,686), response to jasmonic acid (GO:0,009,753), and response to wounding (GO:0,009,611) were enriched in CT (LT vs NT). The DEGs were enriched in starch and sucrose metabolism and glutathione metabolism in both groups, and α-linolenic acid metabolism was enriched only in CT (LT vs NT). DEGs in these processes, including glutathione S-transferases (GSTs), 13S lipoxygenase (LOX), and jasmonate ZIM-domain (JAZ), as well as transcription factors (TFs), such as the ethylene-responsive transcription factor 53 (ERF53), basic helix-loop-helix 92 (bHLH92), WRKY53, and WRKY54.We hypothesize that these genes play important roles in the response to cold stress in this species.
Our data for wucai is consistent with previous studies that suggest starch and sucrose metabolism increased the content of osmotic substances, and the glutathione metabolism pathway enhance the active oxygen scavenging. These two pathways may participated in response to cold stress. In addition, the activation of α-linolenic acid metabolism may promote the synthesis of methyl jasmonate (MeJA), which might also play a role in the cold tolerance of wucai.
Low temperature is one of the abiotic stresses that directly affect plant growth and development . It also affects the geographical distribution of plants and limits their growth . In response to adverse environmental conditions, plants have evolved a series of physiological and biochemical mechanisms. In plants, there are many receptors for different stress signals, which are involved in various stress responses and form a complex response and regulation network in response to stress . When plants are subjected to low temperature stress, a large number of soluble substances accumulate in their tissues to improve cold resistance. For example, starch will be hydrolyzed into simple sugars and other derived metabolites to increase the osmotic pressure of cells to prevent freezing . Plants accumulate a large amount of reactive oxygen species (ROS) under stress. ROS mainly include superoxide anions (O2·−), hydrogen peroxide (H2O2), hydroxyl radicals (·OH), and singlet oxygen (O21) . In order to maintain the intracellular ROS balance at a harmless concentration, plants have evolved a series of enzymatic and non-enzymatic ROS-scavenging mechanisms . Under stress, plants can regulate the activity of these enzymes or antioxidants through a series of physiological and biochemical mechanisms . In this way, the excessive accumulation of ROS in plants can be reduced and the stress resistance of plants can be improved.
Plant hormones such as cytokinin, abscisic acid, ethylene and jasmonic acid (JA) also play an important direct or indirect role in plant response to abiotic stress. Jasmonic acid and its derivatives, such as MeJA and cis-jasmone, are collectively referred to as jasmonates (JAs). They are fatty acids derived from cyclopentanone. They belong to the oxidized lipid family and are collectively referred to as oxidized lipids . In addition to responding to biological stress, they can also regulate the expression of many genes in response to abiotic stress, such as low temperature and salinity . Studies have shown that MeJA can activate antioxidant metabolic pathways and defense mechanisms in various crops and enhance cold tolerance . In addition, some members of ERF, bHLH, MYB, DREB, and WRKY transcription factor families, such as ERF53 , ICE1, DREB1A , MYB4 , and WRKY19 , are also involved in regulating the expression of cold-stress response genes.
High-throughput RNA-sequencing (RNA-seq) can precisely measure the transcription level and provide gene sequence information at the same time. It has high efficiency, high cost performance with high reliability, and is widely used in plant transcriptome analysis, especially in non-model plants that lack genome sequencing data . Over the last decade, RNA-seq has been used in many plant species to investigate plant responses to cold stress, including winter rapeseed , strawberry , and cotton , and has been confirmed as a powerful tool for plant genetics research .
Wucai (Brassica campestris L. ssp. chinensis var. rosularis Tsen) is a genotype of Chinese cabbage and is native to China . It is one of the main vegetables cultivated in the Yangtze–Huaihe River basin . Wucai is a semi-cold-tolerant vegetable, and the optimum growth temperature is 10 °C-20°C. It is a popular vegetable in autumn and winter for its nutritional value and good taste. However, different genotypes of the same species have different tolerance levels to cold stress . With the expansion of the planting area, it is an important task to select cold-tolerant genotypes of wucai.
The response of different genotypes to cold varies from tolerant to extremely sensitive. On the basis of previous studies, we screened six genotypes of wucai for their cold tolerance at the seedling stage . Through a study of the morphological, physiological, and biochemical properties of wucai under both control and cold stress conditions, cold-tolerant and sensitive genotypes were selected. We further conducted complete transcriptome sequencing on the cold-tolerant and sensitive wucai genotypes, providing valuable insights into the molecular mechanisms of this variety. The results lay a foundation for identifying genes with the potential to improve the cold tolerance of wucai, as well for developing cold-tolerant wucai varieties in future breeding programs.
Phenotypic and physiological differences between various genotypes under cold stress
The morphological characteristics of seedlings of six wucai genotypes were determined to select cold-tolerant and cold-sensitive genotypes (Fig. 1A). Based on measurements of the MDA content in leaves under stress and normal conditions, changes in the leaf lipid hydrogen peroxide production rate were studied. After cold treatment, the MDA content in Sw-1 reached the highest level, whereas that in W18 was at the lowest level (Fig. 1B). Relative electrolyte conductivity (REC) in Sw-3 was the highest, followed by Ta2, Sw-1, ws-2, W18 and ws-1 (Fig. 1C).
The performance index on an absorption basis (PIabs) of W18 was the highest and was significantly higher than that of Sw-1 (Fig. 1D). Vj increased the most in Sw-1, followed by Ta2, Sw-3, ws-2, ws-1 and W18 (Fig. 1E). These indexes indicated that W18 is a cold-tolerant wucai genotype and Sw-1 is a cold-sensitive genotype.
Differential responses of two wucai genotypes to cold stress
The phenotypes and physiological changes of two genotypes were significantly different under low temperature stress (Fig. 2). The REC and MDA content in W18 and Sw-1 were significantly increased under cold stress. However, the increase in W18 was smaller than that in Sw-1 (Fig. 2B, C). Under cold stress, the H2O2 content and O2·− generation rate of the two genotypes were significantly increased compared to those under NT. In addition, the H2O2 content and O2·− generation rate of W18 were lower than those of Sw-1 (Fig. 2D, E). At low temperature, the total antioxidant capacity (T-AOC) of both genotypes increased, but the T-AOC of W18 was higher than that of Sw-1 (Fig. 2F). Compared with the control condition, the PIabs of Sw-1 declined, while the PIabs of W18 had no significant change (Fig. 2G). These results indicated that W18 was more tolerant than Sw-1.
Mapping and quantitative assessment of Illumina sequences
The two comparison groups were named CT (LT vs NT) and CS (LT vs NT). Twelve cDNA libraries (three per treatment group) were constructed from total RNA. As indicated in Table S1, among the million raw reads obtained from the libraries, approximately 580.21 million clean reads were identified, ranging from 46.63 million to 49.78 million reads per library. More than 92% of the clean reads had quality scores higher than the Q30 level (an error probability of 0.1%), and the GC content was approximately 48% in each sample. The qualitied clean reads of each library were mapped to the B. rapa reference genome, and the mapping rates in each library exceeded 89% (Table S2). All of the RNA-seq data in this article have been deposited in the NCBI Sequence Read Archive (SRA) database and are accessible under accession number PRJNA735896.
Different transcriptome profiles of CT and CS under LT and NT conditions
After calculating the gene expression, we found that samples of CT and CS in NT had fewer genes than the samples in LT at the expression level of FPKM (the expected number of fragments per kilobase of transcripts per million mapped fragments) > 10, and CS was less than CT in LT (Fig. S1A). The results of the expression density distribution were similar (Fig. S1B). The relationship among the samples was analyzed by principle component analysis (PCA). The results showed that samples in the same treatment in each genotype clustered together, and there was significant distinction between the samples of the two genotypes and the different treatments (Fig. S1C).
DEG identification and analysis
There were 3536 DEGs in CT (LT vs NT) and 3887 DEGs in CS (LT vs NT). Although there were more DEGs in CS (LT vs NT) than in CT (LT vs NT), there were more downregulated DEGs in CS (LT vs NT) than in CT (LT vs NT). A total of 1846 DEGs (852 up-regulated DEGs and 994 down-regulated DEGs) were shared between CT (LT vs NT) and CS (LT vs NT), and 1690 DEGs were specific to CT (LT vs NT) (Fig. 3A). All DEGs were hierarchically clustered with respect to gene expression patterns and were evaluated using log10RPKMs (reads per kilobase per million mapped reads) of the two groups (Fig. 3B and Fig S2).
Although different with respect to cold tolerance, the two groups shared common genes that regulated cold tolerance, including dehydration-responsive element-binding proteins (DREBs), cold-regulated protein (COR15A), GSTs, and beta-glucosidase (BGLU) (Table S3). Several transcription factors, such as bHLH61, WRKY25, and MYB29 were also identified in both groups.
GO enrichment analysis of DEGs
Based on the corrected P-values, we selected the 10 most enriched GO terms in each category (Fig. 4, Table S4 and Table S5). In the biological process category, response to salt stress (GO:0,009,651), response to cold (GO:0,009,409), and defense response (GO:0,006,952) were common in both groups, and response to cadmium ion (GO:0,046,686), response to jasmonic acid (GO:0,009,753), and response to wounding (GO:0,009,611) were unique to CT (LT vs NT). In the cellular component, common DEGs were significantly enriched in the nucleus (GO:0,005,634), the integral component of membrane (GO:0,016,021), and the plasma membrane (GO:0,005,886), whereas mitochondrion (GO:0,005,739) and endoplasmic reticulum (GO:0,005,783) were enriched only in CT (LT vs NT). Metal ion binding (GO:0,046,872), DNA-binding transcription factor activity (GO:0,003,700), ATP binding (GO:0,005,524), and DNA binding (GO:0,003,677) were the top four in the molecular function category, and protein dimerization activity (GO:0,046,983) and iron ion binding (GO:0,005,506) were enriched only in CT (LT vs NT).
KEGG pathway enrichment analysis of DEGs
To understand the function of DEGs, we mapped the DEGs to the reference specification path in the KEGG database. The core genes were aligned with the KEGG database and were assigned to 20 KEGG pathways. KEGG pathway enrichment results showed that the significantly enriched KEGG pathways in both groups were starch and sucrose metabolism (brp00500) and glutathione metabolism (brp00480) (Fig. 5A, B). Based on a comparison of the KEGG pathways of CT (LT vs NT) and CS (LT vs NT), we were interested in pathways with higher enrichment in CT but lower enrichment in CS. One such way was α-linolenic acid metabolism.
DEGs involved in common KEGG pathways in CT (LT vs NT) and CS (LT vs NT)
The starch and sucrose metabolism network changed significantly in the acclimation to cold stress of both groups (Fig. 6A). We identified 31 DEGs in CT (LT vs NT) and 36 DEGs in CS (LT vs NT), with 23 and 26 up-regulated DEGs, respectively. The expressions of beta-amylase 2 (BAM2), sucrose-phosphate synthase1 (SPS1), BGLUs, and BraA01g006680.3C were the most significant in CT (LT vs NT) (Fig. 6B); these genes are mainly involved in the transformation of fructose, sucrose, and glucose and the hydrolysis of starch to maltose and finally to glucose. The starch, fructose, sucrose, and glucose contents of CT and CS increased significantly under cold stress, and there was no significant difference between the two groups (Fig. 6C-F).
In glutathione metabolism, 22 DEGs in CT (LT vs NT) and 21 DEGs in CS (LT vs NT) were identified, with 18 and 17 up-regulated genes, respectively (Fig. 7A, B). Ascorbate peroxidase 1 (APX1), phospholipid hydroperoxide glutathione peroxidase 6 (GPX6), and most of the GST gene family members were up-regulated. Compared with NT, the contents of reduced glutathione (GSH), oxidized glutathione (GSSG), and total glutathione (GSH + GSSG) in CT and CS were significantly increased under cold stress. The GSH content was increased by 69.86% and 55.54%, respectively, and the GSSG content was increased by 20.66% and 24.89%, respectively. The (GSH + GSSG) content increased in the two groups by 61.94% and 50.73%, respectively. The glutathione content and GSH/GSSG ratio under cold stress were higher than those of the control (Fig. 7C-F).
DEGs involved in the KEGG pathway specific to CT (LT vs NT)
KEGG enrichment analysis showed that the α-linolenic acid metabolism pathway was enriched in the CT group. Based on the metabolism pathway and heat map analysis (Fig. 8A, C), we found that 4-coumarate–CoA ligase (4CLL), allene oxide cyclase (AOC), and 12-oxophytodienoate reductase (OPR) were up-regulated in both groups, but the key gene, LOX, was up-regulated only in CT. JA is catabolized by JA carboxyl methyltransferase (JMT) to form its volatile counterpart MeJA. JMT in CS was down-regulated. The above GO enrichment analysis showed that the response to jasmonic acid was enriched only in CT; therefore, we studied the jasmonic acid signal transduction pathway (Fig. 8B), and the results showed that members of JAZ (TIFYs), an important gene family in this pathway, were up-regulated in CT (LT vs NT) and down-regulated in CS (LT vs NT) (Fig. 8C). Compared with NT, the MeJA content in CT under LT increased by 60.33%, while it was down-regulated by 9.38% in CS (Fig. 8D).
Identification of transcription factors in response to cold stress
The RNA-seq results showed that a total of 206 and 208 DEGs were annotated as TFs in CT (LT vs NT) and CS (LT vs NT), and were classified into 47 and 48 families, respectively. The top 10 TF families are shown in Fig. 9A. Among the recognized TF families, WRKY and bHLH represented the most abundant category, followed by MYB and ERF in CT (LT vs NT) and CS (LT vs NT). The trend of the fold change of common TFs between the two groups was similar. In the TFs unique to CT (LT vs NT), bHLH and WRKY were the top two TF families; most of the bHLHs were down-regulated and the WRKYs were up-regulated (Fig. 9B, C).
Validation of gene expression patterns by qRT-PCR
To validate the data obtained from RNA-seq, qRT-PCR was performed. Nineteen DEGs were selected for qRT-PCR, including two transcription factors (one bHLH and one WRKY), seven members (TIFYs) of the JAZ gene family, seven genes involved in starch and sucrose metabolism (BraA01g006680.3C, BAM2, SPS1, BGLUs), and three other genes (Fig. 10). The qRT-PCR results were strongly correlated with RNA-seq data for both CT (R2 = 0.8202) (Fig. S3A) and CS (R2 = 0.8230) (Fig. S3B), demonstrating the reliability of the RNA-seq expression profile.
In order to resist cold, plants change their morphological, physiological and biochemical properties, which involve molecular changes . In this study, we used the RNA-seq technique to analyze the transcriptome of two wucai genotypes after 24 h cold treatment and 24 h normal treatment to determine the response to cold stress. The DEGs common to both genotypes and unique to CT (LT vs NT) were identified and further analyzed, which could provide insights into the candidate genes and metabolic pathways underlying cold tolerance in wucai.
Based on our team’s research results, six genotypes of wucai were screened for tolerance, and two divergent varieties were chosen for further study. REC was negatively correlated with plasma membrane integrity and was related to the strength of plant stress resistance [25, 26]. MDA is a marker to measure lipid peroxidation in plant cells, which is often used to evaluate plant tolerance to biological or abiotic stresses . REC and MDA are commonly used to measure membrane damage and cell stability . Plant photosynthesis is very sensitive to cold stress, and PIabs is a multi-parametric expression that combines the three main functional steps taking place in PSII (light energy absorption, excitation energy trapping, and conversion of excitation energy to electron transport). PIabs reflect the capture of light energy by PSII reaction center and the ability of photosynthetic electron transfer between the two photosystems . PIabs could be used to select individuals for analysis or to screen wucai genotypes for enhanced cold tolerance, and it is positively correlated with the cold tolerance of plants . The change of relative variable fluorescence Vj reflects the transfer of electrons on the PSII electron receptor side from QA to QB. The change of Vj further showed that PSII receptor was affected under cold stress. Vj reflects the photosynthetic capacity which is negatively correlated with the photosynthetic capacity of plants . The changes of these indexes in wucai were consistent with the above research.
A previous study reported that ROS has a dual function in plants . Under normal conditions, the intracellular ROS level is low, which regulates plant growth and development and stress responses . Under abiotic stress, ROS accumulates and has a toxic effect on cells . The H2O2 content and O2·− generation rate of CS were higher than those of CT under cold stress (Fig. 2D, E). The T-AOC is one of the important indexes reflecting the capacity of the non-enzymatic intracellular antioxidant defense system . The tolerance of species and genotypes to abiotic stress is related to the antioxidant capacity of leaves . The increase in the T-AOC of CT was greater than that of CS under cold stress (Fig. 2F). It was consistent with research results on the resistance of Brassica napus to salt and cold stress, and T-AOC might be a potentially useful phenotypic marker of stress resistance . These results indicated that CT has a stronger ability to protect cells from toxicity and has a stronger defense against cold damage. There was no significant change of PIabs in CT, while a significant decrease was observed in CS (Fig. 2E), indicating that the photosynthetic mechanism of CS was damaged. The changes in these parameters are consistent with W18 having greater cold-tolerance than Sw-1.
Gene function and expression studies are needed to identify the key genes for cold resistance in wucai. Based on a comparison of the DEGs of the two groups, the number of DEGs in CS (LT vs NT) was higher, while the number of up-regulated DEGs in CS (LT vs NT) was fewer than that in CT (LT vs NT) (Fig. 3A). In the DEGs specific to CT (LT vs NT), 959 genes were up-regulated, including TIFY, calmodulin, LOX, and transcription factors, WRKY and MYB, which might be involved in greater cold tolerance in CT. GO enrichment specifically in CT included the response to cadmium ion (GO:0,046,686), response to jasmonic acid (GO:0,009,753), and response to wounding (GO:0,009,611). JA and its derivatives play an important role in the resistance to abiotic stress . Up-regulated DEGs in response to cadmium ion, such as APX1 and GPX6, have been reported to participate in cold response in Brassica juncea by scavenging ROS . The DEGs involved in these three GO might participate in cold tolerance in wucai.
A large number of studies have confirmed that low-temperature stress can induce the accumulation of soluble sugars to response to stress tolerance [37, 38]. Low temperature also promotes increased levels of enzymes, such as sucrose-phosphate synthase (SPS), sucrose synthase (SUS), and invertase (INV), which ultimately increased sucrose levels . The accumulation of soluble sugar in rice and cold-treated Arabidopsis thaliana  is associated with cold tolerance. The key genes involved in the conversion of fructose, sucrose, and glucose, such as SPS1, SUS6, beta-fructofuranosidase, insoluble isoenzyme 5 (CWINV5), and BGLUs, were up-regulated in CT and CS (Fig. 6A, B). In both comparison groups, starch synthesis was not inhibited (Fig. 6C, D), and a significant trend was observed that polysaccharides were degraded to disaccharides and soluble, simple sugars. Starch is hydrolyzed/decomposed into maltose/dextrin by alpha amylase 3 (AMY3) and BAM2 and is then hydrolyzed into glucose by 4-alpha-glucanotransferase (DPE2). It was reported that sucrose and starch metabolism was significantly enriched and the content of soluble sugar was significantly increased after cold stress in Poa pratensis with different cold tolerance . Poa pratensis is a cold tolerant species, and this result was consistent with ours. It might indicate that starch and sucrose metabolism was involved in the response of wucai with different cold tolerance to cold stress.
GSH is involved in ROS clearance in the non-enzymatic mode, but the possibility of direct ROS clearance by GSH is considered low . GSH participates in the ROS scavenging enzymatic reaction by regulating the activities of related enzymes in the antioxidant enzyme system (APX, mono dehydroascorbate reductase, dehydroascorbate reductase, and glutathione reductase) and by indirectly scavenging ROS. In this study, APX1 was significantly up-regulated in the glutathione metabolism pathway in both CT (LT vs NT) and CS (LT vs NT) (Fig. 7B). GSH is mainly involved in regulating the response of GPX (phospholipid hydroperoxide glutathione peroxidase) and GST. GPX is a member of the proline oxidase antioxidant enzyme family, which uses GSH to remove H2O2 and reduce the accumulation of lipids and organic hydroperoxides . GST can also be coupled with GPX to activate GPX and participate in H2O2 scavenging. One transcription (BraA09g026940.3C) of GPX6 was up-regulated in both groups, and the other transcription (BraA08g032660.3C) was down-regulated in CT (LT vs NT) but was not detected in CS (LT vs NT). Members of the GST family, specifically GSTU11, GSTU24, and GSTF12, were significantly up-regulated in both groups. These results suggested that the glutathione metabolism pathway, as an important component of the antioxidant process in plants, might play a role in both genotypes of wucai in response to cold stress. This may be a reason why wucai is a semi-hardy species.
LOX initiates the first step of α-linolenic acid oxidation, JA, and MeJA synthesis . It has been reported that there are 13 LOX genes in Brassica rapa, and LOX gene participate in the response to stress . In the present study, LOX3 was differentially expressed and up-regulated in CT. LOX4 was differentially expressed in CS and down-regulated (Fig. 8C). MeJA alleviates low temperature stress in tomato  and cowpea (Vigna sinensis)  by increasing antioxidant synthesis and the activity of some defense compounds. JMT is a key enzyme involved in the conversion of JA to MeJA . It was down-regulated in CS and was not detected in CT. The MeJA content in CT was significantly increased compared with CS, indicating that MeJA might have a stronger response to cold stress than CS (Fig. 8D). JA-mediated signal transduction pathways play an important role in various stress responses . Previous studies have shown that JA acts as a key upstream signal to regulate the ICE-CBF/DREB1 transcription pathway to prevent frost damage in Arabidopsis thaliana . TIFY is a plant-specific transcription factor family that includes four subfamilies, TIFY, JAZ, PPD, and ZML. Studies have shown that members of the JAZ subgroup are involved in the response to abiotic stress by participating in the JA signal transduction pathway , and MeJA induces a defense response by activating JA-related genes, thus endowing pepper with tolerance to chilling injury . And JAZs encode the plant-specific proteins that are involved in JA signaling and stress tolerance . In this study, JAZ genes identified in CT (LT vs NT) were up-regulated and were down-regulated in CS (LT vs NT) (Fig. 8C). MYC2 regulates downstream transcription factors in JA signal transduction, and its downregulation in CT is lower than that in CS. Therefore, we hypothesize that α-linolenic acid metabolism and the signal transduction pathway of JA might be the main reasons for the difference in cold tolerance between CT and CS. Thus, the molecular mechanisms of these JA-related genes need to be further investigated.
TFs are important upstream regulators that are key to the regulation of gene expression in plants under abiotic stress [54, 55]. Previous studies have reported that overexpression of bHLH92 can improve the cold tolerance of Arabidopsis thaliana . The induction of AtERF53  and expression of WRKY53 and WRKY54  increased biotic and abiotic stress tolerance in transgenic Arabidopsis. These TFs were identified only in CT, indicating that that they may play a role in the cold tolerance of wucai (Fig. 9B). Thus, future research could focus on the role of these transcription factors in metabolic pathways to understand the regulatory mechanism of genes in wucai.
In this study, through phenotypic observation and physiological index measurements, we identified the cold-tolerant (W18) and cold-sensitive (Sw-1) genotypes among the six wucai genotypes. Transcriptome analysis was used to screen DEGs related to different cold tolerance levels of wucai, which has a certain reference value for cold-tolerance breeding of wucai. DEGs in the two groups might response to cold stress by participating in starch and sucrose metabolism and glutathione metabolism. The cold tolerance of CT was mainly related to α-linolenic acid metabolism and the JA signal transduction pathway and involved transcription factors, such as bHLH92, ERF53, WRKY53, and WRKY54. These findings provide genetic resources and a theoretical basis for subsequent studies on the improvement of gene function and cold tolerance of wucai and also have a certain reference value for the analysis of the cold tolerance of wucai and other similar species of winter rapeseed (Brassica juncea and Brassica napus).
Plant material and treatment
We selected six wucai genotypes, ws-1, ws-2, Sw-1, Sw-3, W18, and Ta2, which were provided by the Vegetable Genetics and Breeding Laboratory of Anhui Agricultural University. Seeds were sown in plug trays and then transplanted into pots containing substrate and vermiculite with a volume ratio of 2:1. Seedlings were grown in a growth chamber at 25 °C (day) and 15 °C (night) with 300 μmol·m−2·s−1 photon flux density and 70% relative humidity under a 16 /8 h (light/night) photoperiod. For cold treatment, seedlings at the 7-leaf stage were transferred to the following conditions: 9 °C (day) and 4 °C (night) with 300 μmol·m−2·s−1 photon flux density and 70% relative humidity under a 16 /8 h light/dark photoperiod, whereas the control plants were cultured under the normal environmental conditions mentioned above. After 3 d cold treatment, biochemical parameters, such as REC, the MDA content, PIabs, and Vj, were measured and compared with the control.
Measurement of MDA, REC, O2·−, and H2O2
The MDA content was assessed following the methods previously described by Mohammadi et al. with minor modifications . Fresh leaves (0.2 g) were homogenized in 2 mL 10% trichloroacetic acid (TCA) and centrifuged at 4000 rpm for 10 min. Then, 2-mL supernatant was mixed with 2 mL thiobarbituric acid (TBA) (0.6%), heated at 100 °C for 15 min, and cooled at room temperature. The absorbance of each aliquot was measured at 450, 532, and 600 nm. The MDA content was calculated using the equation:
MDA (μmol·g−1 FW) = 6.45 × (A532 − A600) − 0.56 × A450.
Relative electrolyte leakage was determined according to the method reported by Bajji et al. with minor modifications . The leaves were perforated with a hole punch with a radius of 9.5 mm. Three discs were placed in a 20 mL tube containing 10 mL of ultrapure water and then placed in a thermostated water bath at 25 °C for 30 min. After this, the conductivity (EC1) of water was measured using a Thermo OrionSTARA HB conductivity meter (Thermo Orion., Waltham, MA, USA). The tube was heated in boiling water for 30 min, then cooled to room temperature, and the conductivity (EC2) was measured again. The final REC was equal to the percentage of EC1/EC2.
The O2·− generation rate and the H2O2 content were measured using Solarbio reagent kits (Cat#BC3595 and Cat#BC1290, Beijing Solarbio Science & Technology Co. Ltd, Beijing, China).
Analysis of Chlorophyll Fluorescence Parameters
PIabs = (RC/ABS) • [φPo/(1–φPo)] [ψo/(1–ψo)] and is a performance index with an absorption basis. Vj = (Fj–Fo)/(Fm–Fo), which is the relative variable fluorescence intensity at the J-step. To measure PIabs and Vj, the leaves were dark-adapted for 30 min by special clips before being illuminated with saturating light (1 s), after which they were measured using a continuous excitation fluorometer Pocket Plant Efficiency Analyzer (Pocket PEA, Hansatech, UK).
RNA extraction, cDNA library construction, and sequencing
A total of four group samples [two materials (leaves of W18 and Sw-1) under two treatments (control and 10 °C/4 °C 24 h)] with three biological replications were used for RNA-seq. The total RNA of samples mentioned above was extracted using a mir-Vana miRNA isolation kit (Ambion, TX, USA) according to the manufacturer’s instructions. The total RNA was quantified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Samples with an RNA Integrity Number ≥ 7.0 and 28S/18S ratio ≥ 0.7 were subjected to subsequent analysis. The libraries were constructed using the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) following the manufacturer’s instructions. Then, these libraries were sequenced on the Illumina HiSeqTM 2500 platform (BiomarkerBiotech, Beijing, China) to generate 125 bp/150 bp paired-end reads.
Sequence assembly, annotation, and identification of the DEGs
Clean reads were obtained by removing adaptor sequences, more than 10% N bases, and low quality (Q ≤ 20) reads with more than 50% bases. The reads were mapped to the B. rapa reference genome  by hisat2 . Read counts per gene were expressed as FPKM, and unigene abundance differences between the samples were calculated based on the ratio of the FPKM values and false discovery rate (FDR). DESeq software  was used to standardize the counts of each sample gene (use basemean value to estimate the expression), calculate the difference multiple, and use NB (negative binomial distribution test) to test the difference significance of the reads number. Finally, screen the differential protein coding genes according to the difference multiple and difference significance test results. Genes with FDR ≤ 0.05 and |log2 (Fold Change) |≥ 2 were considered DEGs.
Hierarchical clustering of the DEGs was analyzed by HemI software . Statistical enrichment of the DEGs was carried out in KEGG pathways using KOBAS software , and log2 (Fold Change) values were used in the gene expression heatmap.
Measurement of carbohydrate, glutathione metabolism products, and MeJA content
The contents of sucrose, fructose, and starch were determined by anthrone colorimetry  with minor modifications. Glucose, GSH, and GSSG contents were measured using a biochemical reagent kit (Cat#BC2500, Cat#BC1175, and Cat#BC1180, respectively; Beijing Solarbio Science & Technology Co. Ltd, Beijing, China). MeJA was determined by enzyme-linked immunosorbent assay.
Identification of transcription factors
TFs were identified by analyzing InterProScan domain patterns in sequences with high coverage, and the sensitivity was analyzed using PlantTFcat tools .
A total of 15 transcripts were selected to verify the RNA-seq analysis. The transcript-specific primers used in this study were designed using Primer6 and were then synthesized by General Biosystems (Chuzhou, China). The primers used for qRT-PCR are listed in Table S6. The qRT-PCR reactions were performed with the AceQ qPCR SYBR GREEN Master Mix (Vazyme Biotechnology Co., Nanjing, China). BnaActin  was selected as the internal standard to calculate the relative expression levels as follows: FC = 2−△△CT . Three biological repeats for each sample for each gene were performed.
Availability of data and materials
The raw RNA-Seq datasets are available in the Sequence Read Archive of National Center for Biotechnology Information (https://dataview.ncbi.nlm.nih.gov/object/PRJNA735896?reviewer=fnp3ih4ojrq5ohq5vmstrr5koe; accession number: PRJNA735896).
Allene oxide cyclase
C-repeat binding factor
Cold sensitivity in low temprature
Cold tolerance in low temprature
Beta-fructofuranosidase, insoluble isoenzyme
Cytochrome P450 superfamily protein
Differentially expressed genes
Dehydration-responsive element-binding proteins 1
Ethylene-responsive transcription factor
False discovery rate
The expected number of fragments per kilobase of transcripts per million mapped fragments
Phospholipid hydroperoxide glutathione peroxidase
- H2O2 :
JA carboxyl methyltransferase
Kyoto encyclopedia of genes and genomes
Serine/threonine protein kinase 2
- O2·− :
- O2 1 :
Princinple component analysis
Pearson correlation coefficient
- PIabs :
The performance index on an absorption basis
Quantitative real-time PCR
Relative electrolyte conductivity
Reactive oxygen species
Reads per kilobase per million mapped reads
Total antioxidant capacity
Pearce RS. Plant freezing and damage. Ann Bot. 2001;87(4):417–24.
Guo X, Xu S, Chong K. Cold Signal Shuttles from Membrane to Nucleus. Mol Cell. 2017;66(1):7–8.
Mittler R, Vanderauwera S, Suzuki N, Miller G, Tognetti VB, Vandepoele K, Gollery M, Shulaev V, Van Breusegem F. ROS signaling: the new wave? Trends Plant Sci. 2011;16(6):300–9.
Krasensky J, Jonak C. Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory networks. J Exp Bot. 2012;63(4):1593–608.
Qi J, Song C-P, Wang B, Zhou J, Kangasjarvi J, Zhu J-K, Gong Z. Reactive oxygen species signaling and stomatal movement in plant responses to drought stress and pathogen attack. J Integr Plant Biol. 2018;60(9):805–26.
Mittler R. ROS Are Good. Trends Plant Sci. 2017;22(1):11–9.
Farooq MA, Niazi AK, Akhtar J, Saifullah, M Farooq, Z Souri, N Karimi, Z Rengel. Acquiring control: The evolution of ROS-Induced oxidative stress and redox signaling pathways in plant stress responses. Plant Physiology and Biochemistry. 2019;141:353–369
Wasternack C, Feussner I. The Oxylipin Pathways: Biochemistry and Function. In: Annual Review of Plant Biology, Vol 69. Edited by Merchant SS, vol. 69; 2018: 363–386.
Li J, Zhang K, Meng Y, Hu J, Ding M, Bian J, Yan M, Han J, Zhou M. Jasmonic acid/ethylene signaling coordinates hydroxycinnamic acid amides biosynthesis through ORA59 transcription factor. Plant J. 2018;95(3):444–57.
Pauwels L, Goossens A. The JAZ Proteins: A Crucial Interface in the Jasmonate Signaling Cascade. Plant Cell. 2011;23(9):3089–100.
Hsieh E-J, Cheng M-C, Lin T-P. Functional characterization of an abiotic stress-inducible transcription factor AtERF53 in Arabidopsis thaliana. Plant Mol Biol. 2013;82(3):223–37.
Tang K, Zhao L, Ren Y, Yang S, Zhu J-K, Zhao C. The transcription factor ICE1 functions in cold stress response by binding to the promoters of CBF and COR genes. J Integr Plant Biol. 2020;62(3):258–63.
Wu R, Wang Y, Wu T, Xu X, Han Z. MdMYB4, an R2R3-Type MYB Transcription Factor, Plays a Crucial Role in Cold and Salt Stress in Apple Calli. Journal of the American Society for Horticultural Science. 2017;142(3):209.
Niu C-F, Wei W, Zhou Q-Y, Tian A-G, Hao Y-J, Zhang W-K, Ma B, Lin Q, Zhang Z-B, Zhang J-S, et al. Wheat WRKY genes TaWRKY2 and TaWRKY19 regulate abiotic stress tolerance in transgenic Arabidopsis plants. Plant, Cell Environ. 2012;35(6):1156–70.
Ma Y, Shukla V, Merewitz EB. Transcriptome analysis of creeping bentgrass exposed to drought stress and polyamine treatment. Plos One. 2017;12(4):e0175848.
Ma L, Coulter JA, Liu L, Zhao Y, Chang Y, Pu Y, Zeng X, Xu Y, Wu J, Fang Y, et al. Transcriptome Analysis Reveals Key Cold-Stress-Responsive Genes in Winter Rapeseed (Brassica rapa L.). International Journal of Molecular Sciences. 2019;20(5):1071.
Zhang Y, Zhang Y, Lin Y, Luo Y, Wang X, Chen Q, Sun B, Wang Y, Li M, Tang H. A Transcriptomic Analysis Reveals Diverse Regulatory Networks That Respond to Cold Stress in Strawberry (Fragaria x ananassa). Int J Gen. 2019;1-13.
Shen Q, Zhang S, Liu S, Chen J, Ma H, Cui Z, Zhang X, Ge C, Liu R, Li Y, et al. Comparative Transcriptome Analysis Provides Insights into the Seed Germination in Cotton in Response to Chilling Stress. International Journal of Molecular Sciences. 2020;21(6):2067.
Strickler SR, Bombarely A, Mueller LA. DESIGNING A TRANSCRIPTOME NEXT-GENERATION SEQUENCING PROJECT FOR A NONMODEL PLANT SPECIES. Am J Bot. 2012;99(2):257–66.
Zhao M, Yuan L, Wang J, Xie S, Zheng Y, Nie L, Zhu S, Hou J, Chen G, Wang C. Transcriptome analysis reveals a positive effect of brassinosteroids on the photosynthetic capacity of wucai under low temperature. Bmc Genomics. 2019;20(1):1–19.
Yuan L, Wang J, Xie S, Zhao M, Nie L, Zheng Y, Zhu S, Hou J, Chen G, Wang C. Comparative Proteomics Indicates That Redox Homeostasis Is Involved in High- and Low-Temperature Stress Tolerance in a Novel Wucai (Brassica campestris L.) Genotype. International Journal of Molecular Sciences. 2019;20(15):3760.
Min X, Liu Z, Wang Y, Liu W. Comparative transcriptomic analysis provides insights into the coordinated mechanisms of leaves and roots response to cold stress in Common Vetch. Industrial Crops and Products. 2020;158:112949.
Dai Y, Yuan L, Zhang S, Wang J, Xie S, Zhao M, Chen G, Sun R, Wang C. Comprehensive Evaluation for Cold Tolerance in Wucai (Brassica campestris L.) by the Performance Index on an Absorption Basis (PIabs). Agronomy-Basel. 2019;9(2):61.
Shinozaki K, Yamaguchi-Shinozaki K, Seki M. Regulatory network of gene expression in the drought and cold stress responses. Curr Opin Plant Biol. 2003;6(5):410–7.
Chakraborty N, Basak J. Exogenous application of methyl jasmonate induces defense response and develops tolerance against mungbean yellow mosaic India virus in Vigna mungo. Funct Plant Biol. 2019;46(1):69–81.
Armstrong JJ, Takebayashi N, Sformo T, Wolf DE. COLD TOLERANCE IN ARABIDOPSIS KAMCHATICA. Am J Bot. 2015;102(3):439–48.
Shi H, Li R, Cai W, Liu W, Wang C, Lu Y. Increasing Nitric Oxide Content in Arabidopsis thaliana by Expressing Rat Neuronal Nitric Oxide Synthase Resulted in Enhanced Stress Tolerance. Plant Cell Physiol. 2012;53(2):344–57.
Belen Montero-Palmero M, Martin-Barranco A, Escobar C, Hernandez LE. Early transcriptional responses to mercury: a role for ethylene in mercury-induced stress. New Phytol. 2014;201(1):116–30.
Strauss AJ, Kruger GHJ, Strasser RJ, Van Heerden PDR. Ranking of dark chilling tolerance in soybean genotypes probed by the chlorophyll a fluorescence transient O-J-I-P. Environ Exp Bot. 2006;56(2):147–57.
Xia X, Wang Y, Zhou Y, Tao Y, Mao W, Shi K, Asami T, Chen Z, Yu J. Reactive Oxygen Species Are Involved in Brassinosteroid- Induced Stress Tolerance in Cucumber. Plant Physiol. 2009;150(2):801–14.
Gill SS, Tuteja N. Reactive oxygen species and antioxidant machinery in abiotic stress tolerance in crop plants. Plant Physiol Biochem. 2010;48(12):909–30.
Parameshwaran K, Irwin MH, Steliou K, Pinkert CA. D-Galactose Effectiveness in Modeling Aging and Therapeutic Antioxidant Treatment in Mice. Rejuvenation Res. 2010;13(6):729–35.
Ashraf M, Foolad MR. Roles of glycine betaine and proline in improving plant abiotic stress resistance. Environ Exp Bot. 2007;59(2):206–16.
Xu Z, Peng L, Xiao P, Li H, Feng X, Hong X. Exogenous application of poly--glutamic acid enhances stress defense in Brassica napus L. seedlings by inducing cross-talks between Ca2+, H2O2, brassinolide, and jasmonic acid in leaves. Plant Physiology & Biochemistry. 2017;118:460–70.
Ali MS, Baek KH. Jasmonic Acid Signaling Pathway in Response to Abiotic Stresses in Plants. International Journal of Molecular Sciences. 2020;21(2):621.
Verma D, Upadhyay SK, Singh K. Characterization of APX and APX-R gene family in Brassica juncea and B. rapa for tolerance against abiotic stresses. Plant Cell Reports. 2021;1-22.
Rolland F, Moore B, Sheen J. Sugar sensing and signaling in plants. Plant Cell. 2002;14:S185–205.
Azymi S, Sofalian O, Jahanbakhsh GS, Khomari S. Effect of chilling stress on Soluble Protein, sugar and Prolin accumulation in cotton (Gossypium hirsutum L.) genotypes. 2012;4(12):825-30.
Ramon M, Rolland F, Sheen J. Sugar sensing and signaling. The arabidopsis book. 2008;6:e0117–e0117.
Cook D, Fowler S, Fiehn O, Thomashow MF. A prominent role for the CBF cold response pathway in configuring the low-temperature metabolome of Arabidopsis. Proc Natl Acad Sci USA. 2004;101(42):15243–8.
Dong W, Ma X, Jiang H, Zhao C, Ma H. Physiological and transcriptome analysis of Poa pratensis var. anceps cv. Qinghai in response to cold stress. Bmc Plant Biology. 2020;20(1):1–18.
Rao A, Reddy AR. Glutathione Reductase: A Putative Redox Regulatory System in Plant Cells: Sulfur Assimilation and Abiotic Stress in Plants. Berlin, Heidelberg: Springer; 2008.
Bela K, Horvath E, Galle A, Szabados L, Tari I, Csiszar J. Plant glutathione peroxidases: Emerging role of the antioxidant enzymes in plant development and stress responses. J Plant Physiol. 2015;176:192–201.
Mao L, Pang H, Wang G, Zhu C. Phospholipase D and lipoxygenase activity of cucumber fruit in response to chilling stress. Postharvest Biol Technol. 2007;44(1):42–7.
Rai AN, Mandliya T, Kulkarni P, Rao M, Suprasanna P. Evolution and Transcriptional Modulation ofLipoxygenaseGenes Under Heat, Drought, and Combined Stress inBrassica rapa. Plant Mol Biol Report. 2021;39(1):60–71.
Zhang X, Sheng J, Li F, Meng D, Shen L. Methyl jasmonate alters arginine catabolism and improves postharvest chilling tolerance in cherry tomato fruit. Postharvest Biol Technol. 2012;64(1):160–7.
Fan L, Wang Q, Lv J, Gao L, Zuo J, Shi J. Amelioration of postharvest chilling injury in cowpea (Vigna sinensis) by methyl jasmonate (MeJA) treatments. Sci Hortic. 2016;203:95–101.
Shi J, Wang J, Lv H, Peng Q, Schreiner M, Baldermann S, Lin Z. Integrated proteomic and metabolomic analyses reveal the importance of aroma precursor accumulation and storage in methyl jasmonate-primed tea leaves. Horticulture Research. 2021;8(1):1–14.
Piotrowska A, Bajguz A. Conjugates of abscisic acid, brassinosteroids, ethylene, gibberellins, and jasmonates. Phytochemistry. 2011;72(17):2097–112.
Hu Y, Jiang L, Wang F, Yu D. Jasmonate Regulates the INDUCER OF CBF EXPRESSION-C-REPEAT BINDING FACTOR/DRE BINDING FACTOR1 Cascade and Freezing Tolerance in Arabidopsis. Plant Cell. 2013;25(8):2907–24.
Ye H, Du H, Tang N, Li X, Xiong L. Identification and expression profiling analysis of TIFY family genes involved in stress and phytohormone responses in rice. Plant Mol Biol. 2009;71(3):291–305.
Shin S, Park M, Choi J, Kim J. Gene network underlying the response of harvested pepper to chilling stress. J Plant Physiol. 2017;219:112–22.
Zhu D, Cai H, Luo X, Bai X, Deyholos MK, Chen Q, Chen C, Ji W, Zhu Y. Over-expression of a novel JAZ family gene from Glycine soja, increases salt and alkali stress tolerance. Biochem Biophys Res Commun. 2012;426(2):273–9.
Huang Q, Wang Y, Li B, Chang J, Chen M, Li K, Yang G, He G. TaNAC29, a NAC transcription factor from wheat, enhances salt and drought tolerance in transgenic Arabidopsis. Bmc Plant Biology. 2015;15:1–15.
Nakashima K, Takasaki H, Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. NAC transcription factors in plant abiotic stress responses. Biochimica Et Biophysica Acta-Gene Regulatory Mechanisms. 2012;1819(2):97–103.
Jiang Y, Yang B, Deyholos MK. Functional characterization of the Arabidopsis bHLH92 transcription factor in abiotic stress. Mol Genet Genomics. 2009;282(5):503–16.
Kuki Y, Ohno R, Yoshida K, Takumi S. Heterologous expression of wheat WRKY transcription factor genes transcriptionally activated in hybrid necrosis strains alters abiotic and biotic stress tolerance in transgenic Arabidopsis. Plant Physiol Biochem. 2020;150:71–9.
Mohammadi M, Modarres-Sanavy SAM, Pirdashti H, Zand B, Tahmasebi-Sarvestani Z. Arbuscular mycorrhizae alleviate water deficit stress and improve antioxidant response, more than nitrogen fixing bacteria or chemical fertilizer in the evening primrose. Rhizosphere. 2019;9:76–89.
Bajji M, Kinet JM, Lutts S. The use of the electrolyte leakage method for assessing cell membrane stability as a water stress tolerance test in durum wheat. Plant Growth Regul. 2002;36(1):61–70.
Cheng F, Liu S, Wu J, Fang L, Sun S, Liu B, Li P, Hua W, Wang X. BRAD, the genetics and genomics database for Brassica plants. Bmc Plant Biology. 2011;11:1–6.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357-U121.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15(12):1–21.
Deng W, Wang Y, Liu Z, Cheng H, Xue Y. HemI: A Toolkit for Illustrating Heatmaps. Plos One. 2014;9(11):e111988.
Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li C-Y, Wei L. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Research. 2011;39:W316–22.
Laurentin A, Edwards CA. A microtiter modification of the anthrone-sulfuric acid colorimetric assay for glucose-based carbohydrates. Anal Biochem. 2003;315(1):143–5.
Dai X, Sinharoy S, Udvardi M, Zhao PX. PlantTFcat: an online plant transcription factor and transcriptional regulator categorization and analysis tool. Bmc Bioinformatics. 2013;14:1–6.
Yang X, Zhao T, Rao P, Gao K, Yang X, Chen Z, An X. Transcriptome profiling of Populus tomentosa under cold stress. Ind Crops Prod. 2019;135:283–93.
We thank all contributors for their work and would like to thank the reviewers for their valuable comments and suggestions. We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.
This work was supported by Chinese Cabbage Germplasm Bank in Anhui Province (202001n06030004), Germplasm Resource Nursery of Wucai (201901n06030006), Natural science research project of Anhui colleges and Universities (KJ2020ZD11) and the Graduate Innovation Found of Anhui Agricultural University (2021yjs-7).
Ethics approval and consent to participate
Lines ws-1, ws-2, Sw-1, Sw-3, W18 and Ta2 were screened by Chenggang Wang. All seed materials were kept in the Vegetable Breeding Laboratory of the College of Horticulture, Anhui Agricultural University. No specific permits were required for the described field studies. The location is not privately-owned or protected in any way, and the field studies did not involve endangered or protected species. We complied with the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Fig. S1. Regional distribution, density distribution and principal component analysis of gene expression. Fig. S2. Hierarchical clustering DEGs in CS (LTvsNT). Each column represents a comparison group, and each row represents a gene. Fig. S3. qRT-PCR validation of expression profiles obtained by RNA-Seq in CT and CS under cold stress. Table S1. Summary of sequence assembly after illumine sequencing. Table S2. Number of reads sequenced and mapped to the Brassica rapa genome. Table S3. List of the core genes set involved in CS response to cold stress. Table S4. Top30 up GO enrichment analysis in CT(LTvsNT). Table S5. Top30 up GO enrichment analysis in CS(LTvsNT). Table S6. Primer pairs used to detect the expression of selected genes.
About this article
Cite this article
Wang, C., Zhang, M., Zhou, J. et al. Transcriptome analysis and differential gene expression profiling of wucai (Brassica campestris L.) in response to cold stress. BMC Genomics 23, 137 (2022). https://doi.org/10.1186/s12864-022-08311-3
- Brassica campestris
- Cold stress
- α-linolenic acid metabolism