Genome Wide Transcriptional Profile Analysis of Vitis amurensis and Vitis vinifera in Response to Cold Stress

Grape is one of the most important fruit crops worldwide. The suitable geographical locations and productivity of grapes are largely limited by temperature. Vitis amurensis is a wild grapevine species with remarkable cold-tolerance, exceeding that of Vitis vinifera, the dominant cultivated species of grapevine. However, the molecular mechanisms that contribute to the enhanced freezing tolerance of V. amurensis remain unknown. Here we used deep sequencing data from restriction endonuclease-generated cDNA fragments to evaluate the whole genome wide modification of transcriptome of V. amurensis under cold treatment. Vitis vinifera cv. Muscat of Hamburg was used as control to help investigate the distinctive features of V. amruensis in responding to cold stress. Approximately 9 million tags were sequenced from non-cold treatment (NCT) and cold treatment (CT) cDNA libraries in each species of grapevine sampled from shoot apices. Alignment of tags into V. vinifera cv. Pinot noir (PN40024) annotated genome identified over 15,000 transcripts in each library in V. amruensis and more than 16,000 in Muscat of Hamburg. Comparative analysis between NCT and CT libraries indicate that V. amurensis has fewer differential expressed genes (DEGs, 1314 transcripts) than Muscat of Hamburg (2307 transcripts) when exposed to cold stress. Common DEGs (408 transcripts) suggest that some genes provide fundamental roles during cold stress in grapes. The most robust DEGs (more than 20-fold change) also demonstrated significant differences between two kinds of grapevine, indicating that cold stress may trigger species specific pathways in V. amurensis. Functional categories of DEGs indicated that the proportion of up-regulated transcripts related to metabolism, transport, signal transduction and transcription were more abundant in V. amurensis. Several highly expressed transcripts that were found uniquely accumulated in V. amurensis are discussed in detail. This subset of unique candidate transcripts may contribute to the excellent cold-hardiness of V. amurensis.


Introduction
Temperature is one of the primary environmental factors that influences and limits the growth and development of plant species. Low and freezing temperatures that occur in far northern and southern latitudes not only limits the suitable geographical locations for growing crops and horticultural plants, but also affects productivity [1]. Generally the cold-tolerant varieties of different species are used in conventional plant breeding as the main resource for increasing the freezing tolerance of cultivars. However, breeding for increased freezing tolerance without knowledge of the molecular processes responsible for cold hardiness traits limits the potential for crop improvement. Studies which examine signal transduction and gene expression changes during cold stress [1][2][3][4][5][6] not only help to reveal the sensing and regulatory mechanisms important for surviving low temperature stresses in plants, but also provide an approach to characterizing candidate genes for genetic improvement of freezing tolerance of agronomic plants [6].
Many plants show increased freezing tolerance following exposure to low nonfreezing temperatures, a phenomenon known as cold acclimation. During acclimation, plants sense environmental cues, e.g. cold signals, and reorganize their transcriptome in response. A series of biochemical and physiological changes then occur to protect plants from freeze-induced injury. While much is unknown about the actual perception of cold, recent studies indicate that plant cell membranes may play a major role in sensing decreasing temperatures. This perception results in transient increases in cytosolic Ca 2+ levels, triggering a cascading biochemical and molecular reactions [3,5,7]. Abscisic acid (ABA) may also function as a secondary signal to transduce, at least in part, cold signal pathways [5]. Transcription factors also respond to cold signals at the early stage during low temperature exposure [8]. Perhaps the best example to date are the C-repeat binding transcription factors (CBF), which were identified as a key responders to low temperature stress in plants under cold stress [9][10][11][12][13][14]. The expression of CBF is up-regulated by another transcription factor, ICE1 (inducer of CBF expression 1) [15][16][17][18][19]. In the CBF-dependent pathway, the CBF protein recognizes the CRT/DRE cis-element in the promoter regions of COR (cold response) genes, which in turn activate transcription of these downstream genes and leads to chilling and freezing tolerance to plants. In a CBF-independent pathway, the transcription factors HOS10 (a R2R3 myeloblastosis type) play pivotal roles in the regulation of cold-responsive genes and freezing tolerance [20]. Although great progress has been made in recent years in understanding cold signal transduction in plant, the differences in genes or expression patterns in sensitive and tolerant species under cold stress still needs further investigation.
Grape is one of the most important and widely grown fruit crops in the world. The majority of cultivated varieties have been derived from one species, Vitis vinifera. Grapevine is cultivated for the production of fresh fruit, dried raisins, pressed for juice, or fermented to produce wine. Due to the popularity of this crop, grapevines are cultivated across broad geographic regions with many different climates. One of the primary limiting environmental conditions for grapevine cultivation is low temperature stress. If given sufficient time to acclimate to decreasing temperatures in the autumn, many V. vinifera cultivars can tolerate temperatures as low as 215uC in midwinter without suffering lethal injury [21,22]. However, if temperatures drop too rapidly in late fall before acclimation occurs, or temperatures decrease below acclimated hardiness, serious injury can occur to buds and root tissues. Depending on the intensity and duration of freezing temperatures, partial or complete loss of fruit production can occur in the following year. In some cases, freezing temperatures can also result in damage severe enough to cause trunk damage and plant loss. Freezing tolerance of V. vinifera also drops quickly when buds break in early spring. As a result, early spring frost can damage floral primordia and reduce crop yields [22,23]. Some wild Vitis species such as the North American species V. riparia Michx. and the Asian species V. amurensis can tolerate midwinter temperatures of 230uC and lower [22]. Thus these wild Vitis species have been used in grapevine breeding programs for the selection of new freezing tolerant cultivars. Most recently, studies of the molecular and genomic differences between these freezing tolerant species and freezing sensitive V. vinifera have been used to begin understanding the gene regulation of increased freezing tolerance. These studies may also provide candidate molecular elements for transgenic based genetic improvement.
Several CBF transcription factors have been identified in Vitis. Initially, three CBF/DREB family members were identified from V. vinifera and V. riparia [24]. Transcripts of VvCBF 1-3 accumulated quickly after low temperature, drought and exogenous ABA treatments. VvCBF4 was also identified from both V. vinifera and V. riparia [25]. The transcriptional level changes of VvCBF4 were seen to be maintained for several days, different from the transient expression of VvCBF1-3 [25]. No significant difference in expression pattern was observed between V. riparia and V. vinifera. [24,25]. Expression level analysis shown that the VvCBF4 was induced after 4 h at 4uC in leaf, stem and flower of V. vinifera cv. Koshu [26]. Transgenic lines of grape with over-expression of CBF and CBF-like genes; VvCBF2, VvCBF4, VvCBFL or VvZFPL (V. vinifera B-box-type zinc finger protein), show improved freezing tolerance without cold acclimation [21,26]. In V. vinifera cv. Cabernet Sauvignon, the changes of transcript abundance in shoot tips when exposed to chilling (5uC) were evaluated by microarray analysis and the results indicate that chilling affects many calcium signaling transcripts [27].
Despite evidence that CBF genes and their transcriptional cascades are important for freezing tolerance, much about transcriptional regulation mechanisms related to cold acclimation are not clear. In this study, V. amurensis collected from Changbai Mountain (Jilin province in China) with high cold-hardiness was used to investigate responses to low temperature at a transcriptional level. A freezing-sensitive V. vinifera cv. Muscat of Hamburg was also examined as a control to evaluate the specific transcriptome modifications that occur in V. amurensis under cold stress. A well-established whole genome transcriptome analysis method, based on deep sequencing [28] and incorporating published grape genome sequences [29], was utilized to reveal the transcriptional changes that occur after cold treatment.

Ethics statement
No specific permissions were required for these locations/ activities. This study did not involve endangered or protected species.

Plants material and treatment
One-year-old self-rooted seedlings of Vitis vinifera cv. Muscat of Hamburg and Vitis amurensis (collected from Changbai Mountain in Jilin province in China) were grown and maintained in the greenhouse of Institute of Botany, Chinese Academy of Sciences, under a 16-h light/8-h dark photoperiod at 26uC in. Three days before cold treatment, seedlings were transferred into a chamber at 24uC under 16-h light at 6:00 am. Cold treatment was started at 9:00 am with constant light. During the first four hours, the temperature dropped 5uC per hour and was held at 4uC for an additional four hours. The seedlings that were used for control were also transferred into growth chambers but without cold treatment. The shoot apex with one well developed leaf was harvested from three independent replicates. RNAs were isolated for digital expression libraries construction and real-time RT-PCR analysis.

Digital expression libraries construction and Solexa sequencing
Total RNA was isolated from collected samples with or without cold treatment, using plant total RNA isolation kit (Tiandz Inc; Beijing, China). The Gene Expression Sample Prep Kit (Illumina Inc; San Diego, CA, USA) was used for sequence tag preparation according to the manufacturer's protocol, which was also well described by Wu et al. [30]. Briefly, mRNA was purified with biotin-Oligo (dT) magnetic bead adsorption from 6 mg total RNA. First strand cDNA were synthesis with oligo (dT) on the bead. After second-strand cDNA synthesis, double strand cDNA was digested with NlaIII endonuclease, producing a bead-bound cDNA fragment containing sequence from the 39-most CATG to the poly-A tail. These cDNA fragments were purified with magnetic bead precipitation and Illumina adapter 1(GEX adapter 1) was added to newly formed 59 sticky end of cDNA fragments. The junction of GEX adapter 1 and CATG site was recognized by MmeI, which cuts 17 bp downstream of the CATG site, producing 17 bp cDNA sequence tags with GEX adapter 1. The 39 fragments were removed with magnetic bead precipitation and Illumina adapter 2 (GEX adapter 2) was ligated to the new 39 end of the cDNA fragment, which represented the tag library.
The cDNA fragments with GEX adapter 1 and 2 were undergoing 15 cycles of linear PCR amplification by Phusion polymerase (Finnzymes, Espoo, Finland). Resulted 85 base fragments were purified by 6% TBE PAGE Gel electrophoresis. After double strand denaturation, the single chain molecules were fixed onto the Solexa Sequencing Chip (flow cell). Each molecule grew into a cluster sequencing template through in situ amplification, which represented a single tag derived from a single transcript. Four color-labeled nucleotides were added during sequencing that performed by the Beijing Genomics Institute (BGI, www.genomics.org.cn) with Illumina HiSeq 2000 System. The data was submitted into the NCBI SRA database (Accession No. SRP018199). The produced 49 bp sequences contain target tags and 39adaptor. Base-calling were performed using the Illumina Pipeline. After purity filtering and initial quality tests, the reads were sorted and counted for the following analysis.

Sequence annotation
''Clean Tags'' were obtained with Fastx-toolkit (http:// hannonlab.cshl.edu/fastx_toolkit/) by trimming adapter sequences and filtering off adaptor-only tags and low-quality tags (containing ambiguous bases). Sequences alignments were carried out with Bowtie 0.12.8 [31] using Genoscope Grape Genome database (http://www.genoscope.cns.fr/externe/ GenomeBrowser/Vitis/). VBI microbial database (http://vmd. vbi.vt.edu/) and BROAD institute database (http://www. broadinstitute.org/scientific-community/data) were used to exclude the contamination of tags from virus according to Wu et al. [30]. All clean tags were annotated based on transcript sequences of grape reference genes and masked grape genome sequences (exclude the repeating sequences). Only tags perfectly matched or 1 nt mismatch were considered for the further annotation. Sequences encoding proteins of known function were manually categorized into broad functional groups using the Munich Information Centre for Protein Sequences classification as guidance.

Identification of differentially expressed genes
Annotated clean tags for each gene were calculated after alignment and then normalized to TPM (tags per million clean tags, [32,33]). The genes that had TPM less than 10 in both library (NCT or CT) were excluded first. The default value (tag number) of genes that not found in one of the library was 1. Differentially expressed genes (DEGs) during cold treatment in two materials were identified based on a rigorous algorithm developed by Audic and Claverie [34]. P value was used to evaluate the authenticity of differential transcript accumulation [30,34]. Bonferroni corrected p value was applied to control the FDR (false discovery rate) in the multiple comparison and analysis during the identification of DEGs [35]. An ''FDR,0.001and the absolute value of log2-Ratio $1'' was set as the threshold to determine the significance of gene expression difference.

Real-time RT-PCR analysis
Samples were prepared and total RNA was isolated using the same method mentioned above. Real-time RT-PCR was carried out on three independent biological replicates each containing three technical replicates. First-strand cDNA was synthesized from RQ1 RNase-Free DNase (Promega, Madison, Wisconsin, USA)treated total RNA using Superscript III reverse transcriptase (Invitrogen, Carlsbad, CA, USA) and diluted 20 fold as template. Specific primer pairs of selected genes were designed using Primer3 (v. 0.4.0, http://frodo.wi.mit.edu/) and shown in Table  S2. Experiments were carried out using FastStart Universal SYBR Green Master (Roche Diagnostics, Mannheim, Germany) with SteopOneplus TM Real-Time PCR system (Applied Biosystems). Data were analyzed using qbasePLUS software (http://www. biogazelle.com/products). Transcript levels were normalized against the average of the grapevine reference genes: the UBC (ubiquitin-conjugating enzyme, EC922622) and GAPDH (glycer-aldehyde 3-phosphate dehydrogenase, CB973647) according to Reid et al. [36]. The fold change in mRNA expression was estimated using threshold cycles, by the DDCT method.

Digital expression libraries construction and tag sequencing
The shoot apices with one well developed leaf from plant materials were collected and subjected to digital expression library construction. In order to activate plant cold stress responses and trigger changes in gene regulation, a total of 8 hours of cold treatment was used, with four hours gradient cooling from 24 uC to 4uC (the temperature dropped 5uC per hour) and another four hours held at 4uC. A total of four digital expression libraries were constructed from non-cold treatment (NCT) and cold treatment (CT) shoot apex of both V. amurensis and V. vinifera cv. Muscat of Hamburg.
For V. amurensis, 9,423,819 and 9,261,386 tags were sequenced from (non-cold treated) NCT and (cold treated) CT libraries, respectively (Table 1). After filtering out tags containing 'N' and adaptor sequences, tags were align to the VBI microbial database (http://vmd.vbi.vt.edu/) and BROAD institute database (http:// www.broadinstitute.org/scientific-community/data) to exclude any contamination from viruses according to Wu et al. (2010). A total of 9,411,578 and 9,204,750 clean tags were collected from VaNCT and VaCT libraries. Single-copy tags in each library (291,588 in VaNCT and 264,143 in VaCT library) were excluded from further analysis [30]. Finally a total of 9,019,990 and 8,940,607 clean tags were clustered into 202,174 (VaNCT) and 205,283 (VaCT) unique tags. The distribution of unique tag with different copy number in each library were counted and shown in table 1.
For V. vinifera cv. Muscat of Hamburg, 9,590,554 and 9,915,588 tags were generated from NCT and CT libraries respectively (shown in Table 1). After filtering out low quality tags, contamination from viruses, and exclusion of single-copy tags, a total of 9,485,348 and 9,660,003 clean tags were clustered into 182,034 (NCT) and 255,879 (CT) unique tags. The higher number of unique tags observed in the VvCT library may indicate candidate genes related to cold signal response in Muscat of Hamburg. While the difference in unique tag number between cold treated and non-cold treated libraries (73,845 tags) was large in Muscat of Hamburg, the NCT and CT libraries in V. amurensis were much more similar and we identified only 3109 greater unique tags in the VaCT library (Table 1) However, this result is partially explained by the large number of unique tags with copy number between 2 and 5 (127,105 unique tags) in the VvCT library compared with the VvNCT library. These tags are the primary contribution to the significant increase in unique tags observed in Muscat of Hamburg during cold stress.

The saturation of tags in each library
In order to increase the reliability of the expression analysis, the saturation of tags in each library was evaluated by the number of identified genes. The number of tags reached saturation when no new genes are detected. The results are shown in Figure S1. In the four libraries examined, the numbers of identified genes declined rapidly as the number of sequenced tags increased. All the libraries reached a plateau with 6 M tags. No new genes were identified as the tag number closed to 8 M. Considering that more than 9 M of the available tags were generated in each library, the tags were sequenced to saturation and provide a adequate information for the expression pattern analysis.

Alignment of the unique tags to the reference genome
The unique tags from each library were aligned to the published genome and compared with annotated genes from V. vinifera cv. Pinot Noir PN40024 [29] (Table 2). Unique tags that matched the reference or had a single mismatch were counted. For unique tags that only matched one gene, 39.40% of the unique tags in the VaNCT library and 43.61% of unique tags in the VaCT library were identified, while both libraries of Muscat of Hamburg shown higher percentage of aligned unique tags (55.13% in the VvNCT library and 54.68% in the VvCT library). These sets of data were used for the following gene expression pattern analysis. When compared with the total annotated genes in the 126 grape genome, both libraries in V. amurensis had similar identification results (15,549 genes in the VaNCT and 15,604 genes in the VaCT), while the VvCT library had slightly more matched genes (17,126,70.34%) than the VvNCT library (16,188, 66.49%).

Identification of differential expressed genes (DEGs) during cold treatment
Unique tags that perfectly matched reference genes in each library were normalized to tags per million clean tags (TPM) and used to evaluate the expression level of transcripts. To increase the robustness of the data, only the genes that had TPM more than 10 at least in one of the library (NCT or CT library) in each genotype were considered further. The transcripts with at least a two-fold difference in expression during cold treatment in the two genotypes (FDR ,0.001) are displayed in Figure 1. Genes with less than two fold changes between NCT and CT libraries in both materials were excluded from further analysis. The details of DEGs, including original TPM, fold-change, annotation, P value and FDR in both materials are shown in Table S1.
The number of differentially expressed genes with at least a 2 or 5-fold change are shown in Table 3. The results indicate that V. amurensis has a lower number of DEGs during cold stress than Muscat of Hamburg. For V. amurensis, 1,314 genes changed at least 2-fold in the VaCT library, including 893 up-regulated genes and 421 down-regulated genes. In Muscat of Hamburg, 2,307 genes changed at least 2 fold in the VvCT library with 1,333 genes increasing and 974 genes decreasing in expression. A total of 167 and 77 genes were found to be up or down-regulated by at least five fold in the VaCT library, while the numbers increased to 174 and 149 in the the VvCT library. When comparing the DEG (at least 2 fold) between the two grape varieties, 429 genes showed the same expression pattern.
The DEGs with 20-fold or greater expression changes in both grapes are shown in Tables 4 and 5. Although less overall DEGs were found in V. amurensis than in Muscut of Hamburg, more genes (66) were up or down-regulated by more than 20-fold during cold treatment in V. amurensis (Table 4). Fifty-one genes were expressed at greater levels in the VaCT library. These genes could be grouped by gene ontology and were mainly associated with signal transduction (6), transcription (5), transport (4), protein folding (4), metabolism (7) and photosystem (6). Five transcription related genes belong to ERF, MYB and WRKY families, respectively.
Cold regulated protein 27 (COR27, GSVIVT01009490001) was the highest up-regulated gene (332fold) in the VaCT library. Fifteen DEGs were down-regulated greater than 20-fold in the VaCT library. Of these genes, 9 were related to protein folding (Shown in Table 5). GSVIVT01016426001 showed greatest decrease in expression (246-fold) in the VaCT library but has no annotation in the NCBI database. Table 5 also shows other DEGs with large-fold changes (20-fold or more) in the VvCT library. Twenty-three genes were upregulated in the VvCT library and these genes were associated with stress (2), transcription (3), transport (2), photosystem (2), cell cycle (2) and metabolism (4), etc. Three transcription factors belong to the MYB, bZIP and WRKY families, respectively. Twenty-four genes were expressed at greater than 20-fold lower levels in the VvCT library. Similar to the gene ontology observed in the VaCT library, most of the genes were related to protein folding (e.g. heat shock protein and chaperones).
Sixteen DEGs were identified that had 20-fold or greater changes in gene expression and were observed in both V. amurensis and V. vinifera (Table 5

Functional classification of cold stress-related DEGs
The functional classification of DEGs was further examined in both grape varieties to investigate the pattern of transcriptome regulation that occurs during cold stress. The identified DEGs were annotated using BLASTN and BLASTX searches. Genes matching characterized proteins or proteins with putative functions were grouped according to functional categories ( Figure 2). For up-regulated DEGs (Figure 2, right-hand side), genes encoding proteins involved in metabolism comprised the largest functional group in both grape varieties. In V. amurensis, transcription-related DEGs comprised the second largest category. A considerable proportion of the DEGs found in V. amurensis belong o signal transduction and transport categories. Compared with Muscat of Hamburg, a higher proportion of up-regulated DEGs related to metabolism, transcription, signal transduction and transport were found in V. amurensis. Translation and RNA processing related DEGs represented a lower proportion of DEGs in V. amurensis. The second largest category in Muscat of Hamburg was translation related genes (including ribosomal proteins). Genes   Figure 2 (Left-hand side). Many genes responsible for metabolism, transcription, signal transduction, transport and protein folding (including heat shock proteins and chaperones) were significantly down-regulated in both grape varieties. The proportion of genes with reduced expression that classified into translation, cytoskeleton, stress related, protein folding and histone related was higher in the VaCT than in the VvCT library.

Validate the DEGs by real-time RT-PCR analysis
To validate the data from deep sequencing, thirteen upregulated genes with a differential change in expression of 2 to 103 fold were randomly selected from V. amurensis DEGs for real-time RT-PCR analysis. The primers of selected genes are listed in Table S2. UBC and GAPDH were used as reference genes for data normalization according to Reid et al. [36]. The Real-time RT-PCR results in V. amurensis and Muscat of Hamburg are shown in Figure 3. The expression patterns of all detected genes show the same trend using RT-PCR and the Solexa-sequencing method. The list of the genes and the comparison of fold changes between deep sequencing and Real-time RT-PCR in two kinds of materials were shown in Table 7. The scale of the fold change of these genes based on Real-time RT-PCR was dramatically smaller (1-16 fold) than those of deep sequencing based method (2-103 fold, Table 7).

Discussion
V. amurensis shown less sensitive to cold stress than Muscat of Hamburg at transcriptome level In this study we used the Illumina HiSeq 2000 system, a high throughput DNA sequencing platform, to investigate the modification of the transcriptome under cold stress conditions in two grape varieties. Initially, we assumed that V. amurensis, one of the most cold-hardy wild grape species, may quickly response to cold stress and alter the regulation of its transcriptome intensively to overcome the potential for cold damage. In contrast we hypothesized that the cold sensitive variety, Muscat of Hamburg, may show less transcriptional level changes. However, deep sequencing results from cDNA libraries prepared from shoot apices with exposure to a total of 8 hours cold treatment did not support this hypothesis. While many differentially expressed genes were identified in the cold hardy V. amurensis in response to cold stress, a greater number of genes (2,307) were seen to increase in the cold sensitive Muscat of Hamburg.
One potential explanation for the lower number of DEGs identified in V. amurensis is a lack of alignment and annotation of V. amurensis genes to the reference genome. As one of the wild Vitis species, the phylogenetic distance between V. amurensis and PN40024 is larger than V. vinifera cv. Muscat of Hamburg and PN40024. Indeed, higher percentages of the unique tags were found to be aligned to reference genes in Muscat of Hamburg (55.13% in VvNCT and 54.68% in VvCT) than in V. amurensis (39.40% in VaNCT and 43.61% in VaCT library). The similar proportion of matched tags in V. amerensis were also found in previously published research [30], when transcriptome analysis was examined for V. amerensis before and after downy mildew infection ((35.47% in infected library and 36.99% in control). Because of the better alignment results, more genes were identified in Muscat of Hamburg than in V. amurensis (Table 2). In order to confirm that the reduced alignment ability of V. amurensis is responsible for the differences in DEG number between the two kinds of grapes, we carefully checked the DEGs from Muscat of Hamburg in two V. amurensis libraries. The results show that almost all of the DEGs (2,271 from 2,307) in Muscat of Hamburg can be found in at least one of the library in V. amurensis. For V. amurensis, 1,309 DEGs out of 1,314 DEGs can be detected in at least one of the libraries in Muscat of Hamburg. These data indicate that the differences in DEG number between the two grape varieties are unlikely to be due to differences in the proportion of aligned tags. Instead, it suggests that there are very different modifications occurring in the transcriptome of these varieties during cold stress.  V. amurensis responds to cold signal in a different way from Muscat of Hamburg The differences observed in the transcriptome data, including DEG number, gene annotation and functional category, indicate that V. amurensis responds to cold signal in a different way from Muscat of Hamburg. For detected DEGs in V. amurensis and Muscat of Hamburg, very few genes (428 or 87 genes when twoor five-fold was used as threshold) show the same expression pattern during cold treatment (Table 3 and Table S1).
A considerable proportion of DEGs were up-regulated in V. amurensis compared with Muscat of Hamburg under cold stress treatments, suggesting that different genes are responding in the cold tolerant variety. For example, 68% of DEGs were found upregulated in the VaCT library while 58% of DEGs were upregulated in Muscat of Hamburg. Although less overall DEGs were identified in V. amurensis, the functional category analysis indicated that the VaCT library contained higher percentage of up-regulated DEGS involving in several biological processes including metabolism, transcription, signal transduction and transport ( Figure 2). Meanwhile, a lower proportion of downregulated DEGs in these four functional categories were observed in V. amurensis relative to Muscat of Hamburg. It is interesting that a large proportion of DEGs related to translation were identified in the VvCT library (Table 3), implying that Muscat of Hamburg may utilize the synthesis of new proteins as a mechanism for overcoming the damage of cold stress.
The genes which showed the most dramatic difference in expression in both CT libraries also show less similarity between the two genotypes ( Table 4 and 5). For commonly changed genes, COR27 (GSVIVT01009490001) is a 27-kDa protein that responses to cold stress but no functional annotation in Arabidopsis is currently available [37]. Further studies indicated that COR27 is regulated by the circadian clock at warm growth temperatures and cold-induction of COR27 is gated by the clock [37]. Other genes like ACR toxin-sensitivity inducing protein (GSVIVT01025812001), myb TF (GSVIVT01024916001), nitrate transporter (GSVIVT01007624001) and chaperone protein dnaJ (GSVIVT01026228001), are also up-regulated in both kinds of grape. Most of these genes are involved in metabolism, transcription, transport, signal transduction and protein folding, and likely represent fundamental genes that respond to cold signals in Vitis.
Several of the genes which changed in expression in both grape varieties were seen to be expressed at much higher levels in V. amurensis than Muscat of Hamburg. GSVIVT01011652001 was annotated as circadian clock-associated FKF1, which may regulate the day-night rhythm during cold stress in V. amurensis and promote the expression of COR27. GSVIVT01013935001 and GSVIVT01013913001 are similar to ethylene-responsive TF, a type of protein related to CBF genes, that has been seen to be expressed in response to abiotic stress in plants [5,38,39]. GSVIVT01015952001 and GSVIVT01012682001 are annotated as WRKY TF. Genes in this family are related to several plant process including germination, senescence and responses to abiotic stresses such as drought and cold [40,41]. Stilbene synthase 1(STS1, GSVIVT01010589001) was also up-regulated in the VaCT library. STS is widely studied in grape and it is the key enzyme that is responsible for the synthesis of resveratrol in biotic and abiotic stresses [42,43]. It is interesting to observe that this important secondary chemical is expressed in response to cold stress in V. amurensis, potentially linking the biotic defense and antioxidant function of resveratrol, with damage to plant cells from low temperatures. GSVIVT01029173001 is homologous to xyloglucan endotransglucosylase, a gene related to cell wall  modification and has also been observed to play a role in cold temperature response in Arabidopsis [44]. Beta-amylase (GSVIVT01013272001) was also seen to respond to temperature shock and lead to the accumulation of maltose [45].

V. amurensis specific responses
In addition to these genes, a subset of highly expressed genes was only found only in V. amurensis (Table S1) and the majority of the genes relate to the cold stress response. These up-regulated genes are excellent candidate genes for cold stress response and may contribute to the higher cold-hardness of V. amurensis during winter.
Metabolism. The expression of 201 transcripts changed by more than 2-fold and were only observed in the VaCT library. Of these, nine transcripts show at least 5-fold up-regulation and contain more than 20 TPM in at least one of the libraries. Three transcripts were found to be related to starch synthase and degradation, including soluble starch synthase (GSVIVT01004632001), granule-bound starch synthase (GSVIVT01019680001) and alpha-amylase (GSVIVT01020069001). Starch syntheses are a group of important enzymes involved in the synthesis of amylose and amylopectin [46]. Normally, starch synthesis is under strong inhibition even under moderate water deficit condition [47] and none of members in this gene family were up-regulated in the VvCT library (Table  S1). Amylase is responsible for hydrolyzing starch and glycogen into glucose and maltose. The increase of starch-degrading  enzymes during stress condition, which leads to accumulation of soluble sugars, was also found in Arabidopsis [48]. Two transcripts were found to be related to sugar metabolism, including neutral invertase (GSVIVT01034944001) and glucose-6-phosphate 1dehydrogenase (GSVIVT01000913001). In Arabidopsis, both mitochondrial and cytosolic invertase can generate glucose as a substrate for mitochondria-associated hexokinase, contributing to mitochondrial reactive oxygen species homeostasis [49]. Glucose-6-phosphate 1-dehydrogenase (G6PDH) is a key enzyme that regulates the flux of carbon through the pentosephosphate pathway. The transcript of G6PDH was accumulated during salt stress [50] and demonstrates a link between salt and cold stress  response. Other robust up-regulated transcripts were annotated as UDP glycosyltransferase (GSVIVT01028812001), glutathione stransferase (GSVIVT01024290001), carboxylesterase 1 (GSVIVT01010672001) and ATPase subunit 8 (GSVIVT01004977001). UDP glycosyltransferase mediates the transfer of glycosyl residues from activated nucleotide sugars to acceptor molecules [51]. And UDP-glucosyltransferase UGT74E2 in Arabidopsis mediates abiotic stress responses and stress-induced morphological adaptations by regulating auxin homeostasis [52]. Glutathione s-transferase not only functions in detoxification, but also responds to numerous stresses, including those arising from pathogen attack, oxidative stress, and heavy-metal toxicity in plant [53]. Carboxylesterase is an enzyme that is capable of hydrolyzing a wide variety of carboxylic acid esters. Although the functional details are still unknown, a member of this gene family in Arabidopsis (AT1G474800) was detected following cold treatment [54]. ATPase subunit is located in the membrane of mitochondria and constitutes the F0 complex of mitochondrial F-ATPases. The accumulation of ATPase subunit 8 in the VaCT library may indicate its potential roles in mitochondrial membrane stability and cold signal transduction.
Signal transduction. Nine genes that can be classified into signal transduction category were found more abundant (at least 5fold up-regulated) in the VaCT library. Three of them were annotated as serine/threonine protein kinases (GSVIVT01024709001, GSVIVT01035226001 and GSVIVT01015198001). In plant cells, these kinds of kinases accept the information from receptors and convert it into appropriate outputs to regulate several biological processes [55]. A SNF1-type serine/threonine protein kinase of wheat was found to enhance multi-stress tolerance in Arabidopsis [56]. GSVIVT01024979001 and GSVIVT01021283001 were similar to receptor-like protein kinase (RPK), which localize in the cell wall and respond to the external challenges from environment [57]. RPKs that can response to cold stress was also identified in Arabidopsis [58]. GSVIVT01025028001 represents a ras-related protein, which belongs to small GTPase superfamily that participates in vesicular transport [59]. Mitogen-activated protein kinase-kinase-kinase (GSVIVT01022117001), part of the MAPK cascades involved in responses to various biotic and abiotic stresses, hormones, cell division and developmental processes in plant [60], was also up-regulated in the VaCT library. Together with OBP3-responsive protein 1 (GSVIVT01038710001) and choline/ethanolamine kinase (GSVIVT01001426001), these genes may take part in the cold signal transduction in V. amurensis.
Transcription. The DEGs in the transcription-related category were primarily transcription factors. For the 18 TFs that were up-regulated over 5-fold in V. amurensis, all of them had similar expression patterns in Muscat of Hamburg but with lower fold changes. These TFs may play fundamental roles for modifying the transcriptome during plant expose to cold environment. These TFs include ethylene-responsive transcription factor (GSVIVT01013935001, GSVIVT01013913001 and GSVIVT01013931001), WRKY transcription factor (GSVIVT01015952001, GSVIVT01021252001, GSVIVT01024624001 and GSVIVT01012682001), MYB family transcription factor (GSVIVT01024916001), AP2/ERF and B3 domain-containing transcription repressor TEM1-like (GSVIVT01011947001), DRE transcription factor 1 (GSVIVT01002262001), NAC transcription factor (GSVIVT01022354001), GRF domain class transcription factor (GSVIVT01038629001) and others. Members of TFs in these mentioned families were previously identified in response to cold signal and could regulate the expression of downstream genes in plant [8]. V. amurensis also triggers the up-regulated expression of distinct members in identified TF families including WRKY, AP2/ERF, MYB, Homeobox, NAC, et al. (Table 6). For example, six members of WRKY TFs (GSVIVT01011472001, GSVIVT01020136001, GSVIVT01030174001, GSVIVT01001332001, GSVIVT01028147001 and GSVIVT01034968001) were only up-regulated in the VaCT library (Table 6 and Table S1). WRKY TFs often act as repressors as well as activators for transcription, and members of the family play roles in abiotic stresses such as drought and cold [40]. Four AP2/ERF family genes (GSVIVT01016352001, GSVIVT01025100001, GSVIVT01020584001 and GSVIVT01035098001), were also identified with similar expression patterns. In Arabidopsis, CBF1 encodes an AP2 domaincontaining transcriptional activator that binds to the C-repeat/ DRE in response to low temperature and water deficit [10][11][12]. Other members of this gene family are also involved in abiotic stress responses [38]. Further investigation of the V. amurensisspecific TF in AP2/ERF and other TF family members will help to reveal the aspects of transcription that enable this species to tolerate cold stress.
Transport. Two transcripts related to transport were more abundant in the VaCT library, including a nitrate transporter (GSVIVT01008073001) and a hexose transporter (GSVIVT01003181001). The expression of a nitrate transporter gene NRT2 (At1g08090) was induced by cold, mannitol and NaCl treatment in the root of Arabidopsis [61]. Species specific expression differences of nitrate transporter were observed in citrus upon exposure to 4uC [62]. Very high transcript levels of nitrate transporter (431 fold increase) were observed in Poncirus, a cold hardy rootstock for citrus. While much lower levels were observed in Satsuma mandarin [62]. Higher transcription level of nitrate transporter gene indicates the possible role in nitrate acquisition in cold-hardness V. amurensis. A hexose transporter, VvHT5, was highly induced by powdery, downy mildew infection and wounding, suggesting its generalized response to stress in grape [63]. The key regulatory ABA biosynthetic gene, VvNCED1, was activated under these same stress conditions. The promoter region of VvHT5 contains multiple abscisic acid (ABA) response elements, suggesting a role for ABA in regulate the transcription of VvHT5 under stress conditions [63]. In our data, 6-fold up-regulation of hexose transporter (GSVIVT01003181001) was observed in the VaCT library with no obvious change in the VvCT library. And the expression of 9cis-epoxycarotenoid dioxygenase 1 (NCED1, GSVIVT01038080001) was also up-regulated in both grapes (2.9-fold in VaCT and 2.2-fold in VvCT library, Table S1). Future investigation on the role of GSVIVT01003181001 and its regulation by ABA will help to understand the function of sugar metabolism related genes under cold stress.

Stress related
genes. Two transcripts, GSVIVT01025994001 and GSVIVT01021103001, which were annotated as TMV (tobacco mosaic virus) resistance protein N-like and DEHYDRATION-INDUCED 19-like, respectively, were upregulated only in V. amurensis. TMV resistance protein N encodes a TIR-NBS-LRR class protein that confers resistance to TMV [64]. The transcripts of Dehydration-induced 19 in Arabidopsis are regulated by progressive drought stress in an ABA-independent manner [65]. Drought related protein may also respond to cold treatment as a result of the cross-talk between two stress signaling pathways in plants [66]. In fact, four DEGs, which related to drought stress (including AFG16868: dehydrin 1; BAD19956: drought-induced protein RDI; AAO33767: drought-induced protein and NP_177104: early-responsive to dehydration stress protein ERD4), were also found only in the VaCT library with 2 to 5-fold increases during cold treatment (Table S1). Other biotic and abiotic stresses related genes were also changed in VaCT, including disease resistance protein (GSVIVT01019041001,  GSVIVT01005967001,  GSVIVT01001043001  and  GSVIVT01029297001),  salt  induced  protein  (GSVIVT01001340001,  GSVIVT01010794001  and  GSVIVT01029107001) and wound-responsive protein (GSVIVT01019074001 and GSVIVT01009760001). These genes may not only participate in the cold tolerance response but also play roles during the cross-talk between several stress-related networks in V. amurensis.

Conclusion
The reorganization of the transcriptome during cold treatment in V. amurensis and V. vinifera cv. Muscat of Hamburg was evaluated by deep sequencing of short cDNA fragments. Overall, V. amurensis was seen to have less regulation of transcription profiles under cold stress than that of Muscat of Hamburg. However, the proportion of up-regulated DEGs involved in metabolism, transport, signal transduction and transcription were more abundant in V. amurensis. Additionally, several robust changes in transcript level were observed for candidate genes in V. amurensis. These genes may play an important role in the enhanced cold hardiness of V. amurensis and should be the focus of future studies in grapevine.