Physiological and Transcriptional Responses of Industrial Rapeseed (Brassica napus) Seedlings to Drought and Salinity Stress

Abiotic stress greatly inhibits crop growth and reduces yields. However, little is known about the transcriptomic changes that occur in the industrial oilseed crop, rapeseed (Brassica napus), in response to abiotic stress. In this study, we examined the physiological and transcriptional responses of rapeseed to drought (simulated by treatment with 15% (w/v) polyethylene glycol (PEG) 6000) and salinity (150 mM NaCl) stress. Proline contents in young seedlings greatly increased under both conditions after 3 h of treatment, whereas the levels of antioxidant enzymes remained unchanged. We assembled transcripts from the leaves and roots of rapeseed and performed BLASTN searches against the rapeseed genome database for the first time. Gene ontology analysis indicated that DEGs involved in catalytic activity, metabolic process, and response to stimulus were highly enriched. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that differentially expressed genes (DEGs) from the categories metabolic pathways and biosynthesis of secondary metabolites were highly enriched. We determined that myeloblastosis (MYB), NAM/ATAF1-2/CUC2 (NAC), and APETALA2/ethylene-responsive element binding proteins (AP2-EREBP) transcription factors function as major switches that control downstream gene expression and that proline plays a role under short-term abiotic stress treatment due to increased expression of synthesis and decreased expression of degradation. Furthermore, many common genes function in the response to both types of stress in this rapeseed.


Introduction
Throughout their lifecycles, plants are subjected to various external environmental stresses including biotic stresses (such as weeds and diseases) and abiotic stresses (such as drought and salinity). In general, abiotic stress reduces crop yields by more than 50% compared to the less than 10% reduction caused by biotic stress [1]. Drought affects approximately 40% of the world's agricultural land and is considered to be the most serious global agricultural problem [2]. Moreover, it is estimated that 800 million hectares of land have been salinized to some extent [3]. The limiting effect of drought and salinity on plant growth is expected to increase due to global climate change.
Plants adapt to short periods of abiotic stress through physiological regulation. Plant survival increases through osmotic adjustment or the removal of reactive oxygen species (ROS). Proline and sugars play important roles in osmotic adjustment, while superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) mediate ROS removal. In addition to physiological regulation, external stress Int. J. Mol. Sci. 2019, 20, 5604 2 of 21 signals are transduced in plants that induce stimulus-specific changes in gene expression. Plants have several common stress-related genes that allow them to withstand different adverse conditions [1,4].
Rapeseed (Brassica napus) is the third most important oilseed crop after soybean (Glycine max) and palm (Trachycarpus sortunei) [5]. Rapeseed cultivars with low glucosinolate and erucic acid contents are used to produce edible oils, animal feed, and biodiesel, whereas rapeseed cultivars with high erucic acid content have industrial applications. Since drought and salinity severely affect yield in rapeseed, changes in gene expression under different abiotic stress treatments have been studied to identify stress-specific genes for use in genetic engineering. For example, transcriptomic differences in two rapeseed varieties with different levels of drought resistance were identified at seven days after germination [6]. Differentially expressed genes (DEGs) were also identified in the roots of rapeseed seedlings under salinity stress soon after germination [7]. Transcriptomic and epigenomic comparisons have revealed that the expression levels of several drought-responsive genes in rapeseed seedlings depend on methylation patterns in the genome [8].
Transcriptome analysis is a powerful tool for identifying stress-related genes due to the ability to fully analyze gene transcription. Rapeseed (AACC genome) originated from a hybridization between Brassica rapa (AA genome) and Brassica oleracea (CC genome) [9]. However, even though sequencing of the rapeseed genome was completed in 2014 [10], most transcriptome analyses performed to date have been based on the NCBI non-redundant (NR) database [6], B. rapa genome database [7], or both the B. rapa and B. oleracea genomes [8]. Although rapeseed is closely related to other Brassica species, extensive genomic differences exist in the Brassicaceae family [10]. There are several reports of the transcriptome analysis of B. napus under different abiotic stresses such as freezing [11], silicon supply [12], metal stress [13], drought [14], and salt [15]. However, no common abiotic stress-related genes have been identified through the transcriptome analysis of rapeseed under both drought and salinity stress. After elaborating the mechanism behind abiotic stresses, new agronomic and breeding practices could be better applied to increase its stress adaption [16].
In the current study, we analyzed the transcriptome of an industrially important salinity-resistant rapeseed variety with high erucic acid content (Nanyanyou-1; resistant to 0.2-0.3% salinity) under drought (simulated by treatment with 15% (w/v) PEG 6000) and salinity (150 mM NaCl) stress. Our results indicated that proline metabolism played an important role in rapeseed under short-term abiotic stress (3 h). More DEGs were identified under drought vs. salinity stress and in the roots vs. leaf tissue. The findings of this study lay the foundation for developing industrially important rapeseed varieties with tolerance to drought and salinity stress using genetic engineering techniques and effective breeding programs.

Physiological Characteristics of B. napus under Abiotic Stresses
Long-term drought and salinity stress significantly affect the physiological and biochemical functioning of plants. However, there were no obvious differences in superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), H 2 O 2 , or soluble sugar contents in rapeseed seedlings after 3 h of treatment with 150 mM NaCl or 15% (w/v) PEG 6000 ( Figure 1A-E). Malondialdehyde (MDA) and relative water contents were slightly different in the drought and salinity-treated seedlings when compared to the untreated control ( Figure 1F,G), but no difference was detected between the two abiotic treatments. In contrast, both stress treatments triggered a marked increase in proline content ( Figure 1H), which increased by more than 50% under NaCl treatment and was doubled under polyethylene glycol (PEG) treatment. , and proline-PRO (H) contents were determined after 3 h of drought (simulated by polyethylene glycol (PEG)) or salinity (NaCl) treatment. CK represents the control group, (i.e., seedlings treated with ½ Hoagland nutrient solution); PEG represents seedlings treated with 15% (w/v) PEG 6000 plus ½ Hoagland nutrient solution for 3 h. NaCl represents seedlings treated with 150 mM NaCl plus ½ Hoagland nutrient solution for 3 h. Each data point represents the mean of three samples ± SE. Columns with different letters in each graph indicate significant differences based on Duncan's multiple range tests at p < 0.05 among treatments. , and proline-PRO (H) contents were determined after 3 h of drought (simulated by polyethylene glycol (PEG)) or salinity (NaCl) treatment. CK represents the control group, (i.e., seedlings treated with 1 2 Hoagland nutrient solution); PEG represents seedlings treated with 15% (w/v) PEG 6000 plus 1 2 Hoagland nutrient solution for 3 h. NaCl represents seedlings treated with 150 mM NaCl plus 1 2 Hoagland nutrient solution for 3 h. Each data point represents the mean of three samples ± SE. Columns with different letters in each graph indicate significant differences based on Duncan's multiple range tests at p < 0.05 among treatments.

Transcriptomic Sequencing and Unigene Assembly
To further explore the molecular mechanism underlying the physiological responses of rapeseed to drought and salinity stress, the roots and leaves of stress-treated plants were collected, and 18 samples were used to construct transcriptome libraries and sequenced on the Illumina HiSeq platform.
On average, approximately 6.64 Gb were generated per sample (Table 1). After mapping sequenced reads to the B. napus reference genome and reconstructing the transcripts using HISAT, an average of 66.27% mapped reads had met the requirements for further comparison. Ultimately, 42,251 transcripts were obtained from all samples including 4251 transcripts of novel genes that contain features not present in the reference annotation and 7127 long noncoding RNAs (Table 1). The value of Q30 exceeded 95%, suggesting that the sequencing data were reliable.

Analysis of Differentially Expressed Genes
DEGs were defined based on a fold change ≥ 2.00 and adjusted p-value ≤ 0.05. We  (Table S1, fold change ≥ 4.00 also listed here). There were more DEGs under drought treatment than under high-salinity treatment and many more DEGs in roots compared to leaves ( Figure 2). More DEGs were upregulated in roots under drought stress than in roots under salinity stress, but more DEGs were downregulated in roots under salinity stress than in roots under drought stress, highlighting the different expression patterns of DEGs under these two abiotic stresses.
As illustrated in Venn diagrams (Figure 2), 102 DEGs were upregulated and 135 DEGs were downregulated in both leaves and roots under PEG treatment. In contrast, only four DEGs were upregulated and one DEG was downregulated in both leaves and roots under high-salinity treatment. No common DEGs were detected in all tissues under both treatments ( Figure S1). We performed hierarchical clustering analysis based on log2 fold change values to compare the expression patterns of upregulated and downregulated genes between groups ( Figure S2). Basically, the transcriptional profiles between leaves and roots under salinity treatment were similar and also similar to that of roots under drought treatment to the extent. In both roots and leaves, DEGs in the drought-stressed group exhibited greater differences in expression relative to the untreated control than those in the salinity-treated group. In addition, most DEGs in roots showed greater differences in expression relative to the untreated control group than that of DEGs in leaves under both treatments.

GO and KEGG Pathway Enrichment Analysis
Gene ontology (GO) is a common functional classification system used to identify genes with specific functions in three categories: cellular component, molecular function, and biology process. Ultimately, 812 out of 913 DEGs in drought-stressed leaves, 3366 out of 3879 DEGs in drought-stressed roots, 82 out of 95 DEGs in salinity-stressed leaves, and 524 out of 616 DEGs in salinity-stressed roots were assigned to GO categories ( Figure 3). Most genes in each category were assigned to one of three sub-categories (cell part, cell, and organelle from cellular component; catalytic activity, binding, and transporter activity from molecular function; metabolic process, cellular process, and response to stimulus from biological process).  Figure S1).
We performed hierarchical clustering analysis based on log 2 fold change values to compare the expression patterns of upregulated and downregulated genes between groups ( Figure S2). Basically, the transcriptional profiles between leaves and roots under salinity treatment were similar and also similar to that of roots under drought treatment to the extent. In both roots and leaves, DEGs in the drought-stressed group exhibited greater differences in expression relative to the untreated control than those in the salinity-treated group. In addition, most DEGs in roots showed greater differences in expression relative to the untreated control group than that of DEGs in leaves under both treatments.

GO and KEGG Pathway Enrichment Analysis
Gene ontology (GO) is a common functional classification system used to identify genes with specific functions in three categories: cellular component, molecular function, and biology process. Ultimately, 812 out of 913 DEGs in drought-stressed leaves, 3366 out of 3879 DEGs in drought-stressed roots, 82 out of 95 DEGs in salinity-stressed leaves, and 524 out of 616 DEGs in salinity-stressed roots were assigned to GO categories ( Figure 3). Most genes in each category were assigned to one of three sub-categories (cell part, cell, and organelle from cellular component; catalytic activity, binding, and transporter activity from molecular function; metabolic process, cellular process, and response to stimulus from biological process).   To conduct a comprehensive pathway analysis of the DEGs in each group, we performed a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Due to the limitations of the KEGG database, only 80.7% of the DEGs were assigned to KEGG pathways. KEGG pathways were assigned based on a threshold value of p ≤ 0.02. There were 24 significantly enriched pathways for DEGs in CK-L vs. D-L and 39, 3, and 16 pathways for CK-L vs. D-R, CK-R vs. S-L, and CK-L vs. S-R, respectively. Under PEG treatment, DEGs involved in metabolic pathways and biosynthesis of secondary metabolites were highly enriched (Table 2A,B). Under NaCl treatment, the same pathways were enriched in roots as those enriched in roots under drought treatment, while the cutin, suberin, and wax biosynthesis pathways and ABC transporters were markedly enriched in leaves under salinity treatment (Tables 2C,D). To conduct a comprehensive pathway analysis of the DEGs in each group, we performed a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Due to the limitations of the KEGG database, only 80.7% of the DEGs were assigned to KEGG pathways. KEGG pathways were assigned based on a threshold value of p ≤ 0.02. There were 24 significantly enriched pathways for DEGs in CK-L vs. D-L and 39, 3, and 16 pathways for CK-L vs. D-R, CK-R vs. S-L, and CK-L vs. S-R, respectively. Under PEG treatment, DEGs involved in metabolic pathways and biosynthesis of secondary metabolites were highly enriched (Table 2A,B). Under NaCl treatment, the same pathways were enriched in roots as those enriched in roots under drought treatment, while the cutin, suberin, and wax biosynthesis pathways and ABC transporters were markedly enriched in leaves under salinity treatment (Table 2C,D). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using DEGs from four treatment groups. CK-L (CK, leaves) vs. D-L (drought, leaves), CK-R (CK, roots) vs. D-R (drought, roots), CK-L vs. S-L (salinity, leaves), CK-R vs. S-R (salinity, roots) represent DEGs under these two abiotic stresses. Bolded rows represent the most highly enriched metabolic pathways.

Analysis of Differentially Expressed Transcription Factors
To investigate the importance of transcription factors (TFs) in the stress resistance regulatory network, we analyzed differentially expressed TF genes. We identified 78, 587, 3, and 39 TF genes in CK-L vs. D-L, CK-R vs. D-R, CK-L vs. S-L, and CK-R vs. S-R, respectively. Both the numbers and classifications of TF genes differed in roots and leaves under both treatments (Figure 4). Under PEG treatment, the basic region-leucine zipper (bZIP), myeloblastosis (MYB), plant AT-rich sequence-and zinc-binding (PLATZ), and basic/helix-loop-helix (bHLH) TF family genes were enriched in leaves, while the MYB, NAC, bHLH, APETALA2/ethylene-responsive element binding proteins (AP2-EREBP), and WRKY domain contained proteins (WRKY) TF family genes were enriched in roots. Under NaCl treatment, AP2-EREBP and C2C2-CO-like TF genes were enriched in leaves, while NAM/ATAF1-2/CUC2 (NAC), MYB, WRKY, and G2-like TF genes were enriched in roots. Moreover, 435 out of 665 TF genes were upregulated under PEG treatment, whereas similar numbers of TF genes were up-or downregulated under high-salinity treatment (

Analysis of Differentially Expressed Transcription Factors
To investigate the importance of transcription factors (TFs) in the stress resistance regulatory network, we analyzed differentially expressed TF genes. We identified 78, 587, 3, and 39 TF genes in CK-L vs. D-L, CK-R vs. D-R, CK-L vs. S-L, and CK-R vs. S-R, respectively. Both the numbers and classifications of TF genes differed in roots and leaves under both treatments (Figure 4). Under PEG treatment, the basic region-leucine zipper (bZIP), myeloblastosis (MYB), plant AT-rich sequenceand zinc-binding (PLATZ), and basic/helix-loop-helix (bHLH) TF family genes were enriched in leaves, while the MYB, NAC, bHLH, APETALA2/ethylene-responsive element binding proteins (AP2-EREBP), and WRKY domain contained proteins (WRKY) TF family genes were enriched in roots. Under NaCl treatment, AP2-EREBP and C2C2-CO-like TF genes were enriched in leaves, while NAM/ATAF1-2/CUC2 (NAC), MYB, WRKY, and G2-like TF genes were enriched in roots. Moreover, 435 out of 665 TF genes were upregulated under PEG treatment, whereas similar numbers of TF genes were up-or downregulated under high-salinity treatment (Table S2).
As far as proline metabolism, the promoters of DEGs encoding BnOAT (BnaA06g36140D), BnP5CS1 (BnaA03g18760D, BnaA05g05760D, BnaC03g72600D, BnaC04g05620D), and BnP5CS2 (BnaC04g55570D, BnaA04g03460D) were further analyzed (Table S3). The analysis of putative cis-regulatory elements showed that all of these DEGs harbored binding sites for RAV, NAC, bZIP, and MYB. Except for BnaA06g36140D, the other members harbored binding sites for WRKY. Binding sites for AP2-ERF only existed in the promoters of DEGs encoding BnP5CS1.

Validation of Differentially Expressed Genes by Quantitative RT-PCR Analysis
To further verify the reliability of the sequencing data, we randomly selected 10 DEGs (five under PEG treatment, five under NaCl treatment) and subjected them to qRT-PCR analysis. Although the fold change values of these DEGs were not the same as those obtained by sequencing, the trends were consistent ( Figure 5 and Table S4). As far as proline metabolism, the promoters of DEGs encoding BnOAT (BnaA06g36140D), BnP5CS1 (BnaA03g18760D, BnaA05g05760D, BnaC03g72600D, BnaC04g05620D), and BnP5CS2 (BnaC04g55570D, BnaA04g03460D) were further analyzed (Table S3). The analysis of putative cis-regulatory elements showed that all of these DEGs harbored binding sites for RAV, NAC, bZIP, and MYB. Except for BnaA06g36140D, the other members harbored binding sites for WRKY. Binding sites for AP2-ERF only existed in the promoters of DEGs encoding BnP5CS1.

Validation of Differentially Expressed Genes by Quantitative RT-PCR Analysis
To further verify the reliability of the sequencing data, we randomly selected 10 DEGs (five under PEG treatment, five under NaCl treatment) and subjected them to qRT-PCR analysis. Although the fold change values of these DEGs were not the same as those obtained by sequencing, the trends were consistent ( Figure 5 and Table S4).
In Brassica napus, proline levels sharply increase under both drought [26] and high salinity stress [27]. Indeed, in the current study, both drought and salinity stress led to the significant accumulation of proline in rapeseed ( Figure 1H). Therefore, we further analyzed DEGs involved in proline metabolism. We detected two DEGs that were upregulated in leaves under PEG treatment, five DEGs (one upregulated and four downregulated) in roots under drought treatment, and six DEGs (three upregulated and three downregulated) in roots under NaCl treatment (Table S3).
Both the biosynthesis and degradation of proline play important roles in proline accumulation. Pathways involved in proline biosynthesis include the glutamate (Glu) and ornithine (Orn) pathways. In the Glu pathway, glutamate is converted to GSA (glutamate-γ-semialdehyde) via a reaction catalyzed by P5CS (∆ 1 -pyrroline-5-carboxylate synthetase). GSA is further reduced to proline by P5CR (∆ 1 -pyrroline-5-carboxylate reductase). When proline accumulates to a certain level, the expression of P5CS is suppressed, while that of ProDH is induced, leading to proline degradation. In the Orn pathway, ornithine loses its δ-amino group to generate GSA through transamination, a process mediated by OAT (ornithine aminotransferase). GSA participates in the Glu pathway to generate proline. Therefore, the key enzymes in the proline biosynthesis pathway are P5CS, P5CR, and δ-OAT, whereas the key enzyme in the proline degradation pathway is ProDH (proline dehydrogenase) [28,29].
In the current study, only DEGs in the Glu pathway were identified in CK-R vs. S-R and CK-R vs. D-R, whereas DEGs in both pathways were induced in CK-L vs. D-L. No DEGs in the Glu or Orn pathway were enriched in CK-L vs. S-L. Genes encoding BnP5CS proteins in the Glu pathway were upregulated under stress treatment including BnP5CS1 (BnaA03g18760D, BnaA05g05760D, BnaC03g72600D, BnaC04g05620D), and BnP5CS2 (BnaC04g55570D, BnaA04g03460D) ( Figure 6). Similarly, BnaA06g36140D, a gene encoding BnOAT in the Orn pathway, was also upregulated in leaves under drought stress. Furthermore, we identified genes encoding ProDH that were downregulated in roots under both treatments including BnaC02g38230D, BnaA06g39660D, and BnaAnng07910D (Table S3). Thus, when treated with PEG or NaCl, rapeseed seedlings synthesize more proline by increasing the expression of P5CS or OAT, and they also reduce the degradation of proline by suppressing ProDH expression.
In Arabidopsis thaliana, which is a close relative of B. napus, P5CS is encoded by two similar regulatory genes named AtP5CS1 and AtP5CS2, while ProDH is encoded by AtPDH. When treated with ABA or NaCl, Arabidopsis plants accumulate increased levels of proline due to high AtP5CS1 expression and low AtPDH expression [30]. When treated with exogenous H 2 O 2 , proline accumulation increases in plants, which is associated with the activation of the Glu and Orn proline biosynthesis pathways [31]. AtP5CS1 and AtP5CS2 are induced by cold, NaCl, abscisic acid (ABA), desiccation, light, heat, rehydration, and brassinolide treatment [32]. In B. napus, proline accumulation during priming and post-priming germination is associated with the strong upregulation of P5CSA and the downregulation of PDH In Arabidopsis thaliana, which is a close relative of B. napus, P5CS is encoded by two similar regulatory genes named AtP5CS1 and AtP5CS2, while ProDH is encoded by AtPDH. When treated with ABA or NaCl, Arabidopsis plants accumulate increased levels of proline due to high AtP5CS1 expression and low AtPDH expression [30]. When treated with exogenous H2O2, proline accumulation increases in plants, which is associated with the activation of the Glu and Orn proline biosynthesis pathways [31]. AtP5CS1 and AtP5CS2 are induced by cold, NaCl, abscisic acid (ABA), desiccation, light, heat, rehydration, and brassinolide treatment [32]. In B. napus, proline accumulation during priming and post-priming germination is associated with the strong upregulation of P5CSA and the downregulation of PDH [33].
The proline pathway is regulated by several types of TFs. In Arabidopsis, AtP5CS1/2, AtP5CR, and AtOAT contain several TF binding sites 1000 bp upstream of their promoter regions that bind TFs such as MYB, bZIP, AP2/EREBP, WRKY, and RAV TFs [32]. In rice, genes in the Glu pathway are thought to be targeted by many TFs. For example, 24 different classes of TFs have binding sites in the promoters of OsP5CS1/2 and OsP5CR [34]. In Medicago truncatula, MtMYBS1, a MYB TF gene, The proline pathway is regulated by several types of TFs. In Arabidopsis, AtP5CS1/2, AtP5CR, and AtOAT contain several TF binding sites 1000 bp upstream of their promoter regions that bind TFs such as MYB, bZIP, AP2/EREBP, WRKY, and RAV TFs [32]. In rice, genes in the Glu pathway are thought to be targeted by many TFs. For example, 24 different classes of TFs have binding sites in the promoters of OsP5CS1/2 and OsP5CR [34]. In Medicago truncatula, MtMYBS1, a MYB TF gene, is induced by NaCl, PEG, and ABA treatment, and MtMYBS1 enhances the transcription of P5CS [35]. Transgenic Betula platyphylla plants overexpressing BplMYB46 exhibited improved salinity and osmotic tolerance due to the increased expression of P5CS genes [36]. In addition, AP2/EREBP upregulates AtP5CS1 in response to low water potential [37]. In rice, heterologous expression of JERF1, which encodes a tomato ERF protein, increased the accumulation of proline by upregulating OsP5CS [38]. Of course, even in the same TF family, different members could function differently. The analysis of phylogenetic trees showed that reported MtMYBS1 and BplMYB46 had similarities with BnMYBs including BnaA01g32800D and BnaA03g29470D, and BnaA05g30870D. JERF1 had similarities with BnaA07g30130D ( Figure S3). This implies that these TFs may have similar function in the regulation of expression of proline biosynthetic genes. In the current study, we also identified several binding sites for TFs such as MYB, WRKY, and bZIP in the 1000-bp upstream regions of the BnP5CS and BnOAT promoters ( Figure S4). The identified promoter region all had binding sites of MYB and most of them also had AP2-ERF binding sites ( Figure S4). Perhaps in response to abiotic stress, BnP5CSs and BnOAT are upregulated by those TFs in rapeseed to induce the accumulation of proline. Without doubt, the region far away those upstream 1000-bp or 3 untranslated region may have binding sites that also regulate the transcription of proline-related genes.

The Multiple Transcripts in Response to Salinity and Drought Stresses
Based on transcriptome analysis, more DEGs were present in the roots than in leaves after 3 h of abiotic stress treatment (Figure 2 and Figure S2); this result is consistent with previous findings [1]. Using microarray analysis, 624 DEGs were detected in Arabidopsis roots compared to only 285 in leaves in response to a 3-h salinity treatment. Compared to salinity treatment, more DEGs were detected in rapeseed following PEG treatment in the current study. Transcriptome analysis revealed 1700 DEGs in rapeseed after 3 h vs. 24 h of salinity treatment, most of which were downregulated in the whole roots [7]. These findings indicate that DEG responses to abiotic stress differ depending on the plant species, treatment, and growth stage.
Based on KEGG analysis, the most important DEGs in both roots and leaves under drought stress and in roots under salinity stress are involved in the 'metabolic pathway' and 'biosynthesis of secondary metabolites' (Table 2A,B,D). It is not surprising that many DEGs are involved in secondary metabolite biosynthesis due to the important roles of these compounds in plants subjected to stress [39][40][41]. Several ABC transporter genes were differentially expressed in leaves under salinity stress (Table 2C). ABC transporters (ATP-binding cassette transporters) transport organic materials, especially secondary metabolites, for plant development [42]. ABC transporters also play important roles in cutin and wax formation [43]. Several DEGs related to cutin, suberin, and wax biosynthesis were identified in leaves under salinity stress, confirming the importance of secondary metabolites in abiotic stress responses.
MYB and MYB-related TF genes were highly upregulated under both salinity and drought stress ( Figure 4). MYBs, which are characterized by a highly conserved DNA-binding MYB domain are present in all eukaryotes [44]. MYBs form a large gene family involved in development, metabolism, and stress responses [45,46]. MYBs also form homo-and heterodimers with other proteins to regulate various processes in plants [47].
In addition to MYB TF genes, NAC TF genes were highly upregulated by abiotic stress in rapeseed, especially salinity stress ( Figure 4B). Like MYB TFs, NAC TFs form a large family, but are specific to plants [48]. In general, the N-termini of NAC TFs are highly conserved regions that function in DNA binding, and protein-protein interactions. Increasing evidence suggests that NAC TFs are involved in regulating abiotic or biotic stress responses [49].
AP2/EREBP TFs contain DNA-binding AP2 domains, which are also unique to plants [50]. This large TF family participates in diverse stress responses [51]. Plants use abscisic acid (ABA)-dependent and ABA-independent pathways to cope with abiotic stress. MYB TFs function in ABA-dependent pathways, whereas AP2/EREBP TFs function in ABA-independent pathways. Many AP2/EREBP TFs such as CBFs and DREBs have been extensively characterized and manipulated to increase abiotic stress resistance in plants [52].
Numerous DEGs identified in the current study showed the same expression patterns under both abiotic stress treatments (Table S5). There were more DEGs shared in roots than in leaves under the two treatments. Most of these genes were categorized as encoding hypothetical proteins. Several genes involved in proline metabolism or abiotic stress pathways were induced in both the leaves and roots under abiotic stress, whereas those involved in primary metabolism tended to be downregulated under these treatments. These findings promote our understanding of the regulatory mechanisms underlying the drought and salinity stress responses in B. napus and lay the foundation for breeding cultivars with improved tolerance to abiotic stress.

Plant Materials and Growth Conditions
Brassica napus 'Nanyanyou-1' seeds were harvested in Nanjing, Jiangsu Province, China (31~32 • N, 118~119 • E) in 2016 and kept at Nanjing Agricultural University. The seeds were first germinated on wet gauze (soaked with water) in the incubator with a light intensity of 392~415 µmol m −2 s −1 during a daily cycle consisting of 16 h of light at 25 • C and 8 h of darkness at 18 • C. The seedlings were then transferred into a 1 2 Hoagland nutrient solution to conduct a hydroponic experiment under the same culture conditions for nearly 20 days until the fourth leaves had extended. The plants were treated with 1 2 Hoagland nutrient solution containing 15% (w/v) PEG 6000 or 150 mM NaCl for drought and salinity stress treatments, respectively; seedlings treated with 1 2 Hoagland nutrient solution alone were used as the control. Each treatment included three biological replications. Leaves and roots were harvested individually after 3 h of treatment, immediately frozen in liquid nitrogen, and stored at 80 • C until use for RNA extraction and physiological and biochemical analyses.

Physiological and Biochemical Analyses
For water content (WC) determination, the fresh weight (FW) of the whole seedlings was first measured and then dried to a constant weight at 65 • C for 72 h to obtain the dry weight (DW). WC was calculated as follows: WC (%) = (FW − DW)/FW × 100%. MDA (malondialdehyde), Pro (proline), SOD (superoxide dismutase), POD (peroxidase), CAT (catalase), H 2 O 2 , and soluble sugar contents were determined using microdetermination kits (Suzhou Comin Biotechnology Co. Ltd, Suzhou, China). One-way ANOVA using the Least Significant Difference (LSD) method was conducted with SPSS 19.0 software (SPSS Corp., Chicago, IL, USA) to evaluate significant differences among treatments. Figures were drawn using SigmaPlot 10 (Systat Software, Inc., Berlin, Germany) software.

RNA Extraction, cDNA Library Construction, and Sequencing
Eighteen rapeseed plant samples including roots and leaves were sent to Beijing Genomics Institute (BGI) for RNA extraction and cDNA library construction. Total RNA was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA, USA) for each biological replicate and was treated with DNase. The mRNA were purified from total RNA using magnetic Oligo (dT) beads and then fragmented by the fragmentation buffer. The cDNA was synthesized using the mRNA fragments as templates and the library was built. Agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System were used in the quantification and qualification of the sample library. Next, the library was sequenced using Illumina HiSeq 4000 to generate 150 bp, PE type data. The RNA-seq data was uploaded to the NCBI with the accession number PRJNA579479.

Identification and Functional Enrichment Analysis of DEGs
DEseq2 (version: v3.10, Berlin, Germany.) was used to identify DEGs based on the parameters of an adjusted fold change of ≥2.00 and p-value ≤ 0.05 [57]. DEG lists were uploaded and analyzed online, and Venn diagrams were constructed (http://bioinformatics.psb.ugent.be/webtools/Venn/). Hierarchical clustering of the DEGs was performed using MeV-4.9.0 software (http://sourceforge.net/ projects/mev-tm4/). The GO (Gene Ontology) annotation results were visualized, compared, and plotted using the online program WEGO (http://wego.genomics.org.cn) [58]. We also performed KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway classification for functional enrichment of the DEGs [59]. For p-value correction (to control the false discovery rate, FDR), the rigorous Bonferroni correction method was used [60]. DEGs with a corrected p-value and FDR of ≤ 0.001 were defined as significantly enriched for both types of functional enrichment analyses.

Analysis of Putative Cis-Regulatory Elements Involved in Proline Metabolism
The upstream 1000 bp genomic sequence of BnP5CS and BnOAT before the ATG codon were subjected to cis-element analysis. The sequences were analyzed using PLACE (Plant Cis-element Regulatory DNA Elements, http://www.dna.affrc.go.jp/PLACE/) and plantCARE (plant cis-acting regulatory element, http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Verification of Differential Gene Expression by Quantitative Reverse-Transcription PCR
The transcript levels of genes expressed in different tissues or under different treatments were quantified by quantitative reverse-transcription PCR (qRT-PCR) using the same samples chosen for transcriptome analysis. The experiments were performed in an Applied Biosystems 7500 real-time PCR system using a SYBR Premix ExTaq Kit (TaKaRa Code: DRR041A, Japan). The reaction system and procedure used were described previously [63]. Data were processed using the 2 −∆∆CT method. The primers are listed in Table S4.

Conclusions
In this study, we examined the physiological and transcriptional response of industrial rapeseed to drought and salinity treatment. Short-term abiotic treatments caused a significant increase in proline content in the seedlings, while the levels of antioxidant enzymes remained unchanged. GO and KEGG analysis indicated that most of the identified DEGs were involved in the metabolic process, response to stimulus, or biosynthesis of secondary metabolites. Some stress-responsive DEGs were shared between plants subjected to drought and salinity treatments.