Comparative transcriptome profiling of a desert evergreen shrub, Ammopiptanthus mongolicus, in response to drought and cold stresses

Background The molecular mechanisms involved in plant tolerance to either drought or cold have been extensively studied in many plant species. However, few studies have focused on their comparisons especially using non-model plants with strong tolerance to both stresses. Ammopiptanthus mongolicus (Maxim. ex Kom.) Cheng f. is the only evergreen broadleaf shrub grown in the central Asian desert and it has very strong cold and drought tolerance. To provide further insights into plant tolerance, the transcriptome profiles of drought- and cold-treated A. mongolicus seedlings were analyzed using Illumina technology and differentially expressed genes (DEGs) were compared. Results A comprehensive transcriptome of A. mongolicus was sequenced using pooled mRNA extracted from drought-, cold-stressed and unstressed seedlings as well as leaves from naturally grown shrub. These sequences were assembled into 86058 unigenes, of which 51014 unigenes had an annotated function and 2440 encoded transcription factors (TFs). Transcriptome profiles were analyzed in A. mongolicus seedlings after drought and cold treatments at three time points (2, 8 and 24 h). Between 3917 and 6102 unigenes were identified as DEGs at a single time point in both stresses. Among these DEGs 2028 and 2026 DEGs were common across the three time points of drought and cold treatments respectively, and 971 DEGs were co-regulated by both stresses. Functional enrichment analyses identified many common or specific biological processes and gene sets in response to drought and cold stresses. The most pronounced findings are that flavonoid biosynthesis genes were enriched in the DEGs co-up-regulated by both stresses; while membrane protein genes and genes related to chloroplast were abundant in the DEGs specifically up-regulated by drought or cold, respectively. Furthermore, the DREB, ERF, NAC and WRKY TFs were predominantly co-up-regulated by both stresses. Conclusions The present study provides the most comprehensive transcriptome resource and the first dynamic transcriptome profiles of A. mongolicus under drought and cold stresses. This information will deepen our understanding of plant tolerance to drought and cold. The up-regulated DEGs will be valuable for further investigations of key genes and molecular mechanisms involved in the adaptation of A. mongolicus to harsh environments. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-671) contains supplementary material, which is available to authorized users.


Background
Drought and low temperature are two major environmental stresses that greatly affect plant growth and crop production worldwide. This is particularly the case in arid and cold areas. To improve the tolerance of crops and other economic plants to these stresses, it is essential to know how plants respond to these stresses and which genes and biological pathways are involved in the stress tolerance. Over the past decades, thousands of genes and dozens of metabolic and signaling pathways have been identified in plants during drought and cold stresses [1][2][3][4]. Some of these genes have been confirmed to significantly improve plant tolerance to drought, cold and/or salt stresses [5,6]. In addition, many genes have been found to be commonly or specifically regulated under drought and cold in Arabidopsis [7], rice [8] and Brachypodium [9] at transcriptional level. However, similar studies are very limited in other plants, and the molecular mechanisms controlling plant tolerance to both drought and cold remain largely unknown.
Ammopiptanthus is an endangered survivor from the Tethys in the Tertiary Period, and is the only evergreen broadleaf shrub grown in the central Asian desert. The genus Ammopiptanthus (Leguminosae) comprises of two species: A. mongolicus (Maxim. ex Kom.) Cheng f. and A. nanus (M. Pop.) Cheng f. Both species are diploid with 18 chromosomes and have high stress tolerance [10]. The habitats of A. mongolicus are marked by arid climate with an annual precipitation less than 200 mm but a mean annual evaporation up to 3000 mm. Its natural distribution areas are also characterized as sandy or stony soil with high salinity, intense ultraviolet irradiation, and seasonally extreme temperature from about -30°C in winter to more than 40°C in summer. The extreme tolerance of this species to harsh environments makes it invaluable for exploring key stress-tolerant genes and mechanisms, especially those involved in cold and drought tolerance.
The transcriptomes of A. mongolicus roots (treated by 20% PEG-6000 for 72 hours) and seedlings (cultured at normal condition or 4°C for 14 days) were recently sequenced by 454 pyrosequencing [11] and Illumina technology [12], respectively. These studies generated 29056 and 82795 unigenes and identified 32728 cold-regulated genes. However, transcriptome profiling of this species under drought stress and its comparison with cold stress have not been reported in the literatures.
In the present study, we firstly performed a comprehensive transcriptome sequencing of A. mongolicus using pooled mRNA extracted from drought-treated, cold-treated and control seedlings as well as young leaves of plants grown in desert in both summer and winter. We then applied RNA-Seq (RNA sequencing) to investigate seven cDNA libraries derived from the seedling samples exposed to drought and cold stresses for 2 h, 8 h and 24 h, and the non-stressed control seedlings respectively. Finally, we analyzed the DEGs and identified common and specific DEGs during drought and cold treatments. Our study provided the first transcriptome profile of A. mongolicus under drought stress and the comparison of it with the transcriptome profile under cold stress.

Results
Sequencing, de novo assembly and functional annotation of the A. mongolicus transcriptome Illumina sequencing of the pooled mRNA generated 47.29 million (M) clean paired-end reads with Q20 over 98%. De novo assembly of the clean reads resulted in 86058 unigenes with an average length of 756 bp and N50 length of 1279 bp (Table 1). Approximately 41.4% and 22.6% of these unigenes had a length ≥ 500 bp or ≥ 1000 bp, respectively ( Figure 1).
Of the 86058 sequences, 51014 (59.3%) unigenes had at least one significant match in the NCBI protein databases (E-value ≤ 1e-5). Distribution of the annotated unigenes in each database was shown in the Additional file 1. By searching against the PlantTFDB (http://planttfdb.cbi.edu.cn) (E-value ≤ 1e-5) in combination with annotations in the NR (non-redundant) and Swiss-Prot databases, 2440 unigenes were annotated as transcription factors (TFs), representing 4.8% of the annotated unigenes and covering 57 TF families (Table 2). Furthermore, 47486 unigenes had significant hits in the TIGR Plant Transcript Assemblies database (E-value ≤ 1e-5). Most of these hits were from Glycine max, one of its closest relatives, followed by Medicago truncatula, Lotus japonicus, Malus x, Glycine soja, and others. Surprisingly, only 919 (1.9%) and 787 (1.7%) unigenes had significant similarity to the sequences of Arabidopsis thaliana and Oryza sativa, respectively (Additional file 2).
Notably, nearly 1600 of the annotated unigenes were assigned to the class of "Function unknown" (Additional file 3).
To gain insights into function interactions among the unigenes, 26999 of the 86058 unigenes were assigned to 125 KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways. Of the top 10 pathways (Table 3), "Metabolic pathways" were mostly represented, followed by "Biosynthesis of secondary metabolites", "Plantpathogen interaction" and "Plant hormone signal transduction", etc. Unigenes encoding core components in the plant hormone signal transduction pathways were present, including ABA (abscisic acid), BR (brassinosteroid) and ethylene signaling pathways. These pathways were well-documented to participate in the response and adaptation of plants to stressful environments [13,14].

Identification and validation of differentially expressed genes (DEGs)
Transcriptome changes of the A. mongolicus seedlings under drought and cold stresses (2 h, 8 h and 24 h) and   Figure 3). Taken together, 9487 and 6524 unigenes were identified as DEGs with an absolute value of log 2 Ratio ≥ 1 and FDR (false discovery rate) ≤ 0.001 at one or more time points during drought and cold stresses, respectively. To validate the RNA-Seq results, 11 DEGs with different expression patterns were selected for RT-qPCR (real-time quantitative PCR) analysis. These genes encode NAC (NAM, ATAF and CUC) and DREB (dehydration responsive element binding) TFs, HSP (heat shock protein), GLP (germin-like protein), TPP (trehalose-6-phosphate phosphatase) and CCD (carotenoid cleavage dioxygenase). Nine (81.8%) of these genes showed similar expression patterns to those detected by the RNA-Seq data ( Figure 4), while other two genes had similar expression trends at Figure 2 Histogram of GO classification of the annotated unigenes. Unigenes with best BLAST hits were aligned to GO database. In total, 23834 unigenes were assigned to at least one GO term and were grouped into three main GO categories and 44 GO terms. Left Y-axis represented the percentages of unigenes in each main category. Right Y-axis indicated the numbers of unigenes in each GO term. DEGs were common at two of the three time points ( Figure 5). The top 100 DEGs commonly regulated by drought stress were listed in the Additional file 4. Among the common DEGs up-regulated at the three time points, genes known to respond to drought or other stresses were largely represented, such as those encoding ripening-related proteins, LEA (late embryogenesis abundant) proteins, peroxidases, transporters, enzymes in the flavonoid biosynthetic pathways, protein kinases, and ethylene receptors. Many TFs in multiple families were also up-regulated, including the AP2/EREBP (APETALA2/ethylene-responsive element binding protein), NAC, WRKY, MYB (myeloblastosis oncogene), bHLH (basic helix-loop-helix), C2H2 (Cys2/His2-type zinc finger protein 2), bZIP (basic-domain leucinezipper protein), C3H (Cys3/His-type zinc finger protein), GRAS (GAT, RGA and SCR), TCP (TB1, CYC and PCF), COL (constans-like), Dof (DNA binding with one finger), and GRF (general regulatory factor) families, etc. Approximately 17.2% (134) of the common DEGs were functional unknowns. Similarly, the common DEGs down-regulated at the three time points of drought stress had a wide range of functions. For example, genes encoding aquaporin NIP, polyol transporter, gibberellin, small heat shock proteins, calcium-binding proteins, and some bHLH TFs appeared in this group. Interestingly, the unknowns in the down-regulated DEGs (397,  33.5%) were almost twice as those in the commonly upregulated DEGs. Based on GO enrichment analysis, a large number of the DEGs commonly up-regulated in drought were mostly enriched in two similar biological processes: "Response to stimulus" and "Response to stress", followed by "Organic substance transport", "Cellular aromatic compound metabolic process" and "Protein phosphorylation". In contrast, the commonly down-regulated DEGs were significantly enriched in the biological processes "Monocarboxylic acid metabolism", "Small molecule biosynthesis" and "Nucleobase transport". In the cellular component and molecular function categories, the commonly up-and downregulated DEGs were significantly enriched in similar terms, such as "Membrane" and "Oxidoreductase activity" (Additional file 5). KEGG pathway enrichment analysis revealed that genes involved in microbial metabolism in diverse environments and some secondary metabolite degradation were significantly enriched in both the up-and down-regulated DEGs. However, genes in the pathways for nitrogen and phenylalanine metabolisms were specifically enriched in the up-regulated DEGs; while genes involved in Figure 4 Validation of expression patterns of nine DEGs by RT-qPCR assay. Nine DEGs (differentially expressed genes) showed similar expression patterns between RNA-Seq data (A) and RT-qPCR assay (B). These genes included six genes up-regulated under drought and cold stresses: AmNAC1 (NAM, ATAF and CUC), AmNAC2, AmNAC11, AmDREB2C (dehydration responsive element binding 2C), AmDREB2 and AmTPP (trehalose-6-phosphate phosphatase, not significantly down-regulated in C1 library); two genes down-regulated under both stresses: AmHSP (heat shock protein) and AmCCD (carotenoid cleavage dioxygenase); and one gene with different expression trends: AmGLP (germin-like protein). D1, D2 and D3 represented 2, 8 and 24 h after drought treatment; C1, C2 and C3 represented 2, 8 and 24 h after cold treatment. The total RNAs used for the transcriptome profile analyses using RNA-Seq were used in the RT-qPCR assay, including three biological replicates. The relative expression levels of the selected genes were calculated using the 2 -ΔΔCT method. Error bars represented the standard deviation of the mean expression values. metabolic pathways and the pathways for starch and sucrose metabolism as well as carotenoid biosynthesis were only enriched in the down-regulated DEGs (Additional file 6).
As expected, many regulatory genes were specifically up-regulated at the first time point (2 h), such as genes encoding some types of protein kinases, calcium-binding proteins, TFs, as well as enzymes or components in plant hormone (especially ABA) biosyntheses or signaling pathways. After the treatment was prolonged (8 h to 24 h), many of these regulatory genes were no longer differentially expressed. Instead, a large number of functional genes were up-regulated, such as genes encoding osmoprotectant synthetases, dehydration-protective proteins, PR (pathogenesis-related) proteins, hydrolases, and detoxification enzymes.

Differentially expressed genes in response to cold stress
Between 3917 and 4413 DEGs were identified at a single time point during cold treatment, and 2026 DEGs were common across the three time points, including 851 with up-regulation and 1174 with down-regulation ( Figure 5). The top 100 DEGs commonly regulated by cold stress were listed in the Additional file 7. Similar to drought stress, many commonly up-regulated DEGs were annotated to well-known stress-or cold-induced genes. Some of these genes are identical to the drought-induced genes; while others encode different types of proteins, such as fatty acid desaturases, chitinases and proteins or enzymes located in chloroplast. The common DEGs down-regulated at three time points of cold stress also encode various functional and regulatory proteins. Notably, 16.8% (143) and 40.3% (473) of the commonly up-or down-regulated DEGs were unknowns, respectively.
A large number of the DEGs up-regulated under cold stress were significantly enriched in the biological processes "Response to stimulus", "Response to stress", "Defense response" and "Cellular aromatic compound metabolic process". In the down-regulated DEGs, the mostly represented GO term of biological process was also "Response to stimulus", followed by "Lipid biosynthetic process", "Organic acid biosynthetic process" and "Lipid localization". Interestingly, in the cellular component category, the up-regulated DEGs were mostly enriched in "Plastid" and "Chloroplast"; while the downregulated DEGs were not significantly enriched in any terms of this category (Additional file 8). KEGG pathway analysis revealed that the metabolic pathways and the pathways for secondary metabolite biosynthesis or degradation were significantly enriched in both the upand down-regulated DEGs. However, the pathways involved in phenylalanine metabolism and photosynthesis were exclusively enriched in the up-regulated DEGs. In contrast, the pathways related to fatty acid metabolism, carotenoid biosynthesis and several amino acid degradations were only enriched in the down-regulated DEGs (Additional file 9).

Differentially expressed genes co-regulated during drought and cold stresses
Among the 9487 drought-regulated and 6524 coldregulated DEGs with significant expression changes (absolute value of log 2 Ratio ≥ 1 and FDR ≤ 0.001) at one or more time points, 4515 DEGs were co-regulated by both stresses, accounting for 47.6% and 69.2% of the total DEGs, respectively. Of the co-regulated DEGs, 971 DEGs had significant expression changes at all droughtand cold-treated time points, including 334 with up-regulation, 618 with down-regulation and 19 with complex regulation patterns. The results clearly showed an overlap between the DEGs in response to drought and cold. The top 100 co-regulated DEGs were listed in the Additional file 10.
Of the 971 co-regulated DEGs, 69 DEGs (33 up-and 36 down-regulated) were annotated as TFs, including the AP2/EREBP, NAC, WRKY, MYB, bHLH, C2H2, bZIP, GRAS, C3H, Dof, TCP, COL and GRF families. In addition, 46 co-regulated DEGs were annotated as protein kinases and phosphatases, and 17 co-regulated DEGs as synthetases or components in plant hormone biosynthesis or signaling. Besides, many stress-related functional genes appeared in the co-up-regulated DEGs, such as those encoding the enzymes in flavonoid biosynthetic pathways, ripening related proteins, thaumatin-like or osmotin-like proteins, transporters, cell wall-associated proteins, and some types of PR proteins.
The co-up-regulated DEGs were significantly aggregated into two GO biological processes: response to stress and secondary metabolite (mainly aromatic compounds or phenylpropanoids) biosynthesis or metabolism. However, the co-down-regulated DEGs were significantly enriched in carboxylic acid metabolic process and fatty acid or lipid biosynthetic process. In the molecular function category, the co-up-or co-down-regulated DEGs were mainly related to catalytic activity or oxidoreductase activity. Unexpectedly, both the up-and down-regulated DEGs were not significantly enriched in any cellular component terms (Table 5).
According to KEGG pathway analysis, the co-upregulated DEGs were significantly enriched in eight pathways. Most of these pathways are completely or highly overlapped and actually belong to two pathways: microbial metabolism in diverse environments and biosyntheses of secondary metabolites. In contrast, the codown-regulated DEGs were mainly assigned to metabolic pathways and starch/sucrose metabolic pathway. The DEGs involved in carotenoid biosynthesis and bisphenol degradation were also significantly enriched (Table 6). Furthermore, about 30% of the co-regulated DEGs were unknowns or unannotated genes. They may play unique roles in response and adaptation to drought and cold environments of A. mongolicus.
Differentially expressed genes specifically responding to drought or cold stress In total, 4972 and 2009 DEGs were specifically regulated in at least one time point during drought or cold stress, respectively. Of these DEGs, 1057 (445 up-regulated, 566 down-regulated and 46 complicatedly regulated) or 1055 (514 up-regulated, 541 down-regulated) had significant expression changes across all three time points under drought or cold stress, respectively. The top 100 drought-or coldspecific DEGs were listed in the Additional files 11 and 12.
Among the 1057 drought-specific DEGs, 91 DEGs encode TFs, 66 encode protein kinases and phosphatases, and 24 are involved in hormone or cell signaling. Compared to the co-regulated DEGs, more TF gene families were found in the drought-specific DEGs. In the 1055 cold-specific DEGs, however, fewer genes encode regulatory proteins. Some distinct differences appeared between the drought-and cold-specific DEGs. For instance, genes encoding trehalose-6-phosphate synthase, trehalose-6-phosphate phosphatase and early responsive to dehydration protein were only found in the DEGs specifically up-regulated in drought; while genes encoding acidic endochitinase, protochlorophyllide reductase, protein chloroplast import apparatus 2 and chloroplast omega-3 fatty acid desaturase appeared only in the DEGs specifically up-regulated in cold. In addition, about 23.7% of the drought-specific DEGs and 32.3% of the cold-specific DEGs were functional unknowns.
The DEGs specifically up-regulated in drought were mostly enriched in the GO terms "Response to stress", "Response to stimulus", "Membrane" and "Antioxidant activity" (Additional file 13). However, the DEGs specifically up-regulated in cold were mostly enriched in the terms "Aminoglycan catabolic process", "Plastid", "Chloroplast" and "pattern binding" (Additional file 14). On the other hand, the DEGs specifically down-regulated in drought were more enriched in the GO terms "Small molecule biosynthetic process", "Organic acid biosynthetic process" and "Iron ion binding" (Additional file 13); while the DEGs specifically down-regulated in cold were highly enriched in the GO terms "Lipid localization" and "Antioxidant activity" (Additional file 14). These results indicated that distinct molecular mechanisms were involved in A. mongolicus seedlings under drought compared to those in cold stress.

Discussion
A. mongolicus is an emerging model species for studying the response and tolerance of plants to extreme cold and drought environments. By using Illumina technology, 9309 up-regulated and 23419 down-regulated DEGs were previously identified in the cold-acclimated A. mongolicus seedlings [12]. Using an EST approach, 120 abiotic stressresponsive genes were identified in this species, and 82 genes were confirmed to be cold-and/or droughtinducible by qRT-PCR [15]. In the present study, we performed comparative transcriptome profiling of A. mongolicus in response to drought and cold and identified thousands of the DEGs regulated by each stress. In particular, 971 co-regulated DEGs by both stresses and about 1050 DEGs specifically regulated by a single stress were identified, suggesting common and distinct molecular mechanisms underlying the response to drought and cold environments. These results were similar to those found in Arabidopsis [7], rice [8] and Brachypodium [9] despite of some new information also observed in the present study. The co-regulated DEGs could represent a subset of basal drought and cold responsive genes in A. mongolicus and thus may play essential roles in its adaptation to adverse environments; while the specifically regulated DEGs might underline the distinct mechanisms of this species in response to different stresses. Functional enrichment analyses of these DEGs provided clues to the molecular bases of A. mongolicus combating drought and cold environments.
Flavonoids may contribute greatly to the adaptation of A. mongolicus to drought and cold environments Flavonoids constitute a large group of polyphenolic secondary metabolites in plants, which have antioxidant activity and are of prime importance for plant defense against pathogens and UV stress [16]. Recent studies have showed that flavonoids accumulated to high levels in high alpine and polar plants [17] and in A. mongolicus  [19,20]. In the present study, genes encoding almost all of these enzymes were not only moderately or highly expressed in A. mongolicus seedlings in normal condition (average RPKM is 25.3), but also were significantly up-regulated by drought and cold stresses. GO terms and pathways relevant to flavonoids were highly enriched in the co-up-regulated DEGs identified in both drought and cold stresses (Tables 5 and 6). The phenylpropanoid or flavonoid biosynthetic pathways are usually regulated by MYB and/or bHLH TFs [21,22]. Genes encoding these TFs were also included in the co-up-regulated DEGs identified in the present study. Although their exact functions need to be characterized, our results suggest that genes involved in flavonoid biosyntheses may be crucial in the adaptation of A. mongolicus to drought and cold environments. An increased production of reactive oxygen species is a common consequence in plants under abiotic stresses, which usually damage cellular membranes and other cellular components resulting in oxidative stress and eventually cell death [23]. It is likely that a high level of flavonoids may protect A. mongolicus cells from oxidative stress originated from drought and cold, and hence they contribute to stress tolerance.

Membrane proteins may play important roles in drought adaptation of A. mongolicus
In plants, membrane transport and perception systems have essential roles in maintaining cellular homeostasis under adverse environmental stresses through cell-to-cell and/or organ-to-organ communication [24,25]. An increased expression of multiple transporter and channel protein genes in response to drought has been found in some plant species [7,24,26]. Over-expressing some of these genes, such as aquaporin OsPIP2-2 [27] and ABCtransporter AtABCG36 [28], enhanced plant tolerance to drought and/or salt stress. In the present study, the GO term "Membrane" was highly enriched in the DEGs upregulated in drought (Additional files 5 and 13), and approximately 70 membrane protein genes were found with up-regulation under drought stress. Most of these genes were annotated to encode transporters and channel proteins, such as those for ions, sugars, organic acids, water, nitrate and phosphate, suggesting their involvement in drought adaptation of A. mongolicus.
Receptor-like protein kinases (RLKs), another group of membrane proteins, are important signaling components that mediate plant responses to developmental and environmental stimuli [25]. Recent studies have shown that some RLKs, such as Arabidopsis RPK1, GHR1 and AtCRK45, responded to drought, ABA, salt and/or cold stresses and they were implicated in ABA and BR signaling pathways [29,30]. Genes encoding cysteine-rich receptorlike kinases and LRR-RLKs (leucine-rich receptor-like protein kinases) were highly induced in chilling-treated Chorispora bungeana [31]. In the present study, several RLK genes, such as those encoding LRR-RLKs, lectin-like receptor kinase, leucine-rich repeat receptor-like tyrosineprotein kinase and G-type lectin S-receptor-like serine/ threonine-protein kinase, were up-regulated in droughtstressed A. mongolicus seedlings. Although their functions and detailed molecular mechanisms are not clear, an increased transcription of RLK genes suggests they may play important roles in drought adaptation of A. mongolicus.

Chloroplast plays a central role in cold adaptation of A. mongolicus
Physiological studies have confirmed that chloroplast plays a central role in plant cold acclimation and freezing tolerance, which acts as sensor of plants responding to cold and light stress as well as the target of many cold acclimation processes [32]. In the present study, genes related to chloroplast were specifically enriched in the DEGs up-regulated in cold (Additional files 8, 9 and 14), covering almost all aspects of chloroplast. The largest group of these genes was involved in photosynthesis, encoding chlorophyll a-b binding proteins, photosystem II and photosystem I components, cytochrome b6/f complex subunit VIII, thioredoxins, ribulose-1, 5-bisphosphate carboxylase/oxygenase large and small subunits, as well as protochlorophyllide reductase and Mg-chelatase. This result reinforced previous findings in woody or shrub plants with strong cold tolerance, including A. mongolicus [15], rhododendron [33] and Populus euphratica [34], in which several photosynthesis-related genes were up-regulated in cold acclimation. However, it was different from those observed in herbaceous plants like Arabidopsis [35], Brachypodium [9], barley [36] and tomato [37]. In these plants with moderate cold tolerance or cold sensitivity, significant suppression of photosynthesis-related genes was observed under cold stress. An up-regulated transcription of the photosynthesis-related genes might be an important mechanism underlying the maintenance of normal or higher photosynthetic efficiency of A. mongolicus during extreme cold winter.
The second group of the up-regulated genes was annotated to encode early light-induced proteins (ELIPs), 20 kDa chaperonin, chaperone protein DnaJ 11, protease Do-like 8 (DEGP8) and lipoxygenase 4 (LOX 4). ELIPs might protect chloroplasts from light-induced damage or photooxidative damage [38]. Chaperonins or chaperone proteins participate in proper folding of protein substrates, protein disaggregation and degradation [39]. DEGP8 located in the thylakoid lumen [40] is probably involved in the cleavage of photodamaged proteins. LOX 4 is involved in the production of jasmonic acid and functions in biotic stress response [41]. To our knowledge, these genes (except for chaperone protein DnaJ 11) were firstly reported with up-regulation in A.mongolicus under cold stress.
Two other chloroplast-related genes up-regulated in cold stress encode chloroplast omega-3 desaturase and glycerol-3-phosphate acyltransferase. Increasing evidence has indicated that both genes participate in the desaturation of fatty acids in chloroplast lipids and play crucial roles in cold response and tolerance [42,43]. An increased expression of these genes might help to maintain suitable fluidity and stability of chloroplast membranes of A. mongolicus in cold environment. The remaining chloroplast-related genes up-regulated under cold stress may have roles in protein import, ion transport, ATP synthesis, starch degradation, protein synthesis and processing. All these genes could provide structural and physiological mechanisms for the protection of photosynthetic apparatus, and thus for extreme tolerance of A. mongolicus in very cold winter.

Transcriptional regulatory network is involved in the adaptation of A. mongolicus to drought and cold environments
TFs are crucial components in stress signal transduction pathways, which directly control the expression of specific sets of stress-responsive genes [44,45]. Of the A. mongolicus transcriptome, at least 160 TFs in 13 TF families were identified as DEGs during drought-and/or cold treatments, such as AP2/EREBP, NAC, WRKY and bHLH families.
DREBs and ERFs (ethylene responsive factors) are two subfamilies in the AP2/EREBP superfamily [46]. It is wellknown that DREBs play central roles in the regulation of cold-and drought-responsive genes and are important to the related stress tolerance in many plant species [47]. ERFs have been recently found to participate in the response of plants to abiotic stresses, and some ERFs can enhance drought, salt and/or freezing tolerance of their transgenic plants [48,49]. In the present study, 9 DERBs and 8 ERFs were up-regulated in the drought-and coldstressed A. mongolicus seedlings respectively, highlighting the importance of these TFs in the adaptation of this species to drought and cold environments.
The NAC family is another group of TFs involved in the tolerance of plants to abiotic stresses, especially to drought and high salinity [50,51]. Our study identified 8 putative NAC genes up-regulated under drought and cold. Similar results were obtained in the WRKY family, a well-known family involved in biotic-stress tolerance in plants. Recent evidence showed that the WRKY family also functions in response to abiotic stresses [52]. Our findings provide supports to these results, suggesting an important role of the NAC and WRKY families in drought and cold tolerance.
The MYB, bHLH and C2H2 families are also involved in plant abiotic stress tolerance [5,45,53]. In the present study, different members of these families were either up-or down-regulated under drought and cold treatments. This differed obviously from the three TF families discussed earlier which were almost up-regulated during the stress treatments. In summary, sophisticated transcriptional regulation could participate in the adaptation of A. mongolicus to adverse environments, and several families, especially the AP2/EREBP (DREBs and ERFs), NAC and WRKY families might be important in the processes.

Conclusions
A comprehensive transcriptome resource is generated for A. mongolicus and transcriptome profiles under drought and cold stresses are characterized. Comparative transcriptome analysis identified many genes and pathways commonly or specifically regulated in the two stresses, implicating molecular basis of drought and cold tolerance in A. mongolicus. A number of drought-and cold-regulated genes were firstly identified in this study, including some unknown genes. These newly found genes could be important to A. mongolicus in combating extremely cold and drought environments in desert. Our data will facilitate further molecular studies on stress tolerance of A. mongolicus and provide new insights into plant adaptation to harsh environments.

Plant stress treatment and sampling
A. mongolicus seeds were collected from a single shrub grown in the desert area in Dengkou County, Inner mongolia, China. These seeds were surface-sterilized with Clorox solution and soaked in water for three days at 25°C and then sown in pots filled with sand. The pots were placed in a greenhouse at 25°C under long-day condition (16 h light/8 h dark cycle). The seedlings were watered every three days with 1/5 strength of Hoagland's solution. Three months after germination, the seedlings were treated as follows.
For de novo transcriptome sequencing and assembly: the seedlings were placed into a low temperatureprogrammable incubator (BD-PRX-1000A; -10°C to 50°C) for cold treatment at 4°C to -8°C under a 16 h dim light/8 h dark cycle in a gradual cooling way: 4°C for 6 h, 2°C for 14 h, 0°C for 4 h, -2°C for 12 h, -5°C for 8 h, and -8°C for 4 h, in total for 48 h; the cooling rate was 0.5°C per minute. This cold treatment followed by a 3-day recovery period can lead to a slight lethal damage to some leaves of the A. mongolicus seedlings according to our preliminary experiments (unpublished data). The seedlings were carefully removed from sand at 1 h (4°C), 3 h (4°C), 8 h (4°C 6 h + 2°C 2 h), 16 h (4°C 6 h + 2°C 10 h), 24 h (4°C 6 h + 2°C 14 h + 0°C 4 h) and 48 h (cold treatment finished) after the initiation of the cold treatment at 4°C respectively, and washed with pre-cooled water (in a refrigerator at 3°C) and collected. Drought treatment was performed as follows: the seedlings were carefully removed from sand, washed with tap water, placed on filter paper and transferred in the same incubator for dehydration at 25°C under a 16 h dim light/8 h dark cycle. The seedlings were collected at the same time points as cold stress (after 1, 3, 8, 16, 24 and 48 h of dehydration). The untreated seedlings (0 h) were used as control (CK).
For transcriptome profile analysis: the seedlings were placed into the same low temperature-programmable incubator for cold treatment in a gradual cooling way: 4°C for 4 h, 2°C for 12 h, 0°C for 2 h, -2°C for 2 h, -4°C for 2 h, and -6°C for 2 h, in total for 24 h. The drought treatment was performed using dehydration on a filter paper as described earlier. The seedlings were sampled at 2, 8 and 24 h after the initiation of the cold or drought treatment respectively. The untreated seedlings (0 h) were used as control (CK). For each time point, four seedlings were sampled. To minimize the bias of random changes from a single sample, each time point was sampled from three independent stress treatments and the total RNA was pooled equally for further experiment.
Young leaves were sampled from a single A. mongolicus plant grown in the desert area mentioned above in summer (late July) and winter (late December). All samples from stress treatments, controls and the young leaves were immediately frozen in liquid nitrogen and then stored at -76°C for RNA extraction.

RNA extraction and Illumina sequencing
Total RNA was extracted using a modified Trizol method developed by our laboratory (Trizol agent plus 2% β-ME and 3 M KAc with pH 4.8). The quantity of RNAs was verified by an ultraviolet spectrometer (OD 260 /OD 280 ratios with 1.84 to 2.04). RNAs were dissolved in DEPCtreated H 2 O. For de novo transcriptome sequencing and assembly, total RNAs extracted from each sample were pooled equally, and 20 μg of the pooled RNA was used for the isolation of poly (A) + mRNA. For transcriptome profile analysis, poly (A) + mRNA was purified from 10 μg total RNA from each cold-or drought-stressed seedlings and control, respectively.
The mRNAs were isolated by using oligo (dT) magnetic beads and was interrupted into short fragments by adding fragmentation buffer. Taking these short fragments as templates, the first-strand cDNA was synthesized using random hexamer-primer and reverse transcriptase. The second-strand cDNA was synthesized using buffer, dNTPs, RNase H and DNA polymerase I. After the end reparation and ligation of adaptors, the products were amplified by PCR and purified using the QiaQuick PCR extraction kit (Qiagen, Valencia, CA) to create cDNA libraries. The cDNA libraries were sequenced by Beijing Genomics Institute (BGI)-Shenzhen (Shenzhen, China) on an Illumina HiSeq™ 2000 according to the manufacturer's protocols. For de novo transcriptome assembly, the 90 PE (paired-end) strategy was adopted; while for transcriptome profile analysis, the 50 SE (single-end) strategy was used.

Transcriptome assembly
Raw reads from the library were filtered to remove the dirty raw reads if they contain adaptors, over 5% of unknown nucleotides or over 50% of low quality bases (quality value ≤ 5). Clean reads were used for de novo assembly with Trinity, a short read assembling program released recently [54]. Briefly, Trinity firstly combined reads with certain length of overlap to form contigs. These reads were then mapped back to contigs. Finally, Trinity connected the contigs and generated unigenes, i.e. the sequences that cannot be extended on either end.

Functional annotation and classification of the unigenes
The unigenes were aligned to protein databases in the NCBI in the priority order of NR, Swiss-Prot, KEGG and COG by BLASTX algorithm [55], with a cut-off Evalue ≤ 1e-5. Proteins with the highest sequence similarity were retrieved as functional annotations for a given unigene. For unigenes that did not align to any of the above databases, ESTScan software was used to predict their coding regions and corresponding amino acid sequences [56]. To identify putative TFs, all unigenes were blasted against the plantTFDB (http://planttfdb.cbi.edu.cn) with a cut-off E-value ≤ 1e-5, and the top hit for each query was applied. Unigenes with no significant hits in the plantTFDB but with annotations as TFs in the NR and Swiss-Prot databases (E-value ≤ 1e-10) were manually selected and classified into the corresponding TF families. To ascertain the sequence similarity of A. mongolicus to other plant species, total unigenes were searched by BLASTN against the PlantGDB-generated Unique Transcripts (http://www.plantgdb.org/prj/ESTCluster/) (similarity ≥ 80% and E-value ≤ 1e-5).
The GO annotations of the unigenes were obtained by using Blast2GO program [57]. The WEGO software was applied to perform GO functional classification for all annotated unigenes and to show the distribution of gene functions of the species at macro level [58]. After aligning the unigenes to the COG database, the COG functional annotations were acquired. The KEGG pathway annotations were obtained by aligning the unigenes to the KEGG database [59].

Transcriptome profile analyses
The clean reads from each library were mapped to the assembled transcriptome sequences (used as reference genes) using SOAPaligner/SOAP2 [60]. Mismatches with no more than two bases were allowed in the alignment, and ambiguous reads with multi-position match were excluded. Gene expression levels were measured by the RPKM (the numbers of reads per kilobase of exon region in a given gene per million mapped reads) values [61]. If there were two or more transcripts for a gene, the longest one was used to calculate its RPKM. To identify genes responding to drought and cold, the number of reads and RPKM for each coding region in the control and each stress-treated library were determined, and the ratios of RPKM value for each gene in six stresstreated libraries to that in control library were calculated. Statistical significance of the differential expression value for each gene was determined according to Audic and Claverie [62].
Genes were selected as DEGs according to the FDR values ≤ 0.001 and at least a two-fold change (the absolute value of log 2 Ratio ≥ 1) in expression level between two libraries. GO enrichment analysis was performed by mapping all the DEGs to GO terms in the GO database (http://www.geneontology.org/) and the number of unigenes included in each term was calculated. Hypergeometric test [63] was used to identify significantly enriched GO terms in DEGs compared with the assembled transcriptome of A. mongolicus using p-value ≤ 0.05. Pathway enrichment was analyzed by mapping DEGs to the KEGG database (http://www.genome.jp/kegg). The number of unigenes involved in each pathway was calculated. By comparing with the assembled transcriptome of A. mogolicus, the significantly enriched pathways were identified using p-value ≤ 0.05.

Validation of DEGs by RT-qPCR
Eleven genes with different expression patterns were chosen for the validation. The template cDNAs were synthesized from 1.5 μg of total RNAs (pre-treated with DNase I) from the same samples as those for the transcriptome profile analysis using M-MLV Reverse Transcriptase (Promega, USA). The 2 × SYBR-Green I RT-PCR Master Mix (Takara, Japan) was used as labeling agent, and α-tubulin of A. mongolicus was served as an internal reference gene. The reaction mixture (20 μL) contained 10 μL 2× Master Mix, 10 μM forward and reverse primers (Additional file 15) (0.4 μL each) and 1 μL template cDNA (diluted 10 folds with deionized water). The PCR program was as follows: 95°C for 30 sec; the first step: 40 cycles of 95°C for 5 sec, 55 to 60°for 20 sec and 72°C for 20 sec; the second step: 1 cycle of 95°C for 15 sec, 55 to 60°for 1 min, and 72°C for 15 sec. The reaction was performed on a Statagene Mx3000 PCR system. Three independent biological replicates were performed for each sample. The relative expression levels of the selected unigenes were calculated using the relative 2 -ΔΔCT method [64]. Results represent mean standard deviations of the three experimental replicates.