Analysis of the Citrullus colocynthis Transcriptome during Water Deficit Stress

Citrullus colocynthis is a very drought tolerant species, closely related to watermelon (C. lanatus var. lanatus), an economically important cucurbit crop. Drought is a threat to plant growth and development, and the discovery of drought inducible genes with various functions is of great importance. We used high throughput mRNA Illumina sequencing technology and bioinformatic strategies to analyze the C. colocynthis leaf transcriptome under drought treatment. Leaf samples at four different time points (0, 24, 36, or 48 hours of withholding water) were used for RNA extraction and Illumina sequencing. qRT-PCR of several drought responsive genes was performed to confirm the accuracy of RNA sequencing. Leaf transcriptome analysis provided the first glimpse of the drought responsive transcriptome of this unique cucurbit species. A total of 5038 full-length cDNAs were detected, with 2545 genes showing significant changes during drought stress. Principle component analysis indicated that drought was the major contributing factor regulating transcriptome changes. Up regulation of many transcription factors, stress signaling factors, detoxification genes, and genes involved in phytohormone signaling and citrulline metabolism occurred under the water deficit conditions. The C. colocynthis transcriptome data highlight the activation of a large set of drought related genes in this species, thus providing a valuable resource for future functional analysis of candidate genes in defense of drought stress.


Introduction
Water is essential for plant growth in modern agriculture [1]. Drought delays the development of crops, and strongly affects morphology, as well as physiological processes like transpiration, photosynthesis, respiration and translocation of assimilates [2]. Drought avoidance can be achieved through morphological changes in plants, such as decreased stomatal conductance, reduced leaf area, and extensive root systems [3]. Drought tolerance is achieved by physiological and molecular mechanisms, including osmotic adjustment, antioxidant and scavenger compounds [4]. Both strategies involve the induction of specific genes and proteins, such as dehydrins (dehydration-induced proteins), key enzymes for osmolyte biosynthesis, and detoxification enzymes [5,6].
Plant species have developed diverse strategies to adapt and thrive in all kinds of climates and terrains to deal with extreme changes in the environment. These strategies are supported by rich and complex metabolic networks that enable the plant to synthesize a wide range of compounds. Plant responses to abiotic stresses involve interactions and crosstalk between many molecular pathways. High throughput screening techniques such as transcriptome sequencing have been used to study the adaptability of plants to drought [7]. This led to the discovery of many drought related genes. For example, PIP aquaporins were found to fine-tune the environment in response to declining water availability [8]. However, few natural allelic variants have been cloned for drought related traits, so QTL, RNA sequencing, and other trait isolation methods are needed to improve methodology for exploring complex multivariate data [9,10].
The cucurbit family is a large family with several economically important species, such as watermelon (Citrullus lanatus), melon (Cucumis melo), cucumber (Cucumis sativus) and several Cucurbita species with edible fruits [11]. Citrullus colocynthis (L.) Schrad (2n = 2x = 22), the bitter apple, closely related to domesticated watermelon (Citrullus lanatus var. lanatus), is a very droughttolerant perennial herbaceous species in the Cucurbitaceae family [12]. It can survive arid environments by maintaining its water content under severe stress conditions. C. colocynthis is an important medicinal plant and a source of valuable oil [12]. Its seeds were found in several early Egyptian, Libyan and Near Eastern sites from about 4000 BC. This species grows in sandy areas throughout northern Africa, southwestern Asia and the Mediterranean region [11,12]. The species has been used as a model to elucidate the function of genes implicated in the stress response ultimately leading to enhancement of stress tolerance in cucurbit crops through genetic manipulation. Si et al. [13] found dynamic gene expression changes in C. colocynthis root tissues using cDNA amplified fragment length polymorphism (cDNA-AFLP) technique.
Several research groups have used next generation sequencing technologies to study gene expression profiles in species of the cucurbit family. Guo et al [14] used 454 sequencing technology to study the comprehensive profile for watermelon fruit flesh tissues, Grassi et al [15] studied carotenoid pathway regulators in ripening watermelon fruit. The draft genome of watermelon (C. lanatus, 2n = 2x = 22, ,425 Mb) was analyzed by Guo et al [16] using three different watermelon subspecies. Comparative genomic analysis provided an evolutionary scenario for the origin of the 11 watermelon chromosomes derived from a 7-chromosome paleohexaploid eudicot ancestor. The genome sequence of cucumber (Cucumis sativus, 2n = 2x = 14) has been completed, and the genome of melon (Cucumis melo, 2n = 2x = 24) is being sequenced under the Spanish Genomics Initiative (MELONO-MICS) [17,18]. Liu et al [19] used sequencing techniques to identify conserved and novel miRNA in watermelon, while Wincker [20] used comparative analysis of genomes between watermelon and sweet orange to detect the traits related to their domestication.
Here, high-throughput sequencing of the leaf transcriptome from C. colocynthis provides a glimpse at drought related genes in this uniquely drought tolerant cucurbit species. The lack of extensive genomic and functional genomic resources in Citrullus species has hampered research and breeding of the cultivated and economically important C. lanatus. This study should facilitate the identification of valuable multiple genes, needed for complex interactions of plant species with the environment.

Plant Materials and RNA Extraction
C. colocynthis seedlings were grown in Sunshine Mix #8 under a 16 h light/8 h dark photoperiod at 26uC day, 22uC night temperature. Seedlings with 2-3 true leaves (2-3 week old) in 50 ml containers were exposed to drought by withholding water for 0 (day 1, D1), 24 (day 2, D2), 36 (day 3, D3) or 48 hours (day 4, D4). True leaf samples were collected each day at noon, flash frozen in liquid nitrogen and stored at 280uC. RNA was subsequently extracted using the TRIzol method [21].

Preparation of cDNA Library and Sequencing
Illumina sequencing was performed at the HudsonAlpha Institute of Biotechnology (Huntsville, AL) following manufacturer's instructions. RNA-Seq reads were first processed to remove rRNA sequence contamination. First strand cDNA was synthesized with reverse transcriptase and random primers using the small fragment RNAs as template. Second strand cDNA was then synthesized followed by phosphorylation by T4 DNA polymerase. The cDNA fragments were 39 adenylated and ligated to the Illumina's paired-end adapters. Enrichment of cDNA templates was conducted following fifteen cycles of PCR amplification. In total, over 20*4Mp were sequenced for mapping assembly and differential expression analysis. Raw sequence data are available for download at NCBI Sequence Read Archive (SAMN02769576, SAMN02769577, SAMN02769578, and SAMN02769579).

Assembly
Raw sequencing data were filtered Z (0.05) using the CLC Genomics workbench (CLC Bio, Aarhus, Denmark). Paired-end sequences from the four samples were used to construct the de novo C. colocynthis transcriptome assembly with default parameters. Assembly reads were also assembled against the watermelon Table 1. Paired end read statistics before and after trimming of the Citrullus colocynthis transcriptome at 4 time points. (C. lanatus) genome sequences, which were downloaded from the Cucurbit Database (http://www.icugi.org/cgi-bin/ICuGI/index. cgi). Reads were filtered and assembled using the CLC genomics workbench. The parameters used were as follows: 2 points of mismatch cost, 2 points of insertion cost, 2 points of deletion cost, 0.5 as length fraction, 0.95 as similarity fraction. After the de novo assembly and watermelon mapping assembly, we used Trinity to assemble all the contigs with the default parameters.

Identification of Full Length cDNAs
Two methods were used to identify full length cDNAs. First, BLASTX searching (E value: 1e 210 ) was used to detect the matched cDNAs in SwissProt database; second, ESTScan 2.0 was used to identify the translated sequences. Sequences with either start codon (ATG) and stop codon (TAG/TGA/TAA), or sequences with start codon (ATG) and homologue to a known protein with $80% similarity, were chosen as full length cDNAs.

Expression Analysis Using Transcriptome Reference
Pair-end sequencing reads of the four libraries were filtered using CLC Genomics workbench (0.05) before mapping to the references sequences from assembled cDNAs. First, read counts of each unigene were converted to reads per million (RPM).   Secondly, statistical analysis using Kal's test [22] provided in CLC Genomics workbench (P,0.05 and fold change $1.5) was conducted. These transcripts were annotated against the reference sequences.

Gene Ontology Analysis
The functional annotation software BLAST2GO (http://www. blast2go.com/b2ghome) was used to conduct gene ontology (GO) analysis of C. colocynthis genes in this study. The databases used were SwissProt and NCBI. BLAST E-value was set at 1.0e 23 . The major GO analysis was determined by BLAST, mapping, and annotation. Results are presented as a bar chart showing the percent of genes belonging to each group.

qRT-PCR Analysis
For cDNA synthesis, 500 ng of the total RNA for each sample (the same RNA was used for RNA-seq analysis) was used in reverse transcription with ProtoScript First Strand cDNA synthesis kit from New England BioLabs (Ipswich, MA). qRT-PCR was performed with SYBR-Green Supermix from Bio-Rad (Hercules, CA) in an Eppendorf Mastercycler ep realplex (Hauppauge, NY). Table S1 contains gene specific primer sequences. Each reaction contained 10 ml of SYBR-Green supermix, 1 ml of cDNA template, 1 ml forward primer (4 mm), 1 ml reverse primer (4 mm), and 7 ml ddH 2 O. The qRT-PCR program consisted of one cycle at 95uC for 15 sec, followed 40 cycles of 15 sec at 95uC, 15 sec at 55uC, and 30 sec at 72uC. The relative expression data was compared with actin [13] from C. colocynthis. Quantification  Table 3. Validation of the RNA-Seq expression profiles of selected C. colocynthis genes by qRT-PCR.   of the relative transcript levels was performed using the comparative C T method. The induction ratio (IR) was calculated as recommended by the manufacturer and corresponds to 2 2DDCT , where DDCT = (C T, target gene , -C T, actin ) treatment -(C T, target -C T, actin ) control . All experiments were replicated three times.

Transcriptome Assembly Results
The transcriptome of C. colocynthis leaves following 4 days of drought stress was assembled and assessed following paired-end (2*50 bp) Illumina sequencing. The Illumina platform yielded an average of 24 million high-quality reads per sample (Table 1). All sample reads were used to construct a de novo assembly. A reference assembly was constructed using the completely sequenced watermelon genome. A total of 20,581 contigs were generated ( Table 2). The contigs had an average length of 1350 bp and N50 of 1870 bp. BLASTX was used against SwissProt and ESTscan of translated protein sequences to detect full length cDNAs. A total of 5,038 full length cDNAs were detected in our sequencing assembly (Table 2).
Principle component analysis (PCA) was implemented in the CLC workbench. The results (Figure 1) illustrate differential gene expression patterns in C. colocynthis seedlings following four days of withholding water. Gene expression patterns on D1 and D2 were nearly similar, while results from D3 and D4 showed substantial differences as compared to D1 and D2, and gene expression patterns detected on D4 were very different from all other time points as a result of drought stress.

Differentially Expressed Genes Involved in Response to Drought Stress in C. colocynthis
The paired-end reads were mapped to the reference genomes after filtration. Read counts of each unigene were converted to reads per million (RPM). The read number of each cDNA was divided by the total number of reads per day (1, 2, 3, or 4) from the data set, and multiplied by 10 6 . Statistical analysis was conducted using Kal's test (p,0.05 and fold change $1.5). Genes showing non-significant and significant changes in read counts are shown in Table S2. Each sample was compared to the day 1 sample reads for analysis of their significance level of gene expression. The read results indicated that 59 genes showed significant changes at D2, D3, and D4 as compared to D1; 13 genes showed significant differential expression at D2 and D3 as compared to D1; 30 genes  Figure 2 corresponded to the principal component analysis. D3 and D4 transcripts were clustered together, and D1 and D2 transcripts were clustered. Also significant changes were seen at D3 and D4 suggesting that the transcriptional response of many genes was up-regulated during drought. Strong effects were especially observed on D3 and D4.

Gene Ontology (GO) Classification
To functionally categorize significantly changed genes in C. colocynthis under drought treatment, gene ontology analysis by BLAST2GO was performed. C. colocynthis unigenes were categorized in three main GO categorizes: biological process (2672), molecular function (1368) and cellular component (1053). These GO terms were further divided into several sub-categories ( Figure 3). In the biological process category, single organism process genes accounted for more than 20% of the biological process genes. In the molecular function category, more than 40% of genes were associated with a catalytic activity. In the cellular category, more than 35% of the genes were associated with the cellular component.

Validation of Illumina Expression Patterns by qRT-PCR Analysis
To confirm the reliability of the Illumina sequencing read analysis, 8 candidate genes were selected and their expression was compared at D4 and D1 using qRT-PCR. The expression patterns resulting from qRT-PCR showed general agreement with those from the Illumina sequencing analysis (Table 3). Discrepancies with respect to ratio of fold changes between sequencing and qRT-PCR analysis can be attributed to the essentially different algorithm and sensitivity of the two techniques [23]. In the deep-sequencing method the absolute expression rather than relative expression as in qRT-PCR analysis is used. Transcrip- tional qRT-PCR analyses of 16 genes during the drought treatments are shown in Figures 4 and 5. Gene Comp5873 is homeobox-leucine zipper protein, with significant upregulation during all days of drought in C. colocynthis. It is known that the expression of several homeoboxleucine zipper proteins is correlated to stress. Athb-12, a homeobox-leucine zipper domain protein from Arabidopsis, is functionally involved in salt tolerance in yeast [24]. Hahb-4, a homeobox-leucine zipper gene is potentially involved in water stress in sunflower [25]. Comp862 belongs to the glutathione Stransferase (GST) family, which contains heterogeneous, multifunctional dimeric proteins. This gene is highly up-regulated (60x) during drought in C. colocynthis. It is thought that GSTs are involved in cellular detoxification [26]. Comp1108 is a member of the NAC (NAM/ATAF1, 2/CUC2) gene family and NACs are known to be involved in numerous biological processes, including drought stress [27].
Comp10156 is a member of the GID1 (GIBBERELLIN INSENSITIVE DWARF1) family in C. colocynthis with high expression under drought conditions. GID1 is a soluble GA receptor in rice [28]. GA-GID1 complex interacts with DELLA proteins, which are negative regulators of GA signaling pathway. RD22, a known dehydration responsive gene in Arabidopsis, which is mediated by abscisic acid (ABA), may have physiological and molecular significance for processes underlying memory functions of plants in response to ABA and light pulses [29,30]. One RD22-like protein from soybean can alleviate salinity and osmotic stress [31]. RD22-like genes (Comp 6528) with significant up-regulation (.160x) under drought in C. colocynthis were confirmed by qRT-PCR, especially after days of withholding water.
NDR1 (NON-RACE-SPECIFIC DISEASE RESISTANCE1) plays an important role in coordinating broader cellular processes in response to stress and bacterial pathogen infection [32]. The chemical b-aminobutyric acid, which is known to induce resistance in plants, primed the expression of many genes, and NDR1/ NHL10 was one of them [33]. Comp7317 gene, which is a member of NDR1/NH10 showed significant changes at D2 only during the early stage of water deficit stress. GRAS (for GA Insensitive, REPRESSOR of gal-3 (RGA), SCARECROW (SCR)) transcription factors, have a major function in plant development and environmental adaption. These TFs are particularly implicated in the modulation of plant tolerance to stressors as cold, drought, salinity by crosstalks via GA to ABA-dependent and ABA-independent pathways [34]. For example, SCL7 confers salt and drought tolerance in Arabidopsis [35]. SCL14 is involved in the detoxification of xenobiotics and possibly endogenous harmful metabolites [36]. Comp20554, which belongs to the GRAS transcription factor family, showed significant up-regulation especially on D2.
The expression profile of Comp3048, which codes for a methyltransferase, maintained high levels during drought stress, especially on D2. It was found that myo-inositol-O-methyltransferase (Imt1) responded to low temperature stress in transgenic Arabidopsis [37]. Trithorax-like H3K4 methyltransferase from barley is drought inducible [38]. The methylation of myo-inositol catalyzed by myo-inositol methyltransferase (IMT) occurs when plants are under abiotic stress. Over-expression of IMT resulted in improved tolerance to dehydration and salt stress treatment in Arabidopsis [39].
Heat-shock proteins (HSPs) are environmentally induced proteins that enable plants to cope with heat and other environmental stresses. For example, Trichoderma harzianum Hsp70 transgenic Arabidopsis is abiotic stress tolerant [40]. Similarly HSP22 was found to be highly upregulated in C. colocynthis roots during drought conditions [13]. Overexpression of GmHsp90S can decrease damage of abiotic stresses in Arabidopsis [41]. Comp14675 belongs to the HSP family, and showed up-regulation at later stages of D3 and D4 ( Figure 5).
Plant cold shock proteins (CSP) are very conserved among various plant genera [42]. The first CSP identified was WCSP1 from winter wheat, which did accumulate in response to low-temperature stress [43]. Similarly in C. colocynthis CSP (Comp13927) was up-regulated during later stages of water deficit.
Comp2553, one MA3 domain-containing protein gene, showed marked changes on D3 and D4 under drought conditions ( Figure 5). Loss-of-function of ECIP1 (one MA3 domain-containing protein) resulted in enhanced ethylene response but altered salt response [44] in Arabidopsis.
Overexpression of plant BL-1 in Arabidopsis resulted in the attenuation of cell death induced by biotic stresses (pathogens) and abiotic stresses such as heat, cold, drought, salt and chemicalinduced oxidative stresses [45]. BL-1 might function to control the level of the ''pro-survival and pro-death'' signals under multiple stress conditions in plants [46]. For example, cucumber BAX inhibitor-1 is a conserved cell death suppressor induced by cold stress and a negative regulator of programmed cell death (PCD) [47]. Comp372 encodes one BL-1 gene, which was up-regulated by drought at later stages (D3 and D4).
Comp8117, which is a homologue of a Populus EST (CU233481.1), is a drought stress related gene, up-regulated at D4. Comp10586, which encodes a member of the MYB transcription factor family, is up-regulated during the late stage of drought stress. Some MYB members have been shown to regulate plant responses to biotic and abiotic stress conditions [48]. For example, AtMYB96 acts through the ABA signaling pathway to induce pathogen resistance by promoting salicylic acid biosynthesis, and thus regulating stomata movement, drought tolerance and disease resistance in Arabidopsis [49,50]. MYB88 might function directly or indirectly, as positive regulator of stressresponsive genes [51]. TaMYB30-B from wheat did improve drought stress tolerance in transgenic Arabidopsis [52]. Another family of transcription factors, the WRKY family, named for the WRKY domain of about 60 amino acids, contains a highly conserved WRKYGQK heptapeptide at its N-terminus and a zinc-finger-like motif at its C-terminus [53]. WRKY transcription factors are involved in multiple aspects of plant growth, development and stress [54]. Several TaWRKY in wheat with roles in the abiotic stress response acted in an ABA-dependent manner [55]. Here, gene Comp19751encodes a WRKY gene, which showed remarkable up-regulation at D4. MADS-box family members function in reproductive development and stress [56]. For example, OsMADS25 and OsMADS27 transcripts accumulate in response to osmotic stress [57]. Comp6823, which encodes for a MADS-box gene in C. colocynthis, is up-regulated markedly at the last stage of drought (D4; Figure 5).

Analysis of the Drought Stress Signaling Transcriptome in C. colocynthis
Drought stress signal transduction consists of several pathways including ionic and osmotic homeostasis signaling, detoxification response and growth regulation pathways [58]. Genes detected in C. colocynthis leaves during water deficit are listed in Table 4.
In the ion signaling pathway, calcium binding proteins are well known for their involvement in both biotic and abiotic stress response pathways [59]. The calcium ion (Ca 2+ ) as a secondary messenger in plants is sensed by calmodulins (CaMs)/CaM-like protein (CMLs), the caldineurin B-like proteins (CBLs) and Ca2+dependent protein kinases (CDPKs). CaM binds to CaM-binding proteins (CBPs), which function in different pathways under biotic and biotic stresses [60]. A total of 10 Ca 2+ binding proteins were detected in the C. colocynthis transcriptome.
Protein phosphorylation is a central theme in the cell's response to stress. The MAP kinase cascade in transcript levels consist of a number of protein kinases, such as two-component histidine kinase, MAPKKK, MAPKK, MAPK etc. [61]. Here we detected 16 MAP kinases in the C. colocynthis transcriptome.
Detoxification signaling can ameliorate the damage in plants under stresses [63], as noted in many other plant species. A total of 32 detoxification proteins were detected in the C. colocynthis transcriptome.
Molecular mechanisms regulating gene expression in response to drought stress have been studied by analyzing the functional transcription factors in ABA-dependent and ABA-independent pathways [6]. Several regulatory proteins in ABA-dependent orindependent pathway were detected, among which NAC, MYB/C and leucine-rich repeat proteins (LRR). These regulatory proteins can further modulate many responsive transcription factors. Several functional proteins such as heat shock protein (HSP) 70, HSP22, grpE like protein, respiratory burst oxidase D (RBOHD), VIRE2-Interacting protein2 (VIP2), ABA transporter-like protein, synaptobrevin-related protein, translocon outer envelope of chloroplast (Toc34-1), beta-amylase, TIP1 (TIP GROWTH DEFECTIVE 1), and RD22 were detected.

Analysis of Phytohormone Signaling Mediators in C. colocynthis
Phytohormones play important roles in regulating plant responses under biotic and abiotic stress. Elaborate phytohormone signaling networks mediate the adaptability of plants to different environmental conditions [64]. Many phytohormones such as ABA, salicylic acid (SA), jasmonic acid (JA), auxin, ethylene and gibberellic acid (GA) are being studied for their role in abiotic stress responses [33,65,66]. It is known that downstream signaling proteins for auxin, GA, JA, and ABA are subjected to ubiquitindependent degradation [67]. Putative phytohormone signaling genes detected in C. colocynthis during the drought response are listed in Table 5.
Ethylene signaling pathway components were ordered into a hypothetical linear pathway based on both genetic (epistasis) analysis and biochemical interactions [68]. Almost all of the ethylene signaling homologous members (118 ethylene response factors or ERFs) were detected in the C. colocynthis transcriptome. APETALA2/ethylene responsive factor (AP2/ERF) transcription factors are well-known for mediating stress responses and development in plants [69].
The auxin response factor (ARF) family contains transcription factors that bind to auxin-responsive elements (AREs) in the promoters of primary auxin-responsive genes. Aux/IAAs are early auxin-response proteins that bind ARFs, therefore inhibiting ARE-mediated gene transcription. Aux/IAAs are involved in ubiquitin-mediated degradation, which is catalyzed by SCF E3 ubiquitin ligase. TIR1 can stimulate Aux/IAA proteolysis by binding auxin to this protein [70]. All of the major auxin signaling related transcription factors found in the C. colocynthis transcriptome are shown in Table 5. Similar to the auxin pathway, a novel family of transcriptional regulators, the jasmonate ZIM-domain (JAZ) proteins play a target part as Aux/IAA. In addition, COI1 plays a similar role as TIR1, while MYC and R2R3-MYB transcription factors work as ARFs [71]. Several JA signaling pathway related genes exist in the C. colocynthis transcriptome.
In the SA pathway, NPR1 (non-specific disease resistance 1) is a key regulator in SA-dependent defense signaling [72]. Similarly WRKY and TGA play major roles as transcriptional regulators in the SA pathway. We detected 2 NPR1 related proteins, 46 WRKYs and 6 TGAs in C. colocynthis, which might function in the C. colocynthis SA pathway.
Drought triggers the production of the phytohormone ABA, which in turn causes stomatal closure and induces expression of stress-related genes. The soluble PYR/PYL/RCAR receptors function at the apex of a negative regulatory pathway to directly regulate PP2C phosphatases, which in turn directly regulate SnRK2 kinases [73]. Several of the core transcription factors in the ABA pathway are listed in Table 5.

Cellular Metabolism under Drought Stress in C. colocynthis
Global gene expression analyses have shown substantial downregulation of many photosynthetic genes under drought not only in Arabidopsis [76], but also several other species such as indica rice [77,78]. Similarly, many photosystem I and II, chlorophyll a, b binding protein, and oxygen evolving enhancer protein genes [79], showed down-regulation during water deficit stress in C. colocynthis (Table 6).
Citrulline, a non-protein amino acid intermediate in the arginine biosynthetic pathway, has been found to accumulate in leaves of drought tolerant watermelon under water deficit conditions [80,81,82]. Thus several factors related to the response of this model drought tolerant species to stress have been identified (Table 7, Figure S1). Citrulline metabolic genes (carbamoyl-phosphate synthetase, acetyl glutamate synthase, acetylornithine aminotransferase, aminoglutamate decarboxylase, acetylornithine deacetylase, and glutamate dehydrogenase) were found to be significantly up-regulated during drought.
One of the major research goals is to understand the molecular mechanisms underlying drought tolerance in plants. It is clear that drought triggers a wide variety of responses in C. colocynthis. Down regulation of many photosynthetic genes was observed especially at the later stages of drought. Up regulation of many transcription factors, stress signaling factors, detoxification genes, and genes involved in phytohormone signaling occurred throughout the water deficit experiment. Systematic approaches using genomic analyses should lead to the discovery of additional stress factors and provide us with a better understanding of the stress tolerance mechanism of this drought tolerant plant species. Figure S1 Citrulline metabolism pathway in C. colocynthis.

Author Contributions
Conceived and designed the experiments: ZW FD. Performed the experiments: ZW FD. Analyzed the data: ZW FD HH LG. Contributed reagents/materials/analysis tools: FD LG SM. Contributed to the writing of the manuscript: ZW FD.