Transcriptome Profiling of Haloxylon persicum (Bunge ex Boiss and Buhse) an Endangered Plant Species under PEG-Induced Drought Stress

Haloxylon persicum is an endangered western Asiatic desert plant species, which survives under extreme environmental conditions. In this study, we focused on transcriptome analysis of H. persicum to understand the molecular mechanisms associated with drought tolerance. Two different periods of polyethylene glycol (PEG)-induced drought stress (48 h and 72 h) were imposed on H. persicum under in vitro conditions, which resulted in 18 million reads, subsequently assembled by de novo method with more than 8000 transcripts in each treatment. The N50 values were 1437, 1467, and 1524 for the control sample, 48 h samples, and 72 h samples, respectively. The gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis resulted in enrichment of mitogen-activated protein kinase (MAPK) and plant hormone signal transduction pathways under PEG-induced drought conditions. The differential gene expression analysis (DGEs) revealed significant changes in the expression pattern between the control and the treated samples. The KEGG analysis resulted in mapping transcripts with 138 different pathways reported in plants. The differential expression of drought-responsive transcription factors depicts the possible signaling cascades involved in drought tolerance. The present study provides greater insight into the fundamental transcriptome reprogramming of desert plants under drought.


Introduction
H. persicum (Bunge ex Boiss and Buhse), the white saxaul belongs to the family Amaranthaceae, is an extremely drought-tolerant tree species with sand dunes as the habitat. It grows up to a height of 4 m. In the United Arab Emirates (UAE), H. persicum is generally found in gravel plains and sand dunes and has a great potential in landscaping. The tree is evergreen in the habitat with the leaves retrogressed as succulent branches of the tree. The plant has an extensive root system, making it useful for stabilizing sequencing (NGS), efficient technologies are available for identifying the gene regulatory mechanism in plants. Deep ribonucleic acid (RNA) sequencing studies are achievable for all plant species by exploring transcriptome profiles of model and non-model organisms, under various conditions. These high-throughput sequencing techniques are widely used to create transcriptome data, which enables researchers to identify the differential expression profile of genes in plants under different treatment conditions [17][18][19]. In recent times, transcriptomic approaches have gained much attention for identifying the gene regulatory network associated with drought-stress tolerance in plants [20,21]. In this study, polyethylene glycol (PEG) 6000 (10%) was supplemented in culture medium to induce drought stress in H. persicum under in vitro conditions. Understanding the transcriptome profile under PEG-induced drought stress will enable us to identify the differential expression profile of genes under stress. Although some differences may occur between the responses to PEG-induced drought and the drought under field conditions, the present study will evidently provide greater insight into the regulatory networks associated with drought tolerance in these xerophytic plant species that can adapt to harsh climates of the arid region.
Genes 2020, 11, x FOR PEER REVIEW 3 of 19 exploring transcriptome profiles of model and non-model organisms, under various conditions. These high-throughput sequencing techniques are widely used to create transcriptome data, which enables researchers to identify the differential expression profile of genes in plants under different treatment conditions [17,18,19]. In recent times, transcriptomic approaches have gained much attention for identifying the gene regulatory network associated with drought-stress tolerance in plants [20,21]. In this study, polyethylene glycol (PEG) 6000 (10%) was supplemented in culture medium to induce drought stress in H. persicum under in vitro conditions. Understanding the transcriptome profile under PEG-induced drought stress will enable us to identify the differential expression profile of genes under stress. Although some differences may occur between the responses to PEG-induced drought and the drought under field conditions, the present study will evidently provide greater insight into the regulatory networks associated with drought tolerance in these xerophytic plant species that can adapt to harsh climates of the arid region

Plant Material and Drought Treatment
Two-year-old seedlings of H. persicum maintained in the greenhouse conditions, used for the study were provided by Environmental Agency, Abu Dhabi (UAE). Apical shoot bud explants collected from the seedlings were surface sterilized by washing initially with running tap water and then with 1% Tween 20 solution to remove any adhering contaminants. Under sterile conditions, the explants were treated with 20% sodium hypochlorite for 10 min and washed 4-5 times with sterile distilled water. The sterilized explants were inoculated onto a Murashige and Skoog (MS) medium and incubated at 25 ± 2 °C under 16 h photoperiod (3678 Lux) and 8 h dark conditions.
The one-month-old adventitious shoot buds of H. persicum cultures were treated with PEG 6000 10% (−0.60 MPa) in MS media for drought induction under in vitro conditions. The shoot buds were maintained at 16 h photoperiod and 8 h dark conditions at 25 ± 2 °C. Samples were collected at 48 h and 72 h of treatment under PEG media. The MS medium without PEG (−3.2 MPa) served as the control. Each treatment had three replicates and the samples were collected after the preferred incubation period. The water potential of the drought-treated and control samples were measured using a WP4C water potential meter (Meter, Pullman, WA, USA). The relative water content (RWC) was measured using the following equation:

Plant Material and Drought Treatment
Two-year-old seedlings of H. persicum maintained in the greenhouse conditions, used for the study were provided by Environmental Agency, Abu Dhabi (UAE). Apical shoot bud explants collected from the seedlings were surface sterilized by washing initially with running tap water and then with 1% Tween 20 solution to remove any adhering contaminants. Under sterile conditions, the explants were treated with 20% sodium hypochlorite for 10 min and washed 4-5 times with sterile distilled water. The sterilized explants were inoculated onto a Murashige and Skoog (MS) medium and incubated at 25 ± 2 • C under 16 h photoperiod (3678 Lux) and 8 h dark conditions.
The one-month-old adventitious shoot buds of H. persicum cultures were treated with PEG 6000 10% (−0.60 MPa) in MS media for drought induction under in vitro conditions. The shoot buds were maintained at 16 h photoperiod and 8 h dark conditions at 25 ± 2 • C. Samples were collected at 48 h and 72 h of treatment under PEG media. The MS medium without PEG (−3.2 MPa) served as the control. Each treatment had three replicates and the samples were collected after the preferred incubation period. The water potential of the drought-treated and control samples were measured using a WP4C water potential meter (Meter, Pullman, WA, USA). The relative water content (RWC) was measured using the following equation [22]: where, FW is fresh weight, DW is dry weight and TW is turgid weight. DW is calculated after drying the samples for 24 h at 80 • C.

RNA Isolation and Transcriptome Library Preparation
The samples were collected from the PEG treatment at different periods of drought induction for RNA isolation and transcriptome profiling. The purification of RNA was carried out using a RNeasy mini kit (Qiagen, Germantown, MD, USA). The RNA concentration and purity were estimated using Nanodrop spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The RNA integrity of the samples was checked using a Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Briefly, 200 ng of RNA was used for mRNA isolation, fragmentation and priming. The first strand cDNA was synthesized in the presence of Actinomycin D, followed by second strand synthesis. The double stranded cDNA was purified using HighPrep PCR magnetic beads (Magbio Genomics Inc., Gaithersburg, MD, USA). Three replicates of each samples were pooled to form a single RNA sequencing library. The library was prepared with Illumina-compatible NEBNext ® UltraTM Directional RNA Library Prep Kit (New England BioLabs, Ipswich, MA, USA) based on the manufacturer's instructions. Transcriptome sequencing was carried out using Illumina HiSEQ2500 at Genotypic Technology (Bangalore, India). The purified RNA was stored at −80 • C for real-time fluorescent quantitative PCR validation.

Illumina Raw Data Processing and Transcriptome Assembly
The raw illumina paired end (PE) read quality was assessed using FastQC tool [23]. Reads were trimmed by removing 3' adapter, trimming ambiguous bases towards the read end and removal of PE reads, which had more than 30% of low quality (<q30) bases. Reads over 50 bp in length were retained after trimming process. Transcriptome de novo assembly was carried out using Trinity software, with the default k-mer value of 25 [24]. A Linux-based server with 256 GB of RAM and 20 cores was used for the transcriptome assembly. The redundant transcripts were clustered using CD-HIT program [25] to reduce redundancy and to obtain unigenes with sequence similarity (90%) and identity (95%) between the assembled transcripts. The clustered transcripts were used for further annotation with a minimum length of 300 bp.

Differential Gene Expression (DGE) Analysis
The fragments per kilobase per million fragments (FPKM) approach was used to calculate the gene expression. Transcripts from all the samples (control, T1 and T2) were pooled and clustered based on the similarity between transcripts using CD-HIT. The reads from different samples were aligned to the clustered transcriptome using the bowtie2 program [30], a pair of reads which supported each transcript were counted, and the significant gene expression difference was calculated using the DESeq tool [31] (FDR less than 0.01 were considered as a significant result). A heat map for differential expression was generated using the R program (https://www.r-project.org/).

qPCR Validation
The RNA-sequence differential gene expression was confirmed using qPCR. Ten genes that showed the significant differential expression during the DGE analysis were used for qPCR confirmation. The total RNA was converted into cDNA using an Affinity Script qPCR cDNA synthesis kit (Agilent Technologies, Santa Clara, CA, USA) as per the manufacturer's protocol. In brief, 250 ng of RNA from each sample was taken for cDNA synthesis and the first strand cDNA was synthesized using universal oligo dT primers. The cDNA was diluted to 20 ng/µL and 1 µL of was used for each qPCR reaction. Specific qPCR primers were designed using the Primer3plus program (Supplementary Table S1). The expression levels of selected genes were analyzed using SYBR Green chemistry (Brilliant II SYBR Green qPCR master mix (Agilent Technologies, Santa Clara, CA, USA)) in Stratagene mx3005P instrument (Agilent Technologies, Santa Clara, CA, USA). The amplification cycling conditions were as follows: initial denaturation for 95 • C for 10 min, followed by 40 cycles of 95 • C for 30 s, 60 • C for 30 s. The dissociation curve analysis was performed after amplification for primer specificity; the conditions were as follows: 95 • C for 1 min, 55 • C for 30 s and 0.2 • C/s increment up to 95 • C (continuously collect fluorescence from 55-95 • C). Three replicates were included for each sample and the mean Ct value of the technical replicates was used to calculate the expression level. The relative quantification method (2ˆ-∆∆Ct ) was used to calculate the fold change using GAPDH as an internal control. The statistical significance between the treated samples were conducted by Student's t-test using SPSS 16.0.

Physiological Changes in H. persicum
In the present study, in vitro drought stress studies in H. persicum was conducted by using 10% PEG 6000 (Figure 2a). Physiological parameters, such as fresh weight, dry weight, tissue water potential, and relative water content, were measured in all the treatment samples. The fresh weight and dry weight of the tissues showed significant changes under PEG treatment. There was a significant decline in fresh weight and dry weight observed when the shoot buds were subjected to PEG stress for 48 and 72 h. The fresh weight of the sample under control treatment was 0.505 g, which reduced significantly to 0.455 and 0.430 g in samples that are cultured on PEG medium for 48 h and 72 h, respectively. The dry weight of the tissues also showed gradual reduction upon increasing the PEG treatment period. The dry weight was 0.05 g under control condition, which reduced significantly to 0.0426 and 0.0353 g in tissues treated with PEG for 48 h and 72 h respectively (Figure 2b).
Relative water content and the tissue water potential (MPa) were monitored to confirm the influence of PEG treatment in imposing drought under in vitro condition. The relative water content of the samples showed significant changes under PEG treatment. While a maximum of 90% moisture content was observed in the control samples, under PEG treatments the relative water content was significantly reduced to 84% after 48 h. Due to the continued incubation of the cultures in the PEG treatment, the water content declined significantly to 81% after 72 h (Figure 2c). PEG 6000 treatment induced a significant reduction in tissue water potential. While the control samples without PEG

Data Processing and Transcriptome Assembly
The total RNA extracted from the control, T1, and T2 samples were converted into cDNA and sequenced using Illumina HiSeq2500 (151 × 2 bp chemistry). The transcriptome sequencing of the control, T1, and T2 libraries generated 19, 22.4, and 21.7 million reads ( Table 1). The sequenced raw reads were deposited in the NCBI-SRA database (control SAMN08281247, T1 sample SAMN08281244, and T2 sample SAMN08281248). Read trimming resulted in 15.2, 18 and 21.7 million high quality reads (Q > 30) for the control, T1, and T2 libraries, respectively. The length distribution of assembled transcriptome is shown in Figure 3 for the control, T1, and T2 sample, respectively. The adenine (A), thymine (T), guanine (G) and cytosine (C) content of each sample are listed in the Table  1 and Supplementary Figures S1a,b,c. The initial assembly of the control samples resulted in 110,448 contigs; these contigs were further clustered into 87,522 transcripts with a N50 value of 1437 bp and GC percentage of 38.71, which represents 85.6 MB of the control sample transcriptome. De novo assembly of both the T1 and T2 libraries resulted in 105,017 and 123,307 contigs; the clustering of these contigs generated 82,440 and 97,825 transcripts. The assembly for the T1 sample generated 81.3 MB of transcriptome with the GC content of 38.92% and N50 value of 1467 bp. Similarly, the T2 sample generated 98.9 MB of transcriptome with a N50 value and GC content of 1524 bp and 38.54%, respectively (Table 1).

Data Processing and Transcriptome Assembly
The total RNA extracted from the control, T1, and T2 samples were converted into cDNA and sequenced using Illumina HiSeq2500 (151 × 2 bp chemistry). The transcriptome sequencing of the control, T1, and T2 libraries generated 19, 22.4, and 21.7 million reads ( Table 1). The sequenced raw reads were deposited in the NCBI-SRA database (control SAMN08281247, T1 sample SAMN08281244, and T2 sample SAMN08281248). Read trimming resulted in 15.2, 18 and 21.7 million high quality reads (Q > 30) for the control, T1, and T2 libraries, respectively. The length distribution of assembled transcriptome is shown in Figure 3 for the control, T1, and T2 sample, respectively. The adenine (A), thymine (T), guanine (G) and cytosine (C) content of each sample are listed in the Table 1 and Supplementary Figure S1a-c. The initial assembly of the control samples resulted in 110,448 contigs; these contigs were further clustered into 87,522 transcripts with a N50 value of 1437 bp and GC percentage of 38.71, which represents 85.6 MB of the control sample transcriptome. De novo assembly of both the T1 and T2 libraries resulted in 105,017 and 123,307 contigs; the clustering of these contigs generated 82,440 and 97,825 transcripts. The assembly for the T1 sample generated 81.3 MB of transcriptome with the GC content of 38.92% and N50 value of 1467 bp. Similarly, the T2 sample generated 98.9 MB of transcriptome with a N50 value and GC content of 1524 bp and 38.54%, respectively (Table 1).

KEGG Analysis for Pathway Identification
The KEGG database was used to identify 138 important pathways involved in various developmental processes in plants. Some of the vital pathways, such as plant-hormone signal transduction pathways, plant pathogen interaction, MAPK signaling, phenylpropanoid biosynthesis, and ubiquitin-mediated proteolysis, were mapped with a significant number of expressed transcripts (Figure 7). From the assembled transcriptome data, we could find all the expressed genes mainly involved in the MAPK pathways. In total, 38 key genes involved in the MAPK pathway were identified (Supplementary Figure S2). The number of transcripts involved in the MAPK pathway varied from 83 to 94 transcripts in different sample groups. From the expressed transcriptome, 41 enzymes were found to be involved in plant-hormone signal transduction pathways. A significant number of transcripts were mapped against each plant-hormone-mediated signal cascade process (Supplementary Figure S3).

KEGG Analysis for Pathway Identification
The KEGG database was used to identify 138 important pathways involved in various developmental processes in plants. Some of the vital pathways, such as plant-hormone signal transduction pathways, plant pathogen interaction, MAPK signaling, phenylpropanoid biosynthesis, and ubiquitin-mediated proteolysis, were mapped with a significant number of expressed transcripts (Figure 7). From the assembled transcriptome data, we could find all the expressed genes mainly involved in the MAPK pathways. In total, 38 key genes involved in the MAPK pathway were identified (Supplementary Figure S2). The number of transcripts involved in the MAPK pathway varied from 83 to 94 transcripts in different sample groups. From the expressed transcriptome, 41 enzymes were found to be involved in plant-hormone signal transduction pathways. A significant number of transcripts were mapped against each plant-hormone-mediated signal cascade process (Supplementary Figure S3).

KEGG Analysis for Pathway Identification
The KEGG database was used to identify 138 important pathways involved in various developmental processes in plants. Some of the vital pathways, such as plant-hormone signal transduction pathways, plant pathogen interaction, MAPK signaling, phenylpropanoid biosynthesis, and ubiquitin-mediated proteolysis, were mapped with a significant number of expressed transcripts (Figure 7). From the assembled transcriptome data, we could find all the expressed genes mainly involved in the MAPK pathways. In total, 38 key genes involved in the MAPK pathway were identified (Supplementary Figure S2). The number of transcripts involved in the MAPK pathway varied from 83 to 94 transcripts in different sample groups. From the expressed transcriptome, 41 enzymes were found to be involved in plant-hormone signal transduction pathways. A significant number of transcripts were mapped against each plant-hormone-mediated signal cascade process (Supplementary Figure S3).

In Silico Differential Gene Expression (DGE) Analysis
Transcripts from the control, T1, and T2 samples were pooled and clustered using CD-HIT program, which resulted in 229,224 clustered transcripts. These clustered transcripts were considered as a reference for in silico DGE analysis. The read counts for all three samples were generated by aligning reads to the clustered transcripts using bowtie2, and normalized gene expression quantification was estimated using the DESeq program. Differentially expressed metabolically important enzymes were identified using KEGG database comparison (Figure 8). The DGE analysis between the control and T1 condition resulted in 1803 upregulated and 2501 downregulated transcripts (p-value < 0.01, Log2 fold change ±1). In the T2 samples, 3046 upregulated and 1582 downregulated transcripts were identified (p-value < 0.01, Log2 fold change ±1) (Figure 9). Approximately, 1153 transcripts showed significant differential expression in both T1 and T2 samples. Functional annotation for differentially expressed transcripts was performed by aligning against the Uniprot plant protein database. The differential expression of the transcription factors expressed under different treatments were represented as heat maps ( Figure 10). Compared to other TFs, the transcripts of the NAC transcription factor showed higher expression under different treatment levels.
Genes 2020, 11, x FOR PEER REVIEW 11 of 19 Transcripts from the control, T1, and T2 samples were pooled and clustered using CD-HIT program, which resulted in 229,224 clustered transcripts. These clustered transcripts were considered as a reference for in silico DGE analysis. The read counts for all three samples were generated by aligning reads to the clustered transcripts using bowtie2, and normalized gene expression quantification was estimated using the DESeq program. Differentially expressed metabolically important enzymes were identified using KEGG database comparison (Figure 8). The DGE analysis between the control and T1 condition resulted in 1803 upregulated and 2501 downregulated transcripts (p-value < 0.01, Log2 fold change ±1). In the T2 samples, 3046 upregulated and 1582 downregulated transcripts were identified (p-value < 0.01, Log2 fold change ±1) (Figure 9). Approximately, 1153 transcripts showed significant differential expression in both T1 and T2 samples. Functional annotation for differentially expressed transcripts was performed by aligning against the Uniprot plant protein database. The differential expression of the transcription factors expressed under different treatments were represented as heat maps (Figure 10). Compared to other TFs, the transcripts of the NAC transcription factor showed higher expression under different treatment levels.

RT-PCR
Quantitative real-time PCR analysis for the selected genes was performed to confirm the in silico DGE analysis data. From the transcript annotation, ten genes significantly expressed in both T1 (48 h PEG treatment) and T2 (72 h PEG treatment) compared to the control were selected for qPCR analysis. The analysis of ten genes, zinc finger (ZNF), cytochrome P450 (CYT), peroxidase (POX), glutamate decarboxylase (GAD), aldoketoreductase (AKR), calcium-transporting ATPase (CTAPSE), phosphoinositide phospholipase C (PI-PLCs), trehalase (TRE), protein phosphatase-2C (PP2C), and MYB family protein (MYB), showed positive correlations with the in silico DGE analysis data. The fold change trends in both the RNA sequence and the qPCR data are illustrated in Figure 11.

RT-PCR
Quantitative real-time PCR analysis for the selected genes was performed to confirm the in silico DGE analysis data. From the transcript annotation, ten genes significantly expressed in both T1 (48 h PEG treatment) and T2 (72 h PEG treatment) compared to the control were selected for qPCR analysis. The analysis of ten genes, zinc finger (ZNF), cytochrome P450 (CYT), peroxidase (POX), glutamate decarboxylase (GAD), aldoketoreductase (AKR), calcium-transporting ATPase (CTAPSE), phosphoinositide phospholipase C (PI-PLCs), trehalase (TRE), protein phosphatase-2C (PP2C), and MYB family protein (MYB), showed positive correlations with the in silico DGE analysis data. The fold change trends in both the RNA sequence and the qPCR data are illustrated in Figure 11.
Genes 2020, 11, x FOR PEER REVIEW 13 of 19 Figure 11. Validation of selected transcript from RNA seq using qPCR. The graph on the top shows the Log2FC of ten selected transcript DGE expressions in the T1 and T2 sample. The bottom graph represents the corresponding qPCR fold change in the T1 and T2 samples. The error bar represents the standard deviation. The difference between the means were calculated by Student's t-test (*: p ≤ 0.05, **: p ≤ 0.01). T1: 48 h PEG treatment; T2: 72 h PEG treatment.

Discussion
In arid regions drought, salinity and thermal stress are the major hostile environmental factors that often limit crop production. Climate change and rapid urbanization are the potential factors that cause the erosion of biodiversity leading to the extinction of many species. H. persicum, an endangered nearly extinct species native to the arid region of the Middle East, has a great potential in urban landscaping considering its ability to combat thermal and drought stress. Therefore, identifying and conserving the drought-responsive genes from such endangered species has considerable significance for enriching the gene pool and future applications in drought studies. Drought stress and oxidative stress are often interrelated and may also exhibit identical and complexly adverse effects, leading to cellular damage. Exploring the expression profile of drought-responsive genes is essential for understanding the molecular mechanisms and metabolic pathways associated with drought-stress tolerance in crops [6,32].
The focus of the present study was to identify the upregulated and downregulated genes involved in the growth and maintenance of adventitious shoot buds developed in vitro under PEGinduced drought. PEG 6000 (100 g L −1 ) was used in the media to create a rapid drought stress through water deprivation [33]. A significant decline in fresh weight, dry weight, RWC, and tissue water potential (Mpa) along with changes in the phenology of the adventitious shoot buds confirmed the PEG-induced drought stress. Transcriptional reprogramming under stress is one of the major mechanisms that plants undergo during stress tolerance. In the present study, transcriptome profiles of the two drought-induced (PEG) samples T1 (48 h) and T2 (72 h) were compared with the control Figure 11. Validation of selected transcript from RNA seq using qPCR. The graph on the top shows the Log2FC of ten selected transcript DGE expressions in the T1 and T2 sample. The bottom graph represents the corresponding qPCR fold change in the T1 and T2 samples. The error bar represents the standard deviation. The difference between the means were calculated by Student's t-test (*: p ≤ 0.05, **: p ≤ 0.01). T1: 48 h PEG treatment; T2: 72 h PEG treatment.

Discussion
In arid regions drought, salinity and thermal stress are the major hostile environmental factors that often limit crop production. Climate change and rapid urbanization are the potential factors that cause the erosion of biodiversity leading to the extinction of many species. H. persicum, an endangered nearly extinct species native to the arid region of the Middle East, has a great potential in urban landscaping considering its ability to combat thermal and drought stress. Therefore, identifying and conserving the drought-responsive genes from such endangered species has considerable significance for enriching the gene pool and future applications in drought studies. Drought stress and oxidative stress are often interrelated and may also exhibit identical and complexly adverse effects, leading to cellular damage. Exploring the expression profile of drought-responsive genes is essential for understanding the molecular mechanisms and metabolic pathways associated with drought-stress tolerance in crops [6,32].
The focus of the present study was to identify the upregulated and downregulated genes involved in the growth and maintenance of adventitious shoot buds developed in vitro under PEG-induced drought. PEG 6000 (100 g L −1 ) was used in the media to create a rapid drought stress through water deprivation [33]. A significant decline in fresh weight, dry weight, RWC, and tissue water potential (Mpa) along with changes in the phenology of the adventitious shoot buds confirmed the PEG-induced drought stress. Transcriptional reprogramming under stress is one of the major mechanisms that plants undergo during stress tolerance. In the present study, transcriptome profiles of the two drought-induced (PEG) samples T1 (48 h) and T2 (72 h) were compared with the control samples without any drought treatment. This resulted in 82,440 assembled transcript (T1) with an average length of 978 bp (control). The N50 lengths were 1437, 1467, and 1524 bp for the control, T1, and T2, respectively. These findings are comparable with the results of RNA-sequence in Haloxylon ammodendron (N50: 1345 bp) under drought stress [34].
During drought stress, significant changes in gene expression occur in plants in order to survive under drought by regulating cellular process [35]. Elucidating the drought-stress tolerance in plants by differential gene-expression methods provides valuable information for identifying probable stress tolerance mechanism. The DGE analysis displayed the differential expression of approximately 1153 transcripts in both T1 and T2 samples, which was relatively higher in number compared to H. ammodendron transcriptome [34]. In the present study, the number of DEGs showed a significant change when the PEG stress was increased from 48 h to 72 h. Similar findings were reported in cassava, where a significant changes observed in DEGs when the PEG stress period was increased from 3 h to 24 h [36]. In addition, our results indicate that the differential expression of downregulated genes showed gradual reduction in number upon increase in the PEG stress period, whereas the upregulated genes showed a steady increase in number under elevated stress period. Overall, the percentage of downregulated genes was higher in both treatment conditions ( Figure 9). These results were in alignment with the findings of Liu et al. (2014), and Ma et al. (2015), who reported that the abundance of downregulated DGEs was significantly higher compared to upregulated genes under PEG stress [16,17].
The transcription factors are regulatory proteins which play a major role in drought tolerance by synchronizing the signaling network and gene expression under stress [37]. NAC transcription factors are ubiquitous in plants, and their expression shows significant changes under abiotic stresses [38]. In our dataset, the expression of the NAC transcription factor (DN33705_c0_g1_i1) gene was upregulated under PEG stress, compared to the control ( Figure 10). The expression level elevated significantly when the PEG treatment period shifted from 48 h to 72 h. The involvement of the NAC TFs in drought tolerance has been demonstrated in various crops like rice [39], and soybean [40]. Dehydration responsive element binding factor (DREB) is an AP2/ERF-type transcription factor present in plants [41]. The overexpression of the OsDREB2B gene in Arabidopsis resulted in improved tolerance to drought stress [42]. In the present study, DREB (DN33605_c0_g1_i2) transcripts displayed significant upregulation under PEG-induced stress; this denotes the probable role of DREB TFs in the downstream process of the drought-tolerance mechanism in H. persicum. Other drought-responsive TFs, such as MYB, WRKY, bZIP, bHLH, and C2H2, showed differential expression under PEG-induced drought stress, as reported in other plant species [43].
The GO enrichment analysis was performed to predict the role of transcripts in the biological, cellular, and molecular process of the cells. In all the three sample groups, more than 50% of the DEGs were involved in various GO terms. Drought-induced genes are functionally classified as two major groups, namely functional proteins and regulatory proteins. Functional proteins are involved in stress tolerance and the regulatory proteins have a major role in stress response through alteration gene expression and signal transduction [7]. The KEGG pathway analysis resulted in mapping the DEGs with 138 different pathways that were involved in various cellular processes. The MAPK cascades are a multigene family of regulatory networks deployed for the transduction of intra-and extra-cellular signals to the nucleus for suitable adjustments of the cell to the stimuli [13]. In this study, 38 potential genes were mapped with the MAPK pathway, which includes the expression of all the three components of MAPK cascade (MAPKKK, MAPKK and MAPK). MAP kinases that are involved in ethylene, jasmonic acid, and abscisic acid-mediated stress responses are significantly mapped in the present study (Supplementary Figure S2).
Phytohormones play a significant role in plant stress regulation. During drought induction, the environmental stress factors induce the secretion of plant hormones which mediates the immediate cell responses by triggering plant hormone signal transduction pathways. Studies show that abscisic acid exhibit an important role in the drought tolerance of plants through stomatal closure and the activation of the stress-responsive genes [7]. In the present study, under PEG-induced stress, the ABA-mediated drought response facilitated by PP2C gene (Transcript Id: T2_DN28928_c0_g1_i1) was mapped, which in turn triggered the ABF (ABA-responsive element binding factor) (DN22654_c0_g1_i1) gene, which plays a major role in stomatal closure ( Figure 12). Similar responses were reported in Phormium tenax [44]. Reports suggests that, the overexpression of abscisic acid stress ripening (ASR) gene family play a significant role in drought and salinity stress [45]. We found upregulation of ASR gene (DN23961_c0_g2_i3) was observed under PEG induced drought stress, which indicates the possible role of ASR genes in ABA-mediated stress tolerance. Increased stress tolerance in response to the overexpression of ASR genes was reported in transgenic plants under various stress conditions [46][47][48][49].
Genes 2020, 11, x FOR PEER REVIEW 15 of 19 was mapped, which in turn triggered the ABF (ABA-responsive element binding factor) (DN22654_c0_g1_i1) gene, which plays a major role in stomatal closure ( Figure 12). Similar responses were reported in Phormium tenax [44]. Reports suggests that, the overexpression of abscisic acid stress ripening (ASR) gene family play a significant role in drought and salinity stress [45]. We found upregulation of ASR gene (DN23961_c0_g2_i3) was observed under PEG induced drought stress, which indicates the possible role of ASR genes in ABA-mediated stress tolerance. Increased stress tolerance in response to the overexpression of ASR genes was reported in transgenic plants under various stress conditions [46][47][48][49]. In the present study, the ethylene-mediated plant-hormone signal transduction pathway was also explored. Ethylene is an important plant-growth regulator which regulates the plant growth under drought condition by undertaking various developmental changes in the plant [50]. Mapping the DEG for ethylene-responsive transcription factor (ERF) 1 in the KEGG pathway denotes the possible activation of ethylene-induced defense mechanism under PEG-induced stress conditions. The role of ERFs in the developmental process such as cuticle biosynthesis was reported earlier in Arabidopsis [51]. The crosslink between the jasmonic acid and ethylene-mediated pathway observed in the MAPK-signaling cascade, shows the possible role of jasmonic acid signaling intermediates in ethylene-induced defense response (Supplementary Figure S2). These findings are in accordance with the observation of Lorenzo et al. (2002), who reported that the ERF1 expression could be regulated by ethylene and jasmonic acid or the synergistic effect of both hormones [52].
The accumulation of reactive oxygen species (ROS) is one of the major modifications to occur under drought or salt stress, causing damages to the cell membranes [53]. Antioxidant enzymes play a crucial role in drought-defense mechanisms by scavenging ROS generated under drought stress [54]. Peroxidase is an important antioxidant enzyme involved in the detoxification of H2O2 molecules generated under stress. In our study, 30 DEGs matching the expression of peroxidase (EC 1.11.1.7) enzyme were identified. The expression of the enzyme was greater in the T2 sample (72 h) compared to the T1 (48 h) treatment. This shows the increased accumulation of peroxidase enzyme with a concomitant increase in PEG-induced drought stress. A similar pattern of expression changes in In the present study, the ethylene-mediated plant-hormone signal transduction pathway was also explored. Ethylene is an important plant-growth regulator which regulates the plant growth under drought condition by undertaking various developmental changes in the plant [50]. Mapping the DEG for ethylene-responsive transcription factor (ERF) 1 in the KEGG pathway denotes the possible activation of ethylene-induced defense mechanism under PEG-induced stress conditions. The role of ERFs in the developmental process such as cuticle biosynthesis was reported earlier in Arabidopsis [51]. The crosslink between the jasmonic acid and ethylene-mediated pathway observed in the MAPK-signaling cascade, shows the possible role of jasmonic acid signaling intermediates in ethylene-induced defense response (Supplementary Figure S2). These findings are in accordance with the observation of Lorenzo et al. (2002), who reported that the ERF1 expression could be regulated by ethylene and jasmonic acid or the synergistic effect of both hormones [52].
The accumulation of reactive oxygen species (ROS) is one of the major modifications to occur under drought or salt stress, causing damages to the cell membranes [53]. Antioxidant enzymes play a crucial role in drought-defense mechanisms by scavenging ROS generated under drought stress [54]. Peroxidase is an important antioxidant enzyme involved in the detoxification of H 2 O 2 molecules generated under stress. In our study, 30 DEGs matching the expression of peroxidase (EC 1.11.1.7) enzyme were identified. The expression of the enzyme was greater in the T2 sample (72 h) compared to the T1 (48 h) treatment. This shows the increased accumulation of peroxidase enzyme with a concomitant increase in PEG-induced drought stress. A similar pattern of expression changes in genes associated with peroxidases have been reported in many different plant species under drought stress [55]. Several other antioxidant enzymes, such as superoxide dismutase (SOD), catalase (CAT) and glutathione s-transferase (GST), displayed an upregulated DEG profile, which confirms the biochemical rearrangements in H. persicum under the PEG-induced drought stress condition.

Conclusions
In this study, we investigated the transcriptome profile of H. persicum, an endangered desert plant species. This is the first report on this plant, and we reported~50,000 annotated transcripts and their possible metabolic pathways under PEG-induced drought stress. In this analysis, we identified the expression pattern of drought-responsive genes using differential gene expression analysis, which resulted in a minimum of 4304 DEGs under 48 h PEG treatment and a maximum of 4523 DEGs under 72 h PEG treatment. Elucidating the pathway results revealed the possible defense mechanism under PEG-induced drought in H. persicum. The present study will provide a better understanding of drought-response mechanism of the desert plant under stress conditions. Some differences may occur between the responses to PEG-induced drought and the drought under field conditions; nevertheless this provides valuable evidence to confirm the transcriptional reprogramming of genes under stress.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/6/640/s1, Supplementary Table S1: List of primers used for qPCR validation; Table S2: Transcriptome annotation file for control sample; Table S3: Transcriptome annotation file for T1 sample (48h of PEG stress); Table S4: Transcriptome annotation file for T2 sample (72h of PEG stress); Figure  Funding: This study was supported by grants from Khalifa Center for Genetic Engineering and Biotechnology (31R004), UAE University.