Transcriptome Analysis in Chinese Cabbage (Brassica rapa ssp. pekinensis) Provides the Role of Glucosinolate Metabolism in Response to Drought Stress

Although drought stress is one of the most limiting factors in growth and production of Chinese cabbage (Brassica rapa L. ssp. pekinensis), the underlying biochemical and molecular causes are poorly understood. In the present study, to address the mechanisms underlying the drought responses, we analyzed the transcriptome profile of Chinese cabbage grown under drought conditions. Drought stress transcriptionally activated several transcription factor genes, including AP2/ERFs, bHLHs, NACs and bZIPs, and was found to possibly result in transcriptional variation in genes involved in organic substance metabolic processes. In addition, comparative expression analysis of selected BrbZIPs under different stress conditions suggested that drought-induced BrbZIPs are important for improving drought tolerance. Further, drought stress in Chinese cabbage caused differential acclimation responses in glucosinolate metabolism in leaves and roots. Analysis of stomatal aperture indicated that drought-induced accumulation of glucosinolates in leaves directly or indirectly controlled stomatal closure to prevent water loss, suggesting that organ-specific responses are essential for plant survival under drought stress condition. Taken together, our results provide information important for further studies on molecular mechanisms of drought tolerance in Chinese cabbage.


Introduction
Because the increasing world population and worldwide climate change affects agriculture in several ways, an intergovernmental panel on climate change has concluded that increased concentrations of greenhouse gases will lead to dry conditions in the subtropics, creating widespread drought stress in agricultural regions [1]. Drought is considered as the dominant factor reducing crop growth and yield, and is expected to cause serious problems in plant growth and crop production over more than 50% of all agricultural lands by 2050 [2]. In addition, it has been shown that drought sensitivity in crops has increased over the past twenty years [3]. Therefore, improving drought tolerance is a major goal for crop breeders. Under drought conditions, plants display various physiological and biochemical responses at the cellular and whole-plant levels, leading to a range of specific and nonspecific phenotypes [4]. These responses include stomatal movement, repression of These physiological changes suggested that the drought stress treatment was effective, and the leaves and roots of Chinese cabbage plants grown under control conditions (Stage 1) and drought stress conditions (Stage 2 and Stage 3) were harvested to unravel the molecular basis of the response.  The data are presented as mean ± standard error (SE). * p < 0.05, ** p < 0.01, and *** p < 0.001 represent the significant differences in comparison with Stage 1 (S1).

RNA-Seq and Identification of Drought-Responsive Genes
During the last decade, several exciting studies have reported the development of new diseaseresistant crops using genetic engineering tools, but efforts to generate drought-tolerant crops has been less successful [18]. This might be due to the complex responses of crops to drought stress. In fact, drought stress induces various morphological, physiological, and biochemical changes in crops controlled by numerous small-effect loci and hundreds of genes [19], indicating that knowledge of the biochemical and molecular responses to drought should contribute to improving plant tolerance under water-limited conditions. Therefore, to obtain a global overview of molecular mechanisms involved in drought stress response, total RNA was extracted from roots and leaves of control and drought-treated plants and used for generation of transcriptome libraries on an Illumina HiSeq™ 2500 sequencing platform. After removing low-quality reads, 46 to 54 million clear reads (4.64 to 5.49 Gb) from each sample were acquired for further analysis (Table 1). All reads were deposited in the National Agricultural Biotechnology Information Center (NABIC, http://nabic.rda.go.kr). To identify the DEGs in the control and drought-treated groups, we constructed six comparison groups, and DEGs were defined as genes with value of |log2 (fold change)| ≥ 1. As shown in Figure 2A, a total of 7535 DEGs (3557 up-regulated and 3978 downregulated) were detected comparing the L1 and L2 libraries. Comparative analysis of the L2 and L3 libraries showed that the expression of 2666 genes was increased and that of 2608 genes was downregulated in L3. Upon drought stress exposure, more DEGs were down-regulated in roots than in leaves (Figure 2A), indicating that the response mechanism against drought stress differs between roots and leaves. A total of 1280 genes were exclusively up-regulated in L2, and the expression levels of 1282 genes were increased in L3 but not in L2. In roots, transcripts of 1135 genes were downregulated solely in R2, whereas 1752 transcripts were down-regulated exclusively in R3 but not in R2 Figure 1. Physiological response to drought stress in Chinese cabbage. (A) Phenotypes of Chinese cabbage plants after exposure to drought stress for four days. S1 (Stage 1), normal water supply; S2 (Stage 2), water withholding for 2 d (soil water content at 20%); S3 (Stage 3), water withholding for 4 d (soil water content at 5%). Relative water content (B) and changes in levels of H 2 O 2 (C), MDA (D), and protein carbonylation (E) after drought treatment were determined. Values are averages from three biological replicates consisting three plants per sample were used. The data are presented as mean ± standard error (SE). * p < 0.05, ** p < 0.01, and *** p < 0.001 represent the significant differences in comparison with Stage 1 (S1).

RNA-Seq and Identification of Drought-Responsive Genes
During the last decade, several exciting studies have reported the development of new disease-resistant crops using genetic engineering tools, but efforts to generate drought-tolerant crops has been less successful [18]. This might be due to the complex responses of crops to drought stress. In fact, drought stress induces various morphological, physiological, and biochemical changes in crops controlled by numerous small-effect loci and hundreds of genes [19], indicating that knowledge of the biochemical and molecular responses to drought should contribute to improving plant tolerance under water-limited conditions. Therefore, to obtain a global overview of molecular mechanisms involved in drought stress response, total RNA was extracted from roots and leaves of control and drought-treated plants and used for generation of transcriptome libraries on an Illumina HiSeq™ 2500 sequencing platform. After removing low-quality reads, 46 to 54 million clear reads (4.64 to 5.49 Gb) from each sample were acquired for further analysis (Table 1). All reads were deposited in the National Agricultural Biotechnology Information Center (NABIC, http://nabic.rda.go.kr). To identify the DEGs in the control and drought-treated groups, we constructed six comparison groups, and DEGs were defined as genes with value of |log2 (fold change)| ≥ 1. As shown in Figure 2A, a total of 7535 DEGs (3557 up-regulated and 3978 down-regulated) were detected comparing the L1 and L2 libraries. Comparative analysis of the L2 and L3 libraries showed that the expression of 2666 genes was increased and that of 2608 genes was down-regulated in L3. Upon drought stress exposure, more DEGs were down-regulated in roots than in leaves (Figure 2A), indicating that the response mechanism against drought stress differs between roots and leaves. A total of 1280 genes were exclusively up-regulated in L2, and the expression levels of 1282 genes were increased in L3 but not in L2. In roots, transcripts of 1135 genes were down-regulated solely in R2, whereas 1752 transcripts were down-regulated exclusively in R3 but not in R2 ( Figure 2B). Further, 467 DEGs from leaves and 240 DEGs from roots were found to be up-regulated in both tissues, whereas 424 DEGs from leaves and 1181 DEGs from roots were found to be down-regulated in both tissues ( Figure 2B). Non-specific lipid-transfer protein 3 (XP_009132002.1), dehydrin Rab18-like (XP_013609020.1), and sucrose synthase 3 (XP_009111303.1) were among those identified to be up-regulated in both roots and leaves, whereas aquaporin TIP2-1 (XP_013653350.1) was down-regulated in both tissues.
Molecules 2018, 23, x FOR PEER REVIEW 4 of 15 ( Figure 2B). Further, 467 DEGs from leaves and 240 DEGs from roots were found to be up-regulated in both tissues, whereas 424 DEGs from leaves and 1181 DEGs from roots were found to be downregulated in both tissues ( Figure 2B). Non-specific lipid-transfer protein 3 (XP_009132002.1), dehydrin Rab18-like (XP_013609020.1), and sucrose synthase 3 (XP_009111303.1) were among those identified to be up-regulated in both roots and leaves, whereas aquaporin TIP2-1 (XP_013653350.1) was downregulated in both tissues. To understand the universal response of Chinese cabbage to drought stress, all screened DEGs in the six comparison groups were subjected to GO enrichment analysis, and 60 functional groups To understand the universal response of Chinese cabbage to drought stress, all screened DEGs in the six comparison groups were subjected to GO enrichment analysis, and 60 functional groups classified into three categories comprising "cellular component," "molecular function," and "biological process" were identified (Table S1). The category "biological process," consisting of 27 functional groups, exhibited the highest number of annotations alongside "cellular process", followed by "metabolic process", and "single-organism process". In addition, the drought-induced DEGs revealed various affected metabolic processes, mainly the organic substance metabolic processes, cellular metabolic processes and primary metabolic processes ( Figure 2C). Photosynthesis is one of the important processes affected by drought, due to changes including decrease in turgor pressure, reduction in CO 2 assimilation and diffusion, and impairment in the photosynthetic apparatus [20][21][22], suggesting that alterations in primary metabolic processes under drought conditions ( Figure 2C) are mediated via drought-induced photosynthesis inhibition.
Further, transcript abundances of genes involved in the abscisic acid (ABA) biosynthesis and degradation pathway were dynamically changed in response to drought stress, whereas a number of genes involved in biosynthesis of auxin, gibberellin, or ethylene were down-regulated ( Figure S1), suggesting that the hormonal regulation of drought stress response in Chinese cabbage is primarily mediated via ABA.

Transcription Factors and MapMan Analysis of Chinese Cabbage Genes Associated with Drought-Stress Response
Transcription factors (TFs) are known as major players in various transcriptional regulatory mechanisms in stress-induced signal transduction pathways. Therefore, understanding the role of stress-inducible TFs is important for engineering stress-tolerant plants by modulating a large set of genes [23]. In Chinese cabbage transcriptome libraries generated in the above mentioned RNA-seq, a total of 3663 TFs mainly distributed into 45 TF families was found ( Figure S2). Among these TFs, the most abundant differentially expressed TF family was the AP2/ERF family ( Figure S3). Based on the difference in copy numbers and similarity of AP2/ERF domains, this family is further divided into five subfamilies, the Apetala2 (AP2), the ethylene responsive element binding factor (ERF), the Related to ABI3/VP1 (RAV), the dehydration-responsive element-binding (DREB) protein, and the Soloist [24], However, our transcriptome libraries did not contain members of the Soloist subfamily ( Figure S4), probably due to its tissue-specific expression. In fact, expression of the Soloist gene was detected only in the bud of Chinese cabbage plants [25]. The ERF subfamily, which plays significant roles in regulating gene expression in response to biotic and abiotic stresses [6], was the most abundant class of the AP2/ERF family differentially expressed in response to drought ( Figure S4).
The bZIP TF family, which contain a basic region for DNA-binding and a leucine zipper motif for dimerization, is one of the largest and most diverse TF families in plants [26], and plays diverse functions in plant growth and development and environmental stress responses, especially in ABA signaling related to stress responses [27]. Out of 136 bZIP TFs identified in B. rapa (BrbZIP) [28], 132 bZIP TFs were found in Chinese cabbage transcriptome libraries, and exhibited differential expression patterns ( Figure 3A). Among them, 64 bZIP TFs were mainly differentially expressed in drought-treated plants. A total of 17 bZIP TFs were increased in both leaves and roots of drought-treated plants, and five bZIP TFs were only up-regulated in leaves. In addition, seven bZIP TFs were only down-regulated in leaves, six bZIP TFs only in roots, and three bZIP TFs were down-regulated in both roots and leaves. Drought stress-responsive gene expression is regulated by TFs operating in ABA-dependent and ABA-independent signaling pathways. Among these TFs, members of the bZIP TF family regulate the transcriptional induction of genes directly or indirectly linked with stress tolerance in plants [29]. In Arabidopsis, accumulated evidence has suggested that group-A bZIP including AtbZIP35, AtbZIP36, AtbZIP37, AtbZIP38, AtbZIP39, AtbZIP40, and AtbZIP66 play crucial roles in activating plant ABA signaling or in abiotic stress response [30,31]. We also found that BrbZIP23 (XP_009112062.1), BrbZIP15 (XP_009101738.1), BrbZIP17 (XP_009123528.1), BrbZIP21 (XP_009135624.1), BrbZIP22 (XP_009145807.1), and BrbZIP10 (XP_009104023.1), belonging to group-A BrbZIPs, were involved in drought response ( Figure 3B). In addition, BrbZIP101 (XP_009111876.1), BrbZIP114 (XP_009138305.1), and BrbZIP115 (XP_009111076.1), belonging to group-S bZIPs, were also found to be induced under drought conditions. Although the physiological function of group-S bZIPs remains poorly understood, the increased transcription levels of group-S bZIPs due to abiotic stresses might suggest their regulatory roles in stress response in plants [32][33][34][35]. Further, poplar PtabZIP1-like (PtabZIP1L), homologous to the Arabidopsis bZIP1 belonging to group-S bZIPs, has been known to be a positive modulator of drought resistance acting through modulation of multiple metabolic pathways including flavonoid biosynthesis [36]. These findings indicate that drought-induced group-S BrbZIPs ( Figure 3C) might be interesting candidate genes for improving drought tolerance.
Molecules 2018, 23, x FOR PEER REVIEW 6 of 15 bZIPs remains poorly understood, the increased transcription levels of group-S bZIPs due to abiotic stresses might suggest their regulatory roles in stress response in plants [32][33][34][35]. Further, poplar PtabZIP1-like (PtabZIP1L), homologous to the Arabidopsis bZIP1 belonging to group-S bZIPs, has been known to be a positive modulator of drought resistance acting through modulation of multiple metabolic pathways including flavonoid biosynthesis [36]. These findings indicate that droughtinduced group-S BrbZIPs ( Figure 3C) might be interesting candidate genes for improving drought tolerance. To validate the gene expression results from RNA-seq, the expression of nine bZIP TFs upregulated in both roots and leaves was analyzed using qRT-PCR. As shown in Figure 3B, the expression pattern of these genes was found to vary in response to drought. In addition, exposure of Chinese cabbage to NaCl or wounding also significantly induced or repressed the expression of these genes, although the changes were relatively lower than those due to drought treatment ( Figure S5). Further, selected bZIP TFs were found to show complex expression patterns in response to drought, NaCl, and wounding, as evidenced by cluster analysis results represented as a heatmap ( Figure 3C). In cluster I, BrbZIP (XP_009111876.1), also known as BrbZIP101 [26] exhibited relative higher transcript accumulation in drought-treated leaves and was transiently down-regulated in response To validate the gene expression results from RNA-seq, the expression of nine bZIP TFs up-regulated in both roots and leaves was analyzed using qRT-PCR. As shown in Figure 3B, the expression pattern of these genes was found to vary in response to drought. In addition, exposure of Chinese cabbage to NaCl or wounding also significantly induced or repressed the expression of these genes, although the changes were relatively lower than those due to drought treatment ( Figure S5). Further, selected bZIP TFs were found to show complex expression patterns in response to drought, NaCl, and wounding, as evidenced by cluster analysis results represented as a heatmap ( Figure 3C). In cluster I, BrbZIP (XP_009111876.1), also known as BrbZIP101 [26] exhibited relative higher transcript accumulation in drought-treated leaves and was transiently down-regulated in response to NaCl and wounding stresses, whereas BrbZIPs in cluster II only showed transient down-regulation due to wounding stress. These results indicate that BrbZIPs have divergent functions in response to environmental stresses, and the increased expression of BrbZIP101 might suggest its role as a marker of drought-stress response.
Functional enrichment analysis with respect to known gene ontologies using MapMan can be performed to obtain a global overview of high-throughput data in the context of metabolic pathways and biological processes [37]. To investigate the metabolic pathways affected by drought stress, MapMan analysis was performed on L2/L1, L3/L1, R2/R1, and R3/R1 comparisons based on DEGs. In leaves and roots of drought-treated plants, several metabolic pathways including those involving cell wall synthesis, light reactions, and the Calvin cycle were found to be down-regulated, whereas those involving fermentation, starch synthesis, and amino acid synthesis were up-regulated ( Figure 4). However, the impact of drought stress on the transcript levels of genes involved in these metabolic pathways in roots was higher than that in leaves. Noticeably, the genes involved in the sulfur-containing secondary metabolites (S-misc) category pathways were up-regulated in the leaves but down-regulated in the roots under drought-stress conditions, indicating differences between organs of Chinese cabbage in response to drought stress.

Alteration in Glucosinolate Metabolism in Rresponse to Drought Stress
In higher plants, a variety of sulfur-containing secondary metabolites are synthesized, which often play a crucial role in the survival of plants under biotic and abiotic stresses [38]. The alteration in sulfur metabolism under drought conditions is consistent with the promotion of root growth [39]. The down-regulation of glucosinolate metabolism in roots upon drought suggests a switch in sulfur metabolism from glucosinolate to other metabolic processes, which are required for enhancing root growth. Interestingly, we found that the genes involved in sulfur-containing secondary metabolites have opposing expression patterns in leaves and roots during response to drought stress ( Figure 4). With respect to the S-misc BIN category, the genes mainly represented were those responsible for glucosinolate biosynthesis including branched-chain amino acid aminotransferase 4 (BCAT4) and methylthioalkylmalate synthase 1 (MAMI1) ( Figure S6, Table S2). To test whether the increased expression of these genes correlated with the accumulation of glucosinolates, we analyzed drought-induced variation in glucosinolate composition and content. As shown in Figure 5, the total glucosinolate content was increased in response to drought stress in leaves of Chinese cabbage, suggesting that the accumulation of glucosinolates in L3 was due to increase transcription of glucosinolate biosynthesis genes.    [40,41]. Similarly, the accumulation of glucosinolates has been observed in Brassica species as a part of the plant's response to drought, although this response was highly variable depending on the duration and intensity of drought [40]. These findings indicate that drought stress significantly influenced the glucosinolate composition and content through modulating the expression of a range of glucosinolate biosynthesis genes.
Glucosinolate degradation products generated by the enzymatic action of myrosinases have attracted attention for several years due to their important role in plant defense against pests and pathogens [42]. In addition, the multiple roles of glucosinolates and their degradation products in heat tolerance, water transport, and transcriptional reprogramming suggest the physiological importance of these compounds in plant response to environmental stresses [43][44][45]. Exogenous application of glucosinolate degradation products, isothiocyanates, nitriles, and thiocyanates has been found to induce stomatal closure, allowing suppression of water loss [46,47]. Here, we hypothesized that drought-induced accumulation of glucosinolates may regulate stomatal closure leading to suppression of water loss. To test this hypothesis, glucosinolate-derived products were extracted from L1 and L3 samples and were applied to leaf disks floated on the opening medium. Application of glucosinolate-derived products significantly reduced stomatal apertures compared with application of mock-control (solvent control) ( Figure 6). In addition, variation in the levels of identified glucosinolate compounds between non-stressed plants (L1) and drought-treated plants (L3) was observed. The compound 4-methoxyglucobrassicin (0.73 ± 0.01 nmol/mg of D.W) was found to be the major constituent in Chinese cabbage cultivars (Chunkwang), followed by progoitrin (0.55 ± 0.01 nmol/mg of D.W) and glucobrassicin (0.42 ± 0.001 nmol/mg of D.W), whereas glucobrassicanapin (1.50 ± 0.02 nmol/mg of D.W) and 4-methoxyglucobrassicin (1.02 ± 0.002 nmol/mg of D.W) were the abundant glucosinolates in L3. Further, drought-stress treatment led to an approximately six-fold increase in levels of gluconapin and 4-hydroxyglucobrassicin, indicating that drought stress increased the transcription level of genes involved in the glucosinolate biosynthetic pathway ( Figure S6) and glucosinolate accumulation in the leaves of Chinese cabbage ( Figure 5), consistent with the protein competition model which predicts the effects of environmental factors on secondary metabolism [40,41]. Similarly, the accumulation of glucosinolates has been observed in Brassica species as a part of the plant's response to drought, although this response was highly variable depending on the duration and intensity of drought [40]. These findings indicate that drought stress significantly influenced the glucosinolate composition and content through modulating the expression of a range of glucosinolate biosynthesis genes.
Glucosinolate degradation products generated by the enzymatic action of myrosinases have attracted attention for several years due to their important role in plant defense against pests and pathogens [42]. In addition, the multiple roles of glucosinolates and their degradation products in heat tolerance, water transport, and transcriptional reprogramming suggest the physiological importance of these compounds in plant response to environmental stresses [43][44][45]. Exogenous application of glucosinolate degradation products, isothiocyanates, nitriles, and thiocyanates has been found to induce stomatal closure, allowing suppression of water loss [46,47]. Here, we hypothesized that drought-induced accumulation of glucosinolates may regulate stomatal closure leading to suppression of water loss. To test this hypothesis, glucosinolate-derived products were extracted from L1 and L3 samples and were applied to leaf disks floated on the opening medium. Application of glucosinolate-derived products significantly reduced stomatal apertures compared with application of mock-control (solvent control) ( Figure 6). The movement of the stomatal guard cells was found to be more responsive to glucosinolatederived products obtained from L3 than to those obtained from L1 (Figure 6), indicating that the accumulation of glucosinolates induced drought tolerance through stimulating stomatal closure. Taken together, the increased concentration of sulfate transported from the roots to the leaves might be required for enhancing the anti-transpiration effect of ABA in the stomata in leaves [48,49], and glucosinolate degradation products function analogously to ABA in stomatal aperture regulation [50]. These findings suggest that sulfur metabolism plays significant roles in drought stress signaling and responses [51] and indicate that variation in glucosinolate metabolism between roots and leaves ( Figure 5) is required for successful defense against drought conditions.

Plant Material and Stress Treatment
Seeds of Chinese cabbage (B. rapa L. subsp. pekinensis) cultivar Chunkwang were obtained from the Sakata Korea Seed Co. Ltd. (Seoul, South Korea), and were germinated and grown in a growth chamber with a long photoperiod (16 h light/8 h dark) at 24 °C. Soil water content was measured using HH2 Moisture Meter (Delta-T Devices, Cambridge, UK). The water content in the soil was controlled to be 60 to 65% of saturated soil water content. For drought stress treatment, 4-week-old plants (S1) were subjected to water withholding for 2 d or 4 d, resulting in a soil water content of 20% (S2) or 5% (S3), respectively. In addition, 4-week-old plants were watered with 250 mM NaCl or effectively wounded using sterile syringe needles to analyze the expression pattern of bZIP TFs. Each treatment was carried out three times, resulting in three biological replicates.

Determination of Relative Water Content, Malondialdehyde (MDA) Content, H2O2 Accumulation, and Protein Carbonyl Content
Relative water content was determined according to Soltys-Kalina et al. [52], and MDA content was determined using the thiobarbituric acid (TBA) reaction as described by Spanic et al. [53]. About 200 mg of ground material obtained from each sample (leaves from S1, S2, and S3) was used. MDA concentration was estimated by subtracting the non-specific turbidity at 600 nm from the absorption at 532 nm, and was expressed as MDA in nmol per mg of fresh weight using an absorbance coefficient of extinction (155 mM −1 cm −1 ).
As described by Spanic et al. [53], the supernatant obtained via 0.1% ice-cooled TCA extraction was used to analyze H2O2 concentration. After reaction with 1 M KI, the absorbance at 390 nm was The movement of the stomatal guard cells was found to be more responsive to glucosinolate-derived products obtained from L3 than to those obtained from L1 ( Figure 6), indicating that the accumulation of glucosinolates induced drought tolerance through stimulating stomatal closure. Taken together, the increased concentration of sulfate transported from the roots to the leaves might be required for enhancing the anti-transpiration effect of ABA in the stomata in leaves [48,49], and glucosinolate degradation products function analogously to ABA in stomatal aperture regulation [50]. These findings suggest that sulfur metabolism plays significant roles in drought stress signaling and responses [51] and indicate that variation in glucosinolate metabolism between roots and leaves ( Figure 5) is required for successful defense against drought conditions.

Plant Material and Stress Treatment
Seeds of Chinese cabbage (B. rapa L. subsp. pekinensis) cultivar Chunkwang were obtained from the Sakata Korea Seed Co. Ltd. (Seoul, South Korea), and were germinated and grown in a growth chamber with a long photoperiod (16 h light/8 h dark) at 24 • C. Soil water content was measured using HH2 Moisture Meter (Delta-T Devices, Cambridge, UK). The water content in the soil was controlled to be 60 to 65% of saturated soil water content. For drought stress treatment, 4-week-old plants (S1) were subjected to water withholding for 2 d or 4 d, resulting in a soil water content of 20% (S2) or 5% (S3), respectively. In addition, 4-week-old plants were watered with 250 mM NaCl or effectively wounded using sterile syringe needles to analyze the expression pattern of bZIP TFs. Each treatment was carried out three times, resulting in three biological replicates.

Determination of Relative Water Content, Malondialdehyde (MDA) Content, H 2 O 2 Accumulation, and Protein Carbonyl Content
Relative water content was determined according to Soltys-Kalina et al. [52], and MDA content was determined using the thiobarbituric acid (TBA) reaction as described by Spanic et al. [53]. About 200 mg of ground material obtained from each sample (leaves from S1, S2, and S3) was used. MDA concentration was estimated by subtracting the non-specific turbidity at 600 nm from the absorption at 532 nm, and was expressed as MDA in nmol per mg of fresh weight using an absorbance coefficient of extinction (155 mM −1 cm −1 ).
As described by Spanic et al. [53], the supernatant obtained via 0.1% ice-cooled TCA extraction was used to analyze H 2 O 2 concentration. After reaction with 1 M KI, the absorbance at 390 nm was determined via comparison to a standard curve prepared using H 2 O 2 . H 2 O 2 concentration in each sample was expressed as µM per mg of fresh weight.
For protein carbonyl determination, 300 mg of frozen dried samples were incubated with protein extraction buffer containing 50 mM Tris-HCl, pH 7.5, 1 mM EDTA, 5 mM MgCl 2 , and 10 mM β-mercaptoethanol. After centrifugation at 12,000 g for 15 min at 4 • C, the protein content in the supernatant was determined following the Bradford method with bovine serum albumin as standard. Protein carbonyl content was determined using a fluorometric protein carbonyl content assay kit (BioVision, Milpitas, CA, USA) following the manufacturer's instructions.

Identification and Functional Annotation of DEGs
Total RNA was extracted from each sample using the FavorPrep Plant Total RNA Mini Kit (Favorgen Biotech Corp, Pingtung, Taiwan), and quantified using an DeNovix DS-11 (DeNovix Inc., Wilmington, DE, USA). Equal amounts of total RNA from different stages of each sample were pooled to generate a cDNA library. Equal amounts of total RNA from different stages of each sample were pooled to generate a cDNA library. After DNaseI treatment, poly(A) mRNA was isolated with oligo(dT) magnetic beads and fragmented to short fragments. Then, the cDNA library was synthesized as described by Choi et al. [54] and sequenced using an Illumina HiSeq™ 2500 sequencing platform. The quality of raw data was controlled by using the FastQC tool, and adapter sequences, empty reads, low-quality reads (with ambiguous sequence, N), and reads with more than 10% Q < 20 bases (i.e., with a base quality of less than 20) were discarded using Trimmomatic v.0.33. The clean reads were mapped to the B. rapa reference sequence (BrapaV2.1, http://brassicadb.org/brad/datasets/ pub/Genomes/Brassica_rapa/V2.0/V2.1/Chr/) using the HISAT2 aligner [55]. Transcript levels were calculated using SAMtools (http://samtools.sourceforge.net/) and relative transcript abundances were analyzed using DEseq [56]. DEGs were determined by combining a p value cutoff of 0.01 and adjusting to |log2 (fold change)| ≥ 1.For functional annotation, gene ontology (GO) enrichment analysis of DEGs was performed using Blast2GO (v.4.1.0) In addition, MapMan metabolism overview maps were drawn as described by Yu et al. [57].

Quantitative Real-Time PCR (qRT-PCR) Analysis
To further verify the expression profiles observed in the RNA-Seq data, nine bZIP TFs were selected for qRT-PCR analysis performed using the SYBR ® Green Real-time PCR Master Mix (TOYOBO, Co., Ltd., Osaka, Japan) in the CFX96TM Real-time system (Bio-Rad, Hercules, CA, USA) with default parameters. The conditions for all reactions were as follows: 95 • C for 60 s, 40 cycles of 95 • C for 5 s, followed by 60 • C for 15 s, and 72 • C for 20 s. Melting curve analysis was performed to confirm the PCR specificity. The relative expression levels of selected bZIP TFs were normalized to expression of the internal reference gene actin [58]. Specific primer pairs used in qRT-PCR are listed in Table S3.

Measurement of Stomatal Aperture
Stomatal apertures were measured as described by Hossain et al. [47]. Briefly, leaf disks were floated on the opening medium (5 mM KCl, 50 µM CaCl 2 and 10 mM MES-Tris, pH 6.15) for 2 h under light conditions. To analyze the effect of glucosinolate-derived products on stomatal movement, glucosinolate-derived products was extracted from 50 mg of frozen dried samples (leaves from S1 and S3) as described by Lola-Luz et al. [61]. Each extract was re-suspended in 5 mL HPLC-grade methanol and filtered using a membrane filter (0.2 µm pore size). Then, leaf disks floated on the opening medium were treated with 5% of extract (50 µL of each extract suspended in 5 mL MeOH was mixed with 1 mL of opening medium). Two hours after treatment, stomatal apertures were captured using an CKX53 microscope (Olympus, Tokyo, Japan) and the width and length of 100 stomatal apertures were measured using the image processing software ImageJ (National Institute of Mental Health, Bethesda, MD, USA).

Conclusions
Because of the complexity of drought as a stress signal, understanding the biochemical and molecular mechanisms of plant response has remained a major challenge in plant biology and breeding. The present study provides an overview of the changes in transcript abundance in leaves and roots under drought conditions. DEGs found to be responsive to drought suggest differential regulation patterns in the leaf and root, and drought-induced group-S BrbZIPs are interesting candidates for improving plant tolerance. Further, we suggest that organ-specific metabolic responses to drought are part of the survival strategies of plants. Taken together, our dataset provides a crucial starting point for future efforts in understanding the biochemical and molecular mechanisms underlying response to drought stress in Chinese cabbage.