Introduction {#s1}
============

Water deficit is one of the most significant abiotic stresses that influence germination, growth, development, and productivity of plants [@pone.0094277-Okcu1]. Drought stress causes stomatal closure, limited gas exchange, and reduced photosynthesis in plants [@pone.0094277-Lawson1]. Over-reduction of photosynthetic electron transport chain induces the generation of reactive oxygen species (ROS), such as singlet oxygen (^1^O~2~), superoxide anion (O~2~ ^−^), hydrogen peroxide (H~2~O~2~), and hydroxyl radical (•OH), which damage cellular structures and macromolecules [@pone.0094277-Dat1]. Plants have enzymatic and non-enzymatic systems to eliminate ROS. Moreover, plants have developed drought-resistance strategies such as succulent leaves, formation of osmophilic globules, stomatal movement, and reduction of leaf water potential [@pone.0094277-Liu1], [@pone.0094277-McDowell1]. The understanding of plant responses to water deficit have improved because of the application of molecular techniques. The expression of numerous genes in response to drought stress was described in previous studies [@pone.0094277-Fan1], [@pone.0094277-Hazen1]. Seki et al. [@pone.0094277-Seki1] classified the genes expressed during stress into two groups: (i) genes encoding proteins involved in signal transduction (protein kinases and transcription factors) and (ii) genes with products, such as late embryogenesis abundant (LEA) proteins, chaperone, osmoprotectants, and detoxification enzymes, that directly protect cells against stress. Genomic and transcriptomic analyses revealed that various transcriptional regulatory systems are involved in stress-responsive gene induction. Several different sets of *cis*- and *trans*-acting factors are known to be induced by drought stress at the molecular level [@pone.0094277-YamaguchiShinozaki1]. A large number of metabolites and proteins have also been reported to be up regulated in response to drought stress [@pone.0094277-Shinozaki1].

*Reaumuria soongorica* (Pall.) Maxim. is an extreme xerophytic semi-shrub and a typical constructive and dominant species in the desert vegetation in China. *R. soongorica* forms the zonal landscape and is widely distributed in northwest China [@pone.0094277-Ma1]. Water supply is one of the main limiting factors in the habitat of this species. The unique adaptive strategies in the morphology and physiology of *R. soongorica*, such as thick cuticle, hollow stomata, and accumulation of some low-molecular-weight metabolites, are attributed to its special distribution area [@pone.0094277-Liu1]. *R. soongorica* is a vascular flowering plant with desiccation-tolerance. *R. soongorica* leaves wither and enter a state of dormancy during dehydration but is revived when water becomes available. Thus, *R. soongorica* is utilized as a valuable non-model species in exploring drought-tolerance mechanisms. Liu et al. [@pone.0094277-Liu1] found sucrose, malate, and proline, which play important roles in osmoregulation in *R. soongorica* during water loss. Studies on the protection mechanism for photosynthetic properties, activity of antioxidant enzyme, and metabolite changes under drought stress have also been conducted [@pone.0094277-Liu2], [@pone.0094277-Bai1]. However, the molecular mechanism of drought tolerance in *R. soongorica* is still poorly understood. A preliminary genome survey, which included genome size, chromosome number, and karyotype, was conducted with *R. soongorica* [@pone.0094277-Wang1]. In our previous study, the mitogen-activated protein kinase gene (*RsMPK2*) was isolated and its participation as possible mediator under different stresses was investigated [@pone.0094277-Liu3]. We also identified and characterized the response of flavonone 3-hydrolase gene (*RsF3H*), a key enzyme-encoding gene involved in flavonoid biosynthesis pathway in drought treatment [@pone.0094277-Liu4]. Recently, the transcriptome of *R. soongorica* was sequenced and analyzed to obtain valuable genetic information for the investigation on the molecular mechanism of drought tolerance [@pone.0094277-Shi1]. Genetic information from the transcriptome of *R. soongorica* was also used to identify differentially expressed genes (DEGs) under drought stress to understand the mechanism of drought tolerance.

In recent years, with the increasing availability of sequence data, expression profiling has been used to identify DEGs in different tissues, organs, developmental stages or mutants and determine the expression patterns induced by stress in a large number of model and non-model organisms. Microarray analysis and suppression subtractive hybridization are powerful tools for isolating differentially expressed cDNAs that mediate complex biological processes in higher eukaryotes [@pone.0094277-Diatchenko1]. These technologies have been used in a number of animals and plants, such as *Aedesaegypti* [@pone.0094277-Lertkiatmongkol1], soybean [@pone.0094277-Asakura1], *Populuscanadensis* [@pone.0094277-Caruso1], *Scylla paramamosain* [@pone.0094277-Chen1], cucumber [@pone.0094277-Zhou1], and navel orange [@pone.0094277-Da1]. Next-generation sequencing technologies offer new approaches for global measurements of gene expression with the advantage of high efficiency and low cost. Digital gene expression (DGE) tag profiling is based on ultra-high-throughput sequencing of cDNA fragments that uniquely tag the corresponding gene, thereby allowing direct quantification of transcript abundance [@pone.0094277-Hong1]. DGE technology allows identification of millions of DEGs without the need for prior annotations and permits the analysis of organisms that lack genomic information. Thus, DGE has been widely utilized to monitor the differences in transcriptional responses among different tissues and organs in response to stress in silkworm [@pone.0094277-Gao1], *Sagittariatrifolia* [@pone.0094277-Cheng1], *Dugesia japonica* [@pone.0094277-Qin1], *Populus* [@pone.0094277-Chen2], cotton [@pone.0094277-Wei1], spruce [@pone.0094277-Albouyeh1], *Brassica napus* [@pone.0094277-Jiang1], and moss [@pone.0094277-Nishiyama1].

In this study, DEGs in leaves of *R. soongorica* in response to PEG-induced drought stress were examined using DGE tag profiling technology. Analysis of gene expression related to stress response provides further insight into the molecular mechanisms of stress tolerance in *R. soongorica*. Based on the putative functions of the identified genes, some important genes may be cloned. Moreover, the cloning of stress tolerance genes and determination of their expression patterns under drought stress may reveal attractive candidate genes and valuable information to improve drought stress tolerance of plants through genetic engineering.

Materials and Methods {#s2}
=====================

Ethics Statement {#s2a}
----------------

*R. soongorica* is widely distributed in northwest China and is not listed as endangered or protected species. No specific permissions are required for sample collection in the Northern Mountain of Lanzhou, and that the field studies did not involve any endangered or protected species.

Plant Materials and Drought Treatment {#s2b}
-------------------------------------

Seeds of *R. soongorica* were collected from the Northern Mountain of Lanzhou, Gansu, China (36°17′N, 103°48′E, 1700--1900 m elevation). Seeds were planted in the breeding base at the foothill. Three-year-old plants were transferred to a greenhouse under the following conditions: 150 µmol m^−2^s^−1^ (fluorescent tubes), photoperiod 16 h light/8 h darkness, day/night temperature of 25°C/20°C. The plants were washed with deionized water to remove surface contaminants. Soil attached to roots was also rinsed off with deionized water. The plants were then cultivated in 0.125×Murashige and Skoog medium solution (without vitamins and organics) with insufflating air for two weeks to accommodate the liquid environment during treatments. The samples were classified into untreated and PEG-treated groups. The untreated samples were maintained in cultivating conditions as previously described. For drought treatment, 15% PEG6000 was added to the culture solution of *R. soongorica*. All plants were treated for 12 h, from 6∶00 to 18∶00 local time. Leaves from each treatment were collected and immediately stored in liquid nitrogen for later use.

DGE Library Preparation and Sequencing {#s2c}
--------------------------------------

Total RNA samples were isolated according to modified CTAB-method [@pone.0094277-Wang2]. Quality and quantity analysis of total RNA, library construction, and sequencing were performed at the Beijing Genome Institute at Shenzhen, China. An extract of 8 µg of total RNA was obtained and treated with oligo (dT) magnetic bead adsorption to purify mRNA. Oligo (dT) was then used as primer in the synthesis of first- and second-strand cDNA. The 5′-ends of tags can be generated by two types of endonucleases, namely, NlaIII or DpnII. The bead-bound cDNA was subsequently digested with restriction enzyme NlaIII, which recognizes and severs CATG sites. Fragments other than 3′-cDNA fragments connected to oligo (dT) beads were washed away, and the Illumina adaptor 1 was ligated to the sticky 5′-end of the digested bead-bound cDNA fragments. The junction of the Illumina adaptor 1 and the CATG site constitutes the recognition site of MmeI, which is a type of endonuclease with separate recognition and digestion sites. MmeI breaks 17 bp downstream of the CATG site, producing tags with adaptor 1. After removing 3′ fragments with magnetic bead precipitation, the Illumina adaptor 2 was ligated to the 3′-ends of tags, thereby acquiring tags with different adaptors of both ends to form a tag library. After linear PCR amplification, fragments were purified by PAGE gel electrophoresis. During the quality control steps, Agilent 2100 Bioanalyzer and ABI Step One Plus Real-Time PCR System were used in the quantification and qualification of the sample library. Finally, the library was sequenced using Illumina HiSeq 2000 sequencer. The available raw data was then deposited in the NCBI Gene Expression Omnibus (GEO) database (<http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE55412>).

Analysis and Mapping of DGE Tags {#s2d}
--------------------------------

Sequencing-received raw image data were transformed by base calling into sequence data. Prior to mapping the readings from the reference database, all sequences were filtered to remove 3′ adaptor sequence, low-quality sequences (tags with unknown sequences 'N'), tags with a copy number of one (probably sequencing error), and tags that were long or too short, leaving clean tags of 21 nt. For annotation, all clean tags were mapped in reference sequences, and a mismatch of only 1 bp was considered. We obtained the combined transcriptome of *R. soongorica* with (suffix is "\_3″) and (suffix is "\_A") (<http://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-1543/>) [@pone.0094277-Shi1], which was used as reference database. Clean tags were mapped in reference sequences. The clean reads mapped in reference database from multiple genes were filtered. The remaining clean tags were used as unambiguous clean tags. The number of unambiguous clean tags for each gene was calculated and normalized to TPM (transcripts per million clean tags).

Identification and Functional Annotation of DEGs {#s2e}
------------------------------------------------

We used a rigorous algorithm method to identify DEGs between two samples. False discovery rate (FDR) was applied to determine the threshold of P-value in multiple tests and analysis. The DEGs were obtained through FDR ≤0.001 and \|log~2~Ratio\| ≥1. Gene ontology (GO) and pathway classification was used to determine the possible functions of all DEGs by mapping the GO (<http://www.geneontology.org/>) and KEGG (<http://www.genome.jp/kegg/>) databases.

Results {#s3}
=======

Analysis of DGE Libraries {#s3a}
-------------------------

The Illumina platform was used to perform high throughput tag-seq analysis on *R. soongorica* to investigate transcriptome responses to PEG-induced drought stress. Total RNA isolated from untreated and PEG-treated groups were named as libraries C and P, respectively. The major characteristics of the two libraries are summarized in [Table 1](#pone-0094277-t001){ref-type="table"}. An average of 7,231,417 total sequence tags per library was obtained with 466,957 distinct tag sequences from the two libraries. Prior to mapping tag sequences to the reference sequences, adaptor tags were filtered (low-quality tags and tags with one copy) and produced approximately 6,969,782 total clean sequence tags per library with 211,695 distinct clean tag sequences. The distributions of the total and distinct clean tag copy numbers showed highly similar tendencies in the two libraries ([Figure 1](#pone-0094277-g001){ref-type="fig"}). Among the distinct clean tags, more than 60% of the distinct clean tags had 2 to 5 copies, 34% of the distinct clean tags had 5 to 100 copies, and 4.5% had copy numbers higher than 100.

![Distribution of total clean tag and distinct clean tag copy numbers from the libraries C and P.\
C represents the control group; P represents the PEG-treated group.](pone.0094277.g001){#pone-0094277-g001}

10.1371/journal.pone.0094277.t001

###### Categorization and abundance of tags.

![](pone.0094277.t001){#pone-0094277-t001-1}

  Summary                                                              C         P
  --------------------------------- ------------------------------ --------- ---------
  Raw tags                                      Total               7266032   7196802
                                            Distinct tags           445588    488325
  Clean tags                                 Total number           7026443   6913120
                                         Distinct tag numbers       212338    211052
  All tag mapping to gene                    Total number           5869765   5356098
                                        Total % of clean tags       83.54%    77.48%
                                         Distinct tag numbers       110536    105096
                                     Distinct tag % of clean tags   52.06%     49.8%
  Unambiguous tag mapping to gene            Total number           4261924   3722859
                                        Total % of clean tags       60.66%    53.85%
                                         Distinct tag numbers        85870     81776
                                     Distinct tag % of clean tags   40.44%    38.75%
  All tag-mapped genes                          Number               39354     38413
                                           \% of ref genes          41.48%    40.49%
  Unambigous tag-mapped genes                   Number               27207     26381
                                           \% of ref genes          28.68%    27.81%
  Unknown tags                               Total number           1156678   1557022
                                        Total % of clean tags       16.46%    22.52%
                                         Distinct tag numbers       101802    105956
                                     Distinct tag % of clean tags   47.94%     50.2%

Clean tags are tags that remained after filtering out dirty tags (low quality tags) from raw data. Distinct tags are different kinds of tags. Unambiguous tags are clean tags remaining after the removal of tags mapped to reference sequences from multiple genes. C represents the control group; P represents the PEG-treated group.

Analysis of Tag Mapping {#s3b}
-----------------------

A reference gene database that included 94,878 sequences of the *R. soongorica* unigene was preprocessed for tag mapping. Among the sequences, genes with a CATG site accounted for 77.89%. To obtain the reference tags, all CATG+17 tags were used as gene reference tags. Finally, 218,791 total reference tag sequences with 182,184 unambiguous reference tags were obtained. Approximately 40.44% and 38.75% of the distinct clean tags were mapped unambiguously in the unigene database, and 47.94% and 50.2% of the distinct clean tags were not mapped to the unigene virtual tag database in libraries C and P, respectively ([Table 1](#pone-0094277-t001){ref-type="table"}).

The sequencing saturation was analyzed in the two libraries to estimate whether the sequencing depth was sufficient for the transcriptome coverage. The genes that were mapped by all clean tags and unambiguous clean tags increased with the total number of tags. However, when the sequencing counts reached three million tags or higher, the number of detected genes was saturated ([Figure S1](#pone.0094277.s001){ref-type="supplementary-material"}). Given that Illumina sequencing can distinguish transcripts originating from both DNA strands and the strand-specific nature of the sequencing tags obtained, we found approximately 27% and 26% of distinct tags, which were mapped insense transcripts of libraries C and P, respectively ([Figure S2](#pone.0094277.s002){ref-type="supplementary-material"}). Approximately 25% and 24% of distinct tags perfectly matched antisense transcripts, which suggested that antisense genes play important roles in the transcriptional regulation of drought response in *R. soongorica*.

Drought-caused Changes in the Global Gene Expression in *R. Soongorica* {#s3c}
-----------------------------------------------------------------------

To obtain the transcriptional changes in PEG-treated *R. soongorica*, a rigorous algorithm method was applied to identify DEGs from normalized DGE data using pairwise comparison between the two groups (C and P). Results showed that 1,325 genes had FDRs ≤0.001, and \|log~2~Ratio\| ≥1, which were declared to be the DEGs between the control and PEG-treated plants. Among these genes, 379 (28.6%) genes were up-regulated and 946 (71.4%) were down-regulated in response to drought stress. The expression of most genes was suppressed after drought treatment.

Functional annotation was performed to assign the genes of *R. soongorica* with GO terms. The main GO terms included "cellular component", "molecular function", and "biological process". A total of 535 DEGs were annotated as "cellular component". The categories of "cell part" and "cell" had the same proportion (88.0%), followed by "intracellular" and "intracellular part", which accounted for 78.5% and 77.9% of assignments, respectively ([Table S1](#pone.0094277.s003){ref-type="supplementary-material"}). In addition, a total of 570 DEGs were annotated to "molecular function" terms. The categories "catalytic activity", "organic cyclic compound binding", and "heterocyclic compound binding" were prominent with 60.5%, 25.3% and 25.3% of assignments, respectively ([Table S1](#pone.0094277.s003){ref-type="supplementary-material"}). Furthermore, a total of 521 DEGs were annotated to "biological process" terms, in which 71.4% and 63.6% of the assignments belonged to "metabolic process" and "cellular process", respectively, followed by "cellular metabolic process" and "primary metabolic process" ([Figure 2](#pone-0094277-g002){ref-type="fig"}).

![Gene classification based on gene ontology (GO) for DEGs in C and P libraries of *R. soongorica*.\
Only the biological processes were used for GO analysis. C represents the control group; P represents the PEG-treated group.](pone.0094277.g002){#pone-0094277-g002}

To characterize the functional consequences associated with drought response, a pathway analysis of the DEGs based on the KEGG database was performed. The results indicated that genes participating in the biosynthesis of secondary metabolites, plant hormone signal transduction, plant-pathogen interaction, amino sugar and nucleotide sugar metabolism, and ubiquitin-mediated proteolysis were differentially expressed in PEG-treated group ([Table S2](#pone.0094277.s004){ref-type="supplementary-material"}).

Genes Associated with the Major Functional Group were Affected by PEG-induced Drought Stress {#s3d}
--------------------------------------------------------------------------------------------

The global gene expression results showed that a great number of genes were significantly affected by PEG-induced drought stress. We classified the DEGs into three groups.

Group A was composed of some regulatory proteins, such as transcription factors, protein kinase, and other signaling molecules ([Table 2](#pone-0094277-t002){ref-type="table"}). Twenty genes encoding protein kinase, including CBL-interacting protein kinase (CIPK), threonine protein kinase, receptor-like protein kinase (RPK), mitogen-activated protein kinase (MAPK) cascade, phosphoenolpyruvate carboxylase kinase, and other kinases, were remarkably changed under drought stress. The results showed that six genes were up-regulated, and 14 genes were down-regulated. Meanwhile, a total of 14 drought-activated or repressed transcription factor genes were identified from several different families including WRKY, NAC, MYC, TCP, and bZIP. In our results, five transcription factors were identified as up-regulated, whereas nine transcription factors were repressed. The expression profiles exhibited that the majority of regulatory proteins were repressed after 12 h of drought treatment.

10.1371/journal.pone.0094277.t002

###### Selected DEGs between the C and P libraries.

![](pone.0094277.t002){#pone-0094277-t002-2}

  Funtional group                            Unigene ID                          Gene annotation                      TPM-C    TPM-P    Log~2~Ratio (C/P)
  ------------------------------------- -------------------- ------------------------------------------------------- -------- -------- -------------------
  **Protein kinases**                                                                                                                  
                                           Unigene 9485_A                CBL-interacting protein kinase                2.99    12.87          2.11
                                         CL2201. Contig1_3               CBL-interacting protein kinase                25.9    59.45          1.20
                                         CL4525. Contig1_A                  threonine protein kinase                   4.98    15.62          1.65
                                           Unigene 2338_A                   threonine protein kinase                   2.85     7.96          1.48
                                          Unigene 14518_A                serine/threonine-protein kinase              78.13    35.73          -1.13
                                          Unigene 32052_3                 receptor-like protein kinase                 2.56     0.01           -8
                                          Unigene 39297_A                 receptor-like protein kinase                 2.13     0.01          -7.73
                                          Unigene 52174_A                 receptor-like protein kinase                 6.26     1.3           -2.27
                                         CL2777. Contig2_3                receptor-like protein kinase                36.43     14.9          -1.29
                                         CL6483. Contig1_3                receptor-like protein kinase                16.94     7.67          -1.14
                                          Unigene 43240_A                 receptor-like protein kinase                76.71    38.19          -1.01
                                          Unigene 21624_3        mitogen-activated protein kinase kinase kinase        7.26     1.59          -2.19
                                          CL1583.Contig2_A       mitogen-activated protein kinase kinase kinase        7.26     1.59          -2.19
                                           Unigene 9581_3              mitogen-activated protein kinase 2             25.62     10.7          -1.26
                                           Unigene 633_3                  AMP-activated protein kinase                24.76     10.7          -1.21
                                          Unigene 21618_3     calcineurin B-like protein-interacting protein kinase   25.33    12.01          -1.08
                                          Unigene 10286_A            phosphoenolpyruvate carboxylase kinase            0.57     4.63          3.02
                                          Unigene 27156_3                 pyruvate, phosphate dikinase                35.44     85.2          1.27
                                         CL4464. Contig2_3                       pyruvate kinase                      85.82    29.36          -1.55
                                          Unigene 50230_A                       histidine kinase                      40.28    17.94          -1.67
  **Transcription factors**                                                                                                            
                                          Unigene 30342_3                  WRKY transcription factor 7                 4.7     14.75          1.65
                                          Unigene 10063_3                WRKY transcription factor 6--1                1.99     0.01          -7.64
                                          Unigene 27699_3                          MYC protein                        14.23    47.74          1.75
                                         CL7944. Contig1_3                NAC domain-containing protein               14.23     4.48          -1.67
                                          Unigene 31521_3                 NAC domain-containing protein               12.52     4.19          -1.58
                                          Unigene 25674_3                transcription factor TCP20-like              13.66     33.7          1.30
                                         CL7341. Contig1_3               transcription factor TCP20-like              45.26    16.92          -1.42
                                          Unigene 22142_3                transcription factor TCP7-like               31.88    15.04          -1.08
                                          Unigene 21994_3                      DNA binding protein                     6.55    14.32          1.13
                                          Unigene 23673_3                 transcription factor bZIP124                51.24    20.83          -1.29
                                           Unigene 1336_3                  transcription factor bZIP12                 4.7     11.57          1.29
                                          Unigene 28584_3                  transcription factor UNE12                 37.86    12.01          -1.66
                                         CL8097. Contig1_A                 transcription factor UNE12                  8.11     3.18          -1.35
                                         CL5168. Contig1_3                 transcription factor UNE12                 12.95     5.35          -1.28
  **Flavonoids biosynthesis**                                                                                                          
                                          CL93. Contig1_A               UDP-glucose: glucosyltransferase              22.06    50.19          1.19
                                         CL1676. Contig1_A           Flavonol sulfotransferase-like protein            3.56     9.4           1.40
                                         CL2370. Contig2_3                  flavonoid 3′-hydroxylase                  10.39     25.6          1.30
                                           Unigene 1473_3                   flavonoid 3′-hydroxylase                  84.96    34.72          -1.29
                                          Unigene 19573_A                  dihydroflavonol 4-reductase                 1.14    17.65          3.95
  **LEA proteins**                                                                                                                     
                                         CL4558. Contig1_3             late embryogenesis abundant protein             5.84    12.73          1.12
                                           Unigene 442_3           late embryogenesis abundant protein group 9         7.54    22.57          1.58
  **Heat shock proteins**                                                                                                              
                                          Unigene 18824_A            mitochondrial small heat shock protein            3.84    10.99          1.51
                                         CL12936. Contig1_A        heat shock protein 17.5 cytosolic class II          1.85     8.25          2.15
  **aquaporin**                                                                                                                        
                                          Unigene 39258_A                           aquaporin                          1.71     9.69          2.50
                                          Unigene 19242_A                           aquaporin                         426.96   176.77         -1.27
  **Proline transporter**                                                                                                              
                                         CL1109. Contig2_A                     proline transporter                     7.54    27.34          1.86
  **Cold shock protein**                                                                                                               
                                           Unigene 233_3                       cold shock protein                     11.53    58.44          2.34
  **low temperature-induced protein**                                                                                                  
                                          Unigene 24454_A                low temperature-induced protein               0.01     2.03          7.67
  **Dehydration-induced proteins**                                                                                                     
                                           Unigene 922_A                  dehydration-induced proteins                10.67    21.99          1.04
  **universal stress protein**                                                                                                         
                                          Unigene 25385_3                   universal stress protein                  93.22    263.56         1.50
  **Resistance protein**                                                                                                               
                                         CL8753. Contig4_3                 disease resistance protein                  1.85     7.23          1.97
                                         CL6875. Contig2_3                 disease resistance protein                  0.01     2.03          7.67
                                         CL5118. Contig1_3                 disease resistance protein                  4.7      1.01          -2.22
                                          Unigene 26781_3       natural resistance-associated macrophage protein       1.28    16.49          3.69
                                          Unigene 28389_3                  nematode resistance protein                36.86    77.68          1.08
  **defensin precursor**                                                                                                               
                                          Unigene 26068_3                      defensin precursor                     20.92    87.08          2.06
  **chitinase**                                                                                                                        
                                          CL630. Contig2_3                     class IV chitinase                      3.7     67.84          4.20
                                         CL5695. Contig1_A                          chitinase                          6.4      41.8          2.71
                                          Unigene 30551_3                      class IV chitinase                     13.52     6.08          -1.15
                                          CL630. Contig3_3                     class IV chitinase                     60.91    27.63          -1.14

C represents the control group; P represents the PEG-treated group. The DEGs was selected through FDR ≤0.001 and \|log~2~Ratio\| ≥1.

Group B was composed of 13 functional protein genes, which included genes involved in flavonoid biosynthesis (6 DEGs), LEA protein (2 DEGs), small heat shock protein (sHSP) (2 DEGs), aquaporin (2 DEGs), and proline transporters (1 DEG) ([Table 2](#pone-0094277-t002){ref-type="table"}). Contrary to the expression of regulatory protein genes, most of the functional protein genes were up-regulated after 12 h of drought treatment.

Group C was composed of 14 proteins that were annotated as low-temperature-induced protein (1 DEG), dehydration-induced protein (1 DEG), defensing precursor (1 DEG), resistance protein (4 DEG), universal stress protein (1 DEG), and proteins involved in chitinase biosynthesis (4 DEGs) ([Table 2](#pone-0094277-t002){ref-type="table"}). Similar to the alteration of functional protein genes, the expression of these 14 genes also exhibited a general upward trend.

Discussion {#s4}
==========

Plants, which are sessile organisms, try to respond and adapt to stress by altering their gene expression patterns when they are subjected to adverse environments [@pone.0094277-Hazen1]. Therefore, the genes that showed different expressions are probably involved in stress response. Shinozaki and Yamaguchi-Shinozaki [@pone.0094277-Shinozaki1] reported that the drought-inducible genes identified can be classified into regulatory and functional proteins. The regulatory proteins include various transcription factors, protein phosphatases, protein kinases, and other signaling molecules that participate in further regulation of signal transduction and downstream gene expression. The functional proteins are composed of molecules, such as chaperones, water channel proteins, LEA proteins, mRNA-binding proteins, sugar and proline transporters, detoxification enzymes, and various proteases. In the present study, the DEGs between the leaves of drought-stressed and unstressed *R. soongorica* were investigated to conduct a global analysis of the transcriptomic response, which can facilitate our understanding of drought tolerance mechanisms. In our results, 1325 DEGs with 379 up-regulated and 946 down-regulated genes were detected. Our results are similar to that of Fan et al. [@pone.0094277-Fan1], who also reported that down-regulated genes were more abundant than up-regulated in soybean leaves after drought exposure. Functional annotation analysis results indicated that many drought-inducible genes with various functions were identified in *R. soongorica*, including a number of regulatory proteins, functional proteins, and other inducible proteins in response to drought stress.

Regulatory proteins are involved in signal transduction pathways and in controlling the expression of functional genes in stress responses. In higher plants, many regulatory genes, such as those encoding protein kinases and transcription factors, are induced by environmental signals or stresses [@pone.0094277-Shinozaki2]. In this study, 20 DEGs encoding protein kinases were detected after drought treatment. Various protein kinases including CIPK protein, threonine protein kinase, MAPK, calcium-dependent protein kinases, and RPK are important regulatory proteins in the regulation of some stress-inducible genes [@pone.0094277-Zhu1]. Several experiments have demonstrated that protein kinase signaling conferred plant resistance to stress by regulating downstream functional genes [@pone.0094277-Boller1], [@pone.0094277-Lacombe1]. In addition, protein kinase also functions in regulating transcription factors [@pone.0094277-Tena1]. A total of 14 genes encoding transcription factors, including WRKY, NAC, MYC, TCP, and basic region leucine zipper (bZIP), were also significantly differentially expressed. Transcription factors are important in regulating plant responses to biotic and abiotic stresses [@pone.0094277-Chen3]. Altering the expression of certain transcription factors can greatly change the expression of a number of stress-related target genes that influence plant stress tolerance in transgenic plants [@pone.0094277-YamaguchiShinozaki1], [@pone.0094277-Shinozaki1]. Thus, the alteration of protein kinases and transcription factors may function in the signal transduction process by regulating the expression of various functional genes in response to drought stress.

Aside from signaling transduction, plants adapt to adverse environment through the activation or up-regulation of genes that directly function in stress resistance [@pone.0094277-Hazen1]. Flavonoids are a class of secondary metabolites with diverse functions in growth, development, reproduction, and are also involved in diverse stress responses in plants. Flavonoids accumulate rapidly and play a protective role when plants are exposed to abiotic stress [@pone.0094277-Castellarin1]. In recent years, several genes encoding important enzymes, such as phenylalanine ammonia-lyase, F3H, dihydroflavonol 4-reductase (DFR), and anthocyanin synthase, which are involved in the flavonoid biosynthetic pathway, were found to be up-regulated in response to different kinds of stresses, such as salinity, UV-B irradiation, water deficit, and nitrogen stress [@pone.0094277-Castellarin1]--[@pone.0094277-Han1]. Increased expression of genes involved in the flavonoid pathway was detected in our previous study on *R. soongorica* [@pone.0094277-Liu4]. In the present study, most DEGs encoding key enzymes involved in flavonoid biosynthesis were also up-regulated, especially DEGs that encode UFGT, DFR, and F3′H. These findings are consistent with the results of other studies conducted on birch, potato, and rice after stress exposure [@pone.0094277-Morales1]--[@pone.0094277-Lenka1]. The up-regulated genes may result in accumulation of flavonoids. Flavonoids may serve as signaling molecules in signaling cascades. Flavonoids also affect cellular function and modulate gene expression in response to adverse conditions [@pone.0094277-Williams1]. Agati et al. [@pone.0094277-Agati1] found that flavonoids directly scavenged the ROS brought about by stress. Therefore, up-regulation of flavonoid pathway genes may enhance the protective role of *R. soongorica*, thus leading to better drought tolerance.

Secondary metabolites and stress-responsive proteins, such as LEA proteins, accumulate naturally in the seeds of some desiccation-tolerant plants, and are also induced in vegetative tissues when subjected to water-limited conditions [@pone.0094277-Battaglia1], [@pone.0094277-Tunnacliffe1]. LEA proteins, largely attributed to their extreme hydrophilic nature, have also been implicated in detoxification, sequestration of ions, maintenance of protein or membrane structure, binding water, and operation as molecular chaperones during dehydration [@pone.0094277-Bray1]. In this study, two DEGs were identified as up-regulated LEA proteins. Our results are in agreement with the findings of Gechev et al. [@pone.0094277-Gechev1], who found high expression of LEA genes during drought and desiccation in leaves of *Haberlea rhodopensis*. We concluded that dehydration-induced expression of LEA transcripts is a common response in resurrection and desiccation sensitive plants [@pone.0094277-Challabathula1]. Improved resistance to drought and salinity was also observed in transgenic rice when LEA was overexpressed [@pone.0094277-Xu1]. Expression of LEA proteins in yeast confers improved resistance to high salinity and freezing [@pone.0094277-Imai1]. These studies demonstrated that LEA proteins functions in binding water, maintenance of protein or membrane structure, or operation as molecular chaperones to enhance tolerance to dehydration, whereas few additional clues to the mechanism of protein function are provided.

Plant HSPs facilitate protein folding or assembly under diverse developmental and adverse environmental conditions. In this study, two small HSPs (sHSPs) encoding DEGs showed up-regulation under drought stress. One sHSP was annotated to mitochondrial sHSP (MT-sHSP), and the other was annotated as cytosolic sHSP17.5. Although the precise functional mechanism of sHSPs is still unclear, in vivo studies have shown that sHSPs function as molecular chaperones in stressful conditions and in normal development [@pone.0094277-Sun1]. Overproduction of sHSP increased drought tolerance in transgenic rice seedlings [@pone.0094277-Sato1]. MT-sHSP protects the electron transport complex I and prevents damage by ROS [@pone.0094277-Banzet1]. Lee et al. [@pone.0094277-Lee1] found that transgenic plants with over-expression of MT-sHSP efficiently maintained membrane stability by scavenging ROS, and were conferred with enhanced tolerance to salinity. Ectopic expression of cytosolic sHSP 17.1 in *Arabidopsis thaliana* exhibited higher resistance to drought [@pone.0094277-Zhang1]. In a study conducted by Sato et al. [@pone.0094277-Sato2], sHSP 17.5 showed molecular chaperone activity in vitro and represented an example of an HSP capable of protecting cells against thermal extremes. Therefore, these results indicated that the sHSPs may play roles in scavenging ROS or functioning as molecular chaperone to respond to drought stress.

The DEGs encoding aquaporins and proline transporter were also identified. Aquaporins are channel proteins of intracellular and plasma membranes that function in transporting water, small solutes (urea, boric acid, and silicic acid), and gases (ammonia and carbon dioxide). A specific role for aquaporins in embolism refilling and recovery of stem axial conductance after drought was proposed in grapevine [@pone.0094277-Lovisolo1]. In our results, two aquaporin genes were identified. One aquaporin gene was up-regulated, whereas the other gene was down-regulated. An overall tendency of down regulation for aquaporin gene was observed. However, up regulation of certain aquaporin transcripts was also detected in rice and Arabidopsis leaves [@pone.0094277-Alexandersson1], [@pone.0094277-Guo1]. Therefore, transcriptional control of aquaporins in drought-stressed leaves appears to be more complex [@pone.0094277-Maurel1]. Interestingly, transport of proline may also be affected by water deficit. A gene encoding proline transporter was up-regulated under drought stress in the present study, which is in agreement with the result of Rentsch et al. [@pone.0094277-Rentsch1], who also observed the induction of a specific proline transporter in Arabidopsis under water deficit. The DEGs encoding aquaporins and proline transporter in *R. soongorica* may serve as an adaptive strategy by regulating water and proline distribution throughout the whole plant.

Moreover, in addition to two dehydration-induced 19 homologs, proteins that can be induced by other stresses were also detected in the study and annotated to low-temperature-induced (LT) protein, cold shock protein (CSP), defensin precursor, chitinase, disease resistance protein, natural resistance-associated macrophage protein, nematode resistance protein, and universal stress protein. LT proteins and CSPs are commonly induced by cold stress. Various genes are induced by both drought and cold stress, suggesting the existence of crosstalk between drought and cold-stress signaling pathways [@pone.0094277-Shinozaki1]. Defensins, chitinase, disease resistance proteins, natural resistance-associated macrophage protein, and nematode resistance protein were all found to be involved in the defense against phytopathogens [@pone.0094277-Thomma1]--[@pone.0094277-Williamson1]. Therefore, the identification of the genes in this study suggested that a crosstalk between drought and pathogen attack. The universal stress proteins can confer the ability to respond and adapt to environmental changes [@pone.0094277-Isokpehi1]. The identification of all these genes suggests crosstalk between responsive genes or pathways and multiple abiotic or even biotic stresses [@pone.0094277-Zhou1]. Moreover, these identified genes may provide gene resources for understanding the drought tolerance mechanism of *R. soongorica*.

Among the identified DEGs, most regulatory proteins were down-regulated, whereas the majority of functional proteins were up-regulated. The different expression patterns shown by the two groups are interesting and worthy to be explored. Therefore, we tried to discuss the probable interaction of regulatory proteins and functional proteins. Most gene regulators exhibit either antecedent or simultaneous change in the expression level when compared with their targets [@pone.0094277-Yu1]. One reason is that it takes time for the regulator gene to express its protein product and affect (directly or indirectly) the transcription of its target gene [@pone.0094277-Ferrell1]. Other reasons may have been caused by the different mRNA half-lives of the regulator and target genes, and the different profiles between regulator genes and functional genes may also indicate a feedback loop, where a target gene can regulate its regulator [@pone.0094277-Zou1]. Given that signaling pathways are complex dynamic events that occur over time, single time point expression profiles in our study are insufficient to elucidate temporal events. The different expression profiles of regulatory and functional proteins observed in this study are interesting. Expression profiling experiment with a series of time points should be performed to elucidate this problem in further studies.

Conclusions {#s5}
===========

Based on the putative function analysis of annotated DEGs, DEGs are classified into three groups. The first group, which includes protein kinases (CIPK and threonine protein kinase) and transcription factors (WRKY, MYC and TCP), regulates signal transduction. The second group includes genes involved in the biosynthesis of flavonoids in scavenging ROS, sHSP, LEA proteins, aquaporins, and proline transporter that represents functional genes that directly enhance drought-stress tolerance. The sHSP maintains protein or membrane structure and act as molecular chaperones that protect cells from injury under drought stress. Aquaporins and proline transporter regulates water or proline distribution throughout the whole plant. The third group contains proteins induced by other stress factors, which suggests that a crosstalk between different stresses. In conclusion, the tolerant ability of *R. soongorica* may be gained through effective signal transduction and enhanced protection of functional proteins to reestablish cellular homeostasis.

Supporting Information {#s6}
======================

###### 

Sequencing saturation analysis of the two libraries. C represents the control group; P represents the PEG-treated group.

(TIF)

###### 

Click here for additional data file.

###### 

Mapping of distinct clean tags in the two DGE libraries. C represents the control group; P represents the PEG-treated group.

(TIF)

###### 

Click here for additional data file.

###### 

Gene classification based on gene ontology (GO) for DEGs in C and P libraries of *R. soongorica*. C represents the control group; P represents the PEG-treated group.

(XLSX)

###### 

Click here for additional data file.

###### 

KEGG analysis for DEGs in C and P libraries of *R. soongorica*. C represents the control group; P represents the PEG-treated group.

(XLSX)

###### 

Click here for additional data file.

We thank the technical support for Illumina sequencing and initial data analysis from Beijing Genome Institute at Shenzhen, P. R. China, and appreciate the constructive and helpful suggestions from the two anonymous reviewers.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.

[^2]: Conceived and designed the experiments: XRL YBL. Performed the experiments: MLL BC. Analyzed the data: MLL. Contributed reagents/materials/analysis tools: XRL YBL XFM. Wrote the paper: MLL YBL.
