Transcriptomic analysis of Eruca vesicaria subs. sativa lines with contrasting tolerance to polyethylene glycol-simulated drought stress

Background Eruca vesicaria subsp. sativa is one of the Cruciferae species most tolerant to drought stress. In our previous study some extremely drought-tolerant/sensitive Eruca lines were obtained. However little is known about the mechanism for drought tolerance in Eruca. Methods In this study two E. vesicaria subs. sativa lines with contrasting drought tolerance were treated with liquid MS/PEG solution. Total RNA was isolated from 7-day old whole seedlings and then applied to Illumina sequencing platform for high-throughput transcriptional sequencing. Results KEGG pathway analysis indicated that differentially expressed genes (DEGs) involved in alpha-Linolenic acid metabolism, Tyrosine metabolism, Phenylalanine, Tyrosine and tryptophan biosynthesis, Galactose metabolism, Isoquinoline alkaloid biosynthesis, Tropane, Piperidine and pyridine alkaloid biosynthesis, Mineral absorption, were all up-regulated specifically in drought-tolerant (DT) Eruca line under drought stress, while DEGs involved in ribosome, ribosome biogenesis, Pyrimidine metabolism, RNA degradation, Glyoxylate and dicarboxylate metabolism, Aminoacyl-tRNA biosynthesis, Citrate cycle, Methane metabolism, Carbon fixation in photosynthetic organisms, were all down-regulated. 51 DEGs were found to be most significantly up-regulated (log2 ratio ≥ 8) specifically in the DT line under PEG treatment, including those for ethylene-responsive transcription factors, WRKY and bHLH transcription factors, calmodulin-binding transcription activator, cysteine-rich receptor-like protein kinase, mitogen-activated protein kinase kinase, WD repeat-containing protein, OPDA reductase, allene oxide cyclase, aquaporin, O-acyltransferase WSD1, C-5 sterol desaturase, sugar transporter ERD6-like 12, trehalose-phosphate phosphatase and galactinol synthase 4. Eight of these 51 DEGs wre enriched in 8 COG and 17 KEGG pathways. Conclusions DEGs that were found to be most significantly up-regulated specifically in the DT line under PEG treatment, up-regulation of DEGs involved in Arginine and proline metabolism, alpha-linolenic acid metabolism and down-regulation of carbon fixation and protein synthesis might be critical for the drought tolerance in Eruca. These results will be valuable for revealing mechanism of drought tolerance in Eruca and also for genetic engineering to improve drought tolerance in crops. Electronic supplementary material The online version of this article (10.1186/s12870-019-1997-2) contains supplementary material, which is available to authorized users.


Background
Drought is one of the most important environmental stresses in the world [1][2][3]. Previous studies indicated that plant responses to drought stress is quite complicated and exploring drought tolerance mechanism in plants is still a challenge [4][5][6][7]. Next generation sequencing (NGS) using Illumina HiSeq2000 covers the transcriptome many folds and allows quantification of the detected transcripts [8,9]. By using next generation sequencing the gene expression is quantified, making it possible for quantitative comparisons [10]. Many researchers have used transcriptome sequencing to study the responses to drought stress in crops [10][11][12][13][14][15][16].
Biological and genetic diversity exists among and within plant species regarding drought tolerance, however many of these adaptive mechanisms are not completely understood. Genome-wide identification of drought-responsive transcriptions using genotypes with contrasting droughttolerance will help revealing interactions among metabolic pathways in response to drought stress [12]. Eruca vesicaria subsp. sativa is one of the Cruciferae species most tolerant to drought stress [17,18]. In our previous study some extremely drought-tolerant/sensitive Eruca lines were obtained [19]. In this study, Illumina sequencing technology was used to identify the differentially expressed genes (DEGs) and their biochemical pathways related to drought tolerance in drought tolerant (DT) versus droughtsusceptible (DS) Eruca genotypes. Expression level of some selected drought-responsive genes was validated by using quantitative RT-PCR. To our knowledge, this is the first report about the whole transcriptome analysis in E. vesicaria subsp. sativa. The results provide transcriptome data valuable for understanding drought tolerance mechanism in Eruca and the differentially expressed genes could be valuable for improving drought tolerance in crops.

Results
Overview of Illumina RNA sequencing and de novo assembly of E. vesicaria subs. sativa transcriptome In total 314,322,210 raw reads after sequencing produced 254,660,002 clean reads including 430,166 unique genes with N50 value of 638 bp and average length of 511.56 bp. The unigene length ranged from 201 to 16, 517 bp ( Table 1).

GO anotation and enrichment
GO (Gene ontology) with Nr annotation indicated that 4953 unigenes were classified into 69 groups that could be categorized into three main classifications: "biological process" (2859), "cellular component" (643) and "molecular function" (1451). In the "molecular function" class the largest number of unigenes (984) were involved in "catalytic activity"; in the "cellular component" class the largest number of unigenes (168) were involved in "cell part"; in the "biological process" class the largest number of unigenes (534) were involved in "metabolic process" (Additional file 1).  Referring to Audic's algorithm [20], FDR (False Discovery Rate) was calculated based on the P-value that corresponds to differential expression tests of transcripts. By using "FDR ≤ 0.05 and the absolute value of log 2 Ratio ≥1" as the threshold, we identified in the DS-MS vs DS-PEG group 3287 unigenes up-regulated and 6318 downregulated; in the DT-MS vs DT-PEG couple, 2299 unigenes up-regulated and 2355 down-regulated (Fig. 4). It was specifically in DT-MS vs DT-PEG but not in DS-MS vs DS-PEG that 2946 DEGs were detected, among which 1584 were up-regulated and 1362 down-regulated, 449 had Nr annotation and 552 annotated to 232 KEGG pathways (data not shown).

Discussion
Comparative study using lines with contrasting droughttolerance is a useful tool for identifying drought-responsive genes. Identification of DEGs exclusively in the tolerant genotype would be valuable for revealing the mechanisms responsible for stress tolerance [21]. In the present study we provided detailed transcriptomic profiles of whole seedlings of two Eruca lines with contrasting drought tolerance.

ABA responses and stress signaling
Drought stress induces ABA accumulation, which in turn leads to stomatal closure to keep the water status in plants [22][23][24]. Evidence indicates that WRKY proteins, including those induced by ABA, are upregulated under drought stress in rice [25]. OsWRKY45 overexpression increased drought tolerance [26], while activated expression of WRKY57 resulted in drought tolerance in Arabidopsis [27]. In this study 19 of 26 DEGs for WRKY transcription factor were significantly up-regulated specifically in the drought tolerant Eruca line under drought stress (hereafter referred as significantly upregulated), among which one for WRKY transcription factor 55-like (c195662 g1 i2) was most significantly upregulated specifically in drought-tolerant Eruca line PI 426649 under PEG-simulated drought stress (hereafter referred as most significantly up-regulated).
Basic helix-loop-helix (bHLH) genes are important in phytohormone signaling. Rice OsbHLH148 confers drought tolerance by interacting with OsJAZ proteins [28]. Arabidopsis bHLH122 and PebHLH35 from Populus euphratica are positive regulators of drought and salt tolerance and osmotic signaling [29,30]. In this study 5 DEGs for transcription factor bHLH were down-regulated and 8 up-regulated, among which one for transcription factor bHLH19-like (c196938 g3 i6) was most significantly up-regulated (Additional file 6).

Antioxidants and ROS modulation
Reactive oxygen species (ROS) often acts as secondary messenger modulating stomatal closure in responses to different stimuli [22,37]. In Arabidopsis, calmodulinbinding transcription activator 1 (AtCAMTA1) is involved in regulation of membrane integrity by inducing ABA responses to drought stress. The camta1 mutants are highly susceptible to drought stress [38]. In this study 10 DEGs for calmodulin-binding protein were all up-regulated, among which one for calmodulin-binding transcription activator 2-like (c195596 g1 i7) was most significantly up-regulated (Additional file 6).
Mos Cysteine-rich receptor-like kinases (CRKs) are regulated by ROS [39]. It was found that CRK45 positively regulated plant responses to drought and salt stresses by inducing expression of ABA-responsive and stress-inducible genes. CRK45 overexpression enhanced drought tolerance in plants [40]. Overexpression of Arabidopsis CRK5 also increased ABA sensitivity and thus enhanced drought tolerance [41]. In this study 10 of 12 DEGs for cysteine-rich receptor-like protein kinase were up-regulated, among which one for cysteine-rich receptor-like protein kinase 17 (c196393 g1 i2) was most significantly up-regulated (Additional file 6).
It has been shown that MAPKs are involved in plant signal transduction in responses to different environmental stimuli [42][43][44][45][46]. MAPK, MAPKK and MAPKK kinase constitute a functional MAPK cascade. Activation of MAPK helps its translocation to nucleus to phosphorylate and activate transcription factors [47]. AtMKK1 in Arabidopsis activates AtMPK3 to transfer abiotic stess signals [48]. NtMEK2 (MAPKK) activates WIPK and SIPK for drought signal transmission [49]. Overexpression of GhMKK1 increased drought tolerance in Nicotiana benthamiana [41]. Overexpression of AtMKK1 in Arabidopsis decreased ROS levels and increased drought tolerance, while AtMKK1 deficiency resulted in elevated ROS and increased sensitivity to drought tolerance [50]. Plants over expressing AtMKK4 accumulated fewer ROS and showed less water-loss under drought stress [51]). In this study one DEG for mitogen-activated protein kinase kinase 4 (c198487 g1 i1) was most significantly up-regulated (Additional file 6).
NAC genes such as SNAC3 contributes to drought resistance and osmotic modulation independent of ABA [52]. SNAC3 could interact with WD domain-containing protein to adjuste ROS in rice. In this study all 13 DEGs for NAC domain-containing protein were up-regulated and one DEG for WD repeat-containing protein 53 (c198319 g1 i2) was most significantly up-regulated (Additional file 6).

Alpha-linolenic acid, jasmonate signaling and cell membrane stability
Plants keep membrane fluidity and integrity by modulating oleic and linolenic acid levels during stresses [53]. Alphalinolenic acid is a major precursor for messengers including jasmonic acid generated by oxidative modifications [54][55][56]. Antisense expression of Arabidopsis omega-3 fatty acid desaturase gene led to reduced salt/drought tolerance in tobacco [57], while over-expression of FAD3 or FAD8 resulted in increased drought tolerance [58]. Lenka et al. [12] found that eight enzymes involved in alpha-linolenic acid metabolism were significantly induced by drought stress in drought-tolerant rice. In this study we also found all DEGs involved in alpha-linolenic acid metabolism up-regulated (Additional file 5).
Jasmonic acid (JA) plays critical role in stomatal closure during drought stress [59]. It is also involved in the production of antioxidants regulating ascorbate and glutathione metabolism [60]. Based on the pathway established by Vick and Zimmerman [61], JA is produced from alpha-linolenic acid, and in this process allene oxide cyclase and OPDA reductase (OPR) play important roles [62][63][64][65][66]. In this study all four DEGs of OPDA reductase and 6 DEGs for allene oxide cyclase were up-regulated, among which two DEGs for allene oxide cyclase (c184414 g1 i1, c184414 g2 i1) were most significantly up-regulated (Additional file 6).
Cuticular waxes play important roles in reducing nonstomatal water loss under stresses [69]. C-5 sterol desaturase catalyzes the incorporation of C-5 double bond into D7-sterols to form D5, 7-sterols. Overexpression of fungal C-5 sterol desaturase increased wax deposition and drought tolerance in tomato [70]. In this study one DEG for C-5 sterol desaturase (c201142 g1 i2) was most significantly up-regulated (Additional file 6).
O-acyltransferase WSD1 is also involved in cuticular wax biosynthesis [71]. Zhang et al. [35] indicated that KCS and WSD, and their up-stream regulators, were upregulated under drought stress. In this study one DEG for O-acyltransferase WSD1-like (c175262 g2 i1) was most significantly up-regulated (Additional file 6).

Trehalose, galactinol and sugar transport
In plants various abiotic stresses lead to sugar accumulation in the vacuole [72][73][74], suggesting that sugar biosynthesis and vacuolar sugar transporters play important roles under these conditions. Sugar transporter ERD6 expression is induced by drought stress [75], while ERD6-like transporter (ESL1) expression is enhanced by drought and ABA treatment [76]. In this study 2 DEGs for sugar transporter ERD 6-like 12 (c197650_g1_i3, c197650_g1_i5) were most significantly up-regulated (Additional file 6).
Trehalose is also important in protecting plants from stresses [77]. Trehalose-6-phosphate (T6P) is suggested to act as signaling metabolite responsing to the environment [78]. In plants T6P is dephosphorylated by trehalose-6phosphate phosphatase to produce trehalose [79]. In this study all 6 DEGs for trehalose-phosphate phosphatase were up-regulated, among which 3 were most-upregulated (Additional file 6).

Down-regulation of carbon fixation and protein synthesis
Lenka et al. [12] found that under drought stress carbon fixation was up-regulated, while in this study we found that in the drought-tolerant Eruca line DEGs for carbon fixation in photosynthetic organisms were all downregulated under drought stress (Additional file 5). Ribosomes are places where proteins are synthesized. Dhindsa and Bewley [86] found decline of amino acid incorporation and polysome population in the droughttolerant moss Tortula ruralis under drought stress. In this study KEGG pathway analysis indicated that DEGs involved in ribosome, ribosome biogenesis, Aminoacyl-tRNA biosynthesis were all down-regulated in drought tolerant Eruca under PEG-simulated drought stress (Additional file 5); Blast analysis also indicated that 9 of 10 DEGs for translation initiation factor, 8 DEGs of ATP synthase, 6 DEGs for ribosomal RNA gene, 5 of 6 DEGs for eukaryotic translation initiation factor, 118 of 121 DEGs for ribosomal protein, were down regulated, among which 38 DEGs for ribosomal protein were most significantly down-regulated (Additional file 6), suggesting that Eruca might also respond to drought stress by down-regulation of protein synthesis and carbon fixation in photosynthetic organisms.
KEGG pathways in which all DEGs were up-regulated, such as Tyrosine metabolism, Phenylalanine, tyrosine and tryptophan biosynthesis, Galactose metabolism, Isoquinoline alkaloid biosynthesis, Tropane, piperidine and pyridine alkaloid biosynthesis, Mineral absorption; KEGG pathways in which most DEGs were up-regulated, such as Porphyrin and chlorophyll metabolism, Ubiquinone and other terpenoid-quinone biosynthesis, Arachidonic acid metabolism, Glutathione metabolism, glycerophospholipid metabolism, Phenylalanine metabolism, Plant-pathogen interaction, were up-regulated (Additional file 5). Blast analyses indicated that all 5 DEGs for dehydration-responsive elementbinding protein, 9 of 10 DEGs for glutathione S-transferase, 8 of 10 DEGs for transcription factor MYB, 4 of 5 DEGs for calcineurin B-like protein, all 16 DEGs for receptor-like protein kinase, 15 of 17 DEGs for zinc finger CCCH domain-containing protein, all 14 DEGs for F-box protein, were up-regulated (Additional file 6). These up-regulated pathways and DEGs might also play important roles in drought tolerance in Eruca.

Conclusions
Based on the transcriptomic analyses we postulated that the most significantly up-regulated DEGs for Ethyleneresponsive transcription factors, WRKY and bHLH transcription factors involved in stress signaling and ABA responses; DEGs for calmodulin-binding transcription activator, cysteine-rich receptor-like protein kinase, mitogen-activated protein kinase kinase and WD repeatcontaining protein involved in antioxidants and ROS modulation; DEGs for OPDA reductase and allene oxide cyclase involved in JA production, C-5 sterol desaturase involved in producing D5, 7-sterols, aquaporin involved in decreasing ion leakage and malondialdehyde, Oacyltransferase WSD1 involved in cuticular wax biosynthesis; trehalose-phosphate phosphatase, galactinol synthase 4 and sugar transporter ERD 6-like 12 involved in osmoprotectant production, up-regulation of DEGs involved in alpha-linolenic acid metabolism, and down-regulation of carbon fixation and protein synthesis might be critical for the drought tolerance in Eruca line PI 426649. These results might be valuable for revealing mechanism of drought tolerance in Eruca and for improving drought tolerance in crops.

Materials
In this work Eruca vesicaria subsp. sativa PI 426649 highly tolerant and PI 426652 highly sensitive to PEGsimulated drought stress [19] were used for transcriptomic analysis.

Tissue sampling, RNA extraction, library preparation and Illumina sequencing
Eruca seeds of PI 426649 and PI 426652 originally from the Agricultural Research Service, USDA were germinated on filter paper immersed in liquid MS medium [87] without sugar or organic components. Seven days later the seedlings were treated for 10 h with 20% PEG-6000/liquid MS and then harvested and frozen imediately in liquid nitrogen and then stored at − 80°C. Drought-tolerant PI 426649 is denoted as 'DT and Drought-sensitive PI 426652 is denoted as 'DS'. Treatment with liquid MS medium is denoted as 'MS' while treatment with PEG is denoted as 'PEG'. Four samples (DT-MS, DT-PEG, DS-MS, DS-PEG) each with two biological replicates were taken for the present study. Total RNA was isolated from the whole seedlings that had been stored at − 80°C by using TRIZOL total RNA extraction reagent (TAKARA) according to the manufacturer's protocol. RNA integrity was verified by 1.5% Agrose gel electrophoresis and confirmed using a 2100 Bioanalyzer analyzer (Agilent, CA, USA). The mRNA enrichment, RNA fragmentation, the first and second strand cDNA synthesis and purfying, sequencing adaptors ligation and PCR amplification were performed as described [14]. The libraries were applied to Illumina sequencing platform (HiSeq 2000, SanDiego, CA, USA) for high-throughput sequencing using a paired-end read protocol with 100 bp of data collected per run.

De novo sequence assembly, clustering and homology search
After sequencing, the raw image data was transformed into sequence data by base calling using CASAVA package provided by Illumina and saved as raw reads of fastq format and then treated with Trimmomatic(v0.30) to get clean reads. Transcriptome de novo assembly was carried out with short reads assembling program-Trinity (r2013-02-25) [88]. Unigenes from each sample's assembly were taken into further process of sequence splicing and redundancy removing to acquire non-redundant unigenes as long as possible by TGICL [89].

Functional annotation
Blastx alignment (E-value < 0.00001) between unigenes and protein databases like Nr (NCBI non-redundant database), KEGG (Kyoto Encyclopedia Of Genes and Genomes, www.genome.jp/kegg/) and COG (Clusters of Orthologous Groups of proteins, http://www.ncbi.nlm. nih.gov/COG/) is performed as described [90]. After Nr annotation, Blast2GO program was used to get GO annotation for the unigenes. WEGO software was used to do GO functional classification for all unigenes and to understand the species distribution of gene functions at the macro level [91]. KEGG pathway annotation was carried out by using bidirectional best hit method (BBH) at KEGG Automatic Annotation Server (KAAS).

Analysis of differentially expressed genes (DEGs), GO and KEGG pathway enrichment
Reads of the drought tolerant and sensitive Eruca accessions with PEG/MS treatment were mapped back to our de novo assembling results using RSEM [92]. To evaluate the gene expression, the number of unique match reads was calculated and then normalized to FPKM (Fragments per Kilo base of transcript per Million mapped reads) and then used to calculate the unigene expression with restrictive conditions of | log 2 Ratio | ≥1.0 and FDR ≤ 0.05. GO enrichment analysis of these DEGs was performed using blast2GO with P-value ≤1 and pathway enrichment analysis was carried out using Path finder software against the KEGG database with Q-value ≤1.
Validation of the differentially expressed genes by qRT-PCR