De Novo Transcriptome and Small RNA Analysis of Two Chinese Willow Cultivars Reveals Stress Response Genes in Salix matsudana

Salix matsudana Koidz. is a deciduous, rapidly growing, and drought resistant tree and is one of the most widely distributed and commonly cultivated willow species in China. Currently little transcriptomic and small RNAomic data are available to reveal the genes involve in the stress resistant in S. matsudana. Here, we report the RNA-seq analysis results of both transcriptome and small RNAome data using Illumina deep sequencing of shoot tips from two willow variants(Salix. matsudana and Salix matsudana Koidz. cultivar ‘Tortuosa’). De novo gene assembly was used to generate the consensus transcriptome and small RNAome, which contained 106,403 unique transcripts with an average length of 944 bp and a total length of 100.45 MB, and 166 known miRNAs representing 35 miRNA families. Comparison of transcriptomes and small RNAomes combined with quantitative real-time PCR from the two Salix libraries revealed a total of 292 different expressed genes(DEGs) and 36 different expressed miRNAs (DEMs). Among the DEGs and DEMs, 196 genes and 24 miRNAs were up regulated, 96 genes and 12 miRNA were down regulated in S. matsudana. Functional analysis of DEGs and miRNA targets showed that many genes were involved in stress resistance in S. matsudana. Our global gene expression profiling presents a comprehensive view of the transcriptome and small RNAome which provide valuable information and sequence resources for uncovering the stress response genes in S. matsudana. Moreover the transcriptome and small RNAome data provide a basis for future study of genetic resistance in Salix.


Introduction
The increasing concern about climate change and energy security has resulted in the focusing on the economic importance of Salix species, given their utility for bioenergy production. However, plants have to tolerate many abiotic and biotic stresses, which are serious threats to agriculture and forestry. These common environmental stresses include salinity, extreme temperatures, drought, chemical toxicity, oxidative stress, pests, and pathogen infection. To exploit their potential for renewable energy, willows need to be kept free of pests and diseases and the yield should be improved without significantly increasing the requirement for fertilizers and water [1].
Previous studies have reveled many Salix species have the capacity to tolerate abiotic and biotic stress. The willow species S. caprea and S. cinera are grown on nutrient-deficient and industrially-contaminated soils; thus, they sever as biological filters for waste water and sludge disposal [2]. Nine different clones of six species of Salix (Salix cordata Muhlenb. non Michaux, S. fragilis L., S. caprea L., S. cinerea L., S. burjatica Nazarov. and S. viminalis L.) and one hybrid (S. x calodendron Wimm.) were exposed to heavy metals in solution culture showed reduced phytotoxicity and increased metal resistance [2]. Two Russian willow species, Salix viminalis and Salix dasyclados have shown strong genetic control of freezing resistance at during autumn [3]. Salix caprea, Salix repens and their F1 hybrids had the ability to tolerate different insects and pathogens [4].
Many willow species are characterized by their biological property and ability to adapt to stressful conditions, including high biomass, tolerance to heavy metals, rapid growth, and tolerance of flooding [5]. Salix matsudana Koidz. is a deciduous, rapidly growing tree that is native to northeastern China, and is one of the most widely distributed and commonly cultivated willow species in China [6]. The Chinese name for S. matsudana is 'Han Liu', where 'Liu' means willow and 'Han' means that the willow is drought tolerant. Salix matsudana Koidz. cultivar 'Tortuosa' is widely propagated and planted as an ornamental tree in China because of its attractive twisted stem. This variant is also a water loving plant, which is often cultivated near the river.
The high-throughput next generation sequencing technology such as RNA-seq and Digital Gene Expression (DGE) is an efficient method to discover genes in an unbiased manner [7][8][9][10]. Measurement of mRNA and miRNA expression levels, and elucidation of the regulatory relationship between miRNA and their corresponding mRNA targets are crucial for understanding many pathways and biological systems including plant development, and biotic and abiotic stress responses [11].
However, there have been a few studies probing the ability of S. matsudana withstand drought stress. In this study, we performed a transcriptome and miRNA sequencing from one-year-old stem shoots of Salix matsudana Koidz. and Salix matsudana Koidz. cultivar 'Tortuosa' using deep sequencing technology to discover the differentially expressed genes and miRNAs in stem shoots of the two varieties of Salix. Our results showed that 196 transcripts, which are mainly related to plant abiotic and biotic stress responses, are commonly up-regulated in S. matsudana. These up-regulated genes in S. matsudana included a set of biotic and abiotic stress response genes, such as pathogen defense, insect resistance, antibiotics tolerance, antioxidant, and hormone related genes. Functional analysis of miRNA targets also showed that many genes were involved in stress resistance in S. matsudana.
The transcriptome and small RNAome reported here on the two Chinese willow cultivars significantly enhance our knowledge of stress response genes in Salix.

Plant materials
Shoots from a natural population of Salix matsudana and Salix matsudana Koidz. cultivar 'Tortuosa' were harvested in spring 2013 at Beijing Botanical Garden of China. To ensure that there absolute minimal differences in the growth conditions, all shoots were sampled within the distribution of 30630 m 2 . Samples were frozen in liquid nitrogen immediately for the extraction of total RNA.
This study was carried out in strict accordance with the recommendations in the Guide for Observational and field studies (http://www.plosone.org/static/publication). All necessary permits were obtained for the described field studies. The sampling of all individuals of Salix. matsudana and Salix matsudana Koidz. cultivar 'Tortuosa' were approved by Xiping Zheng, director of Beijing Botanical Garden of China.

Libraries construction and deep sequencing
For each group of samples, equal numbers of shoots from five individual plants were pooled for the total RNA extraction with the TRIZOL reagent according to the manufacturer' instructions (Invitriogen). Small RNA and transcriptome sequencing, used the total RNA from the different pooled samples. After eliminating traces of genomic DNA by treatment with DNase I, an automated capillary gel electrophoresis was used to assess RNA quality and quantity using a Bioanalyzer 2100 with RNA 6000 Nano Labchips (Agilent Technologies Ireland, Dublin, Ireland). sRNA and Transcriptome sequencing libraries were prepared using an Illumina Small RNA Sample Prep Kitand an Illumina TruSeq RNA Sample Prep Kit, respectively. After the two libraries had been prepared, raw reads generated by using the Illumina Hiseq 2000 were initially processed to get clean reads. Then, all the clean reads were assembled using a de novo assembly program Trinity [12,13].
Gene validation and quantitative real-time PCR analysis Ten selected genes specific PCR primers (Table S2) were designed corresponding to the conserved region of willow cDNA sequences from the sequenced database. PCR was performed in a total volume of 50 ml containing 2.5 mmol/L Mg 2+ , 0.15 mmol/L dNTPs, 0.4 mmol/L of each primer, 15 ng cDNA and 0.5 U PFU DNA polymerase (NEB). The PCR products were purified and ligated into the PMD19-T vector (Takara Bio), and then transformed into E. coli Jm109. After sequencing the positive clones with ABI 3730 (Applied Biosystems, USA), the sequences were validated to the corresponding transcriptome seqences.
Quantitative real-time PCR of mature miRNA was performed following the high-stringency protocol in which a poly A polymerase was used to add a poly A tail. The selected genes real-time quantitative PCR were conducted using the SYBR Green Perfect (Takara) and StepOnePlus TM System (Applied Biosystems). All of the PCR products were sequenced and the dissociation curve was analyzed to verify the amplification specificity. The purified PCR products were used to make standard curve for establishing a quantitative correlation between the CT values and the transcript copy numbers [14]. Each qRT-PCR reaction was repeated at least three times, and each standard curve contains at least 5 points. 5.8 s rRNA and Populus actin gene were used for reference gene for miRNA and transcripome genes expression validation, respectively.

Small RNA prediction and identification of new mi RNAs in willow
Small interfering RNA (siRNA) has their own structure feature, which is a 20-25 nt long double-stranded RNA, and each strand of which is 2 nt longer than the other at the 39 end. According to this, reads were aligned with each other in order to find potential siRNA candidate by using the alignment software Bowtie combined with Perl scripts. To identify new miRNA, sRNAs with non-annotation information were subjected to a secondary structure prediction by using the software RNAfold [15]. Further identification of new miRNA was conducted by Mireap (http:// sourceforge.net/projects/mireap/).

Result
Transcriptome sequencing by RNA-seq and de novo assembly Two RNA-seq libraries were prepared from the total RNA extracted from Salix matsudana ('SM') and Salix matsudana Koidz. cultivar 'Tortuosa' ('SMT') and were subjected to pair-end read (PE) sequencing using Illumina platform. We removed reads with adaptors, reads with unknown nucleotides larger than 5%, and low quality reads to obtain clean PE reads in SM and SMT samples. A total of 60,521,772 clean PE reads consisting of 6,506,875,453 nucleotides (nt) with an average GC content of 45.11% were obtained in the SM sample while 58,648,772 clean PE reads consisting of 6,292,916,298 nt with an average GC content of 45.14% were obtained in the SMT sample (Table 1). An overview of the raw reads data is given in Figure S1. All highquality clean reads were assembled into 106,430 contigs with an average length of 944 bp having a minimum length of 201 bp and a maximum length of 17,226 bp ( Figure 1A). The contigs were further joined into 48,817 unigenes with an average length of 773 bp, a N50 length of 1,468 bp and a N90 length of 294 bp ( Figure 1B).  Functional annotation and classification of the assembled unigenes Gene annotation was performed based on sequence homologies to the databases of NCBI non-redundant protein (Nr), NCBI nucleotide sequence (Nt), Kyoto Encyclopedia of Genes and Genomes (KEGG), Swiss-Prot protein, Protein family, Gene Ontology (GO), and Cluster of Orthologous Groups (COGs) using BLAST2GO. In total, 2,592 (5.2% of all unigenes) unigenes annotated in all databases, and 34,049 (69.74% of all unigenes) significantly matched a sequence in at least one of the public databases (Table 2). For the nr annotations, 29,657 of unigenes (60.75%) were found to have matches in the database. Further analysis of the BLAST data showed that 64.11% of the top hits had strong homology with the E-value ,1.0e 250 while 35.89% of the matched sequences showed moderate homology with the Evalue between 1.0e 25 and 1.0e 250 ( Figure 1C).
The identity distribution pattern showed that 79.44% of the sequences had a similarity higher than 80% ( Figure 1D). The majority of the annotated sequences corresponded to the known nucleotide sequences of plant species with 83.14%, 6.46%, 2.87%, 0.63%, and 0.36% matches with Populus trichocarpa, Ricinus communis, Vitis vinifera, Glycine max, and Medicago truncatula respectively.
The size distribution of the BLAST-aligned coding sequence (CDS) and predicted proteins are shown in Figure S2A, B, respectively. The remaining 41.49% (20,256/48,817) of unigenes that did not match sequences in the databases were analyzed by ESTScan to predict coding regions. An additional 12,740 unigenes (62.89%) also showed orientation in the transcriptome coding sequence ( Figure S2C, D). The sequences without a homologous hit may represent novel genes specifically expressed in willow stem shoots, or they could be attributed to other technical or biological biases such as assembly parameters. Furthermore, some cDNA are non-coding, lineage-specific or highly variable and need to be further verified.
GO annotation is an international classification system that can provide standardized vocabulary for assigning functions of the uncharacterized sequences [16]. A total of 23,249 unigenes (47.62% of all the assembled unigenes) were assigned at least one GO term. GO annotation assignments classified 21,577 unique contigs into 22 subcategories of the biological process category, 11 subcategories of the molecular function category, and 14 subcategories of the cellular component category at level 2 ( Figure 2). Among biological processes, transcript sequences assigned to biological regulation, cellular process, metabolic process and single-organism were the most abundant. The subcategories with the most highly abundant transcripts in the cellular components category included cell, cell part, macromolecular complex, membrane, and organelle. Within the molecular function category, the majority of the GO terms were predominantly assigned to the binding and catalytic activity.
KEGG is a database resource that integrates genomic, chemical and systemic functional information that can facilitate to understand the biological functions of genes in terms of networks [18]. The assembled unigenes were annotated with KEGG Orthology (KO) numbers using BLASTx alignments against KEGG with a cut-off E value of 10 25 . Overall, 5,444 of 48,817 (11.15%) unigenes were matched in the database and were assigned to 242 KEGG pathways. Among the five categories of pathway hierarchy 1, the 'metabolism' represented the largest group, which were sorted into 11 subcategories including carbohydrate metabolism, energy metabolism, amino acid metabolism, lipid metabolism, nucleotide metabolism, and some others ( Figure 3B).

Detection of differentially expressed genes between two transcription libraries
A quality control test for the data assembled from each cDNA library confirmed that they were suitable for statistical analysis for differentially expressed genes identification ( Figure S3). A total of 292 differentially expressed genes were detected between SM and SMT groups when log 2 (fold change) .1 and P,0.005 were used as cutoff values. Further analysis showed that 196 genes involved in stress resistance, pathogen defense, antibiotics tolerance, insect resistance, antioxidant, hormone related and ribosomal RNAs were up regulated ( Figure 4A), whereas 96 genes involved in cell wall synthesis, photosynthesis, stress resistance and antibiotics tolerance were down regulated in SM as compared to SMT. Most of the differentially expressed genes are related to plant adaption stress which includes drought, cold, salt, pathogen, insect, and antibiotics. Among the 196 up regulated genes, 133 genes (67.86%) are related to stress adaption indicating that S. matsudana has a better ability to tolerant harsh environment.
Deep sequencing of small RNA Two libraries were constructed for small RNA deep sequencing, and we identified 12,532,326 and 20,113,588 sequence reads respectively from SM and SMT. Among the two data sets, the most abundant classes of small RNA found were 21 and 24 nt, which concurs with most reports of deep sequencing of small RNA in plants ( Figure 4B) [19][20][21]. The sRNA were classified into eight categories (miRNA, tRNA, rRNA, snRNA, snoRNA, ta-siRNA, repeats, and others) by comparing the acquired to miRBase (http://microrna.sanger.ac.uk/sequences), Rfam (http://rfam. sanger.ac.uk/) and Genbank (http://www.ncbi.nih.gov/Gen bank/) ( Figure S4). Besides the unknown sRNA, miRNA made up 16.25% and 15.94% of the total acquired sequences respectively, which were represented in the main by 21 nt small RNAs. There was little variation in the different sRNA categories between SM and SMT (SD = 2.31%) indicating that sRNA categories are relatively stable in the two willow species.
Known miRNA and comparison of miRNA expression levels between the two libraries A total of 166 known miRNAs belonging to 35 families were obtained in SM and SMT. These miRNA exhibited variable abundance. Deep sequencing analysis is a powerful strategy that allows the identification and quantification of differences in the expression of sRNA. In the 35 families, miR159 was the most accumulated miRNA with a total of 410,956 reads detected between the two libraries (87,687 miRNAs in SM and 323,269 in SMT) ( Table 3). There was considerable expression level diversity between the families probably because the willow one-year-old stem is a highly differentiated organ and the genes are in the process of dynamic expression during development. For instance, miR159, miR166, miR319, and miR396 were sequenced more than one hundred thousand times while miRNAs, miR399, miR408, miR481, miR1446, and miR6459 were detected less than 10 times. Members of the same family were also not expressed equally. The two libraries, SM and SMT, had 36 miRNAs that were differentially expressed, among which 24 were up regulated, and 12 were down regulated in SM compared to SMT.

Identification of new miRNAs
Plant miRNA research has significantly progressed because of deep sequencing technology; however, stricter requirements are needed for identifying a miRNA candidate as a reliable new miRNA. Altogether, we obtained 31 predicted novel miRNAs by using two additional criteria [22]. Among these novel miRNAs, 28 and 30 were expressed in SM and SMT, respectively. The majority of the miRNAs (18 of the 31, 58.06%) have a uridine at 59 terminal, showing a preference of AGO1 [23]. The novel miRNAs also displayed unequal expression levels as known miRNAs between the two libraries (Table S1). As none of the novel miRNAs had homologs in miRBase, we presumed that they were S. matsudana-specific, or were restricted to closely related species.

miRNA Targets predicting and sequences validating
Transcriptome of SM and SMT was used to identify miRNA targets. The 35 known miRNA families had 108 affiliated target genes while the 31 novel miRNAs did not target any genes. Functional annotation of the targets was found to be diverse and included the response to stimulus, signal transduction factors, transcription factors, and genes involved in other biological processes. Among the different expressed targets, 24 miRNAstarget pairs showed reverse expression changing pattern when the results from miRNA profiling and the transcriptome profiling were compared ( Figure 5). Functional annotations showed that the 24 miRNA-targets pairs were involved in auxin signaling, stress tolerance, signal transduction factors, and other proteins were involved in various biological processes.
To validate the reliability of RNA-seq, full-length cDNA sequences of twenty selected genes were isolated by T-A cloning and compared with the assembled sequences. The length of these genes ranged from 192 to 1,960 bp, and the sequence variation was minimal. All of the genes had the sequence pairwise identity bigger than 98% (Table S2), which validated the reliability of RNA-seq. To partially validate the miRNA and its target genes expression levels, nine miRNA-target pairs were randomly selected for a real-time qPCR experiment. The results showed that the expression of miRNA target pairs changed trends as compared with the sequencing results ( Figure 6).

Overview of the deep sequencing datasets
Although many studies have reported the transcriptome analysis of willow, here have been no studies on mRNA and miRNA expression profiles in S. matsudana. In this study, the RNA-seq technology was used to generate transcriptome and small RNAome data and examined the global gene expression profiling of shoot tips from S. matsudana. This work demonstrates that RNA-seq is a useful and effective tool for de novo transcriptome and small RNA assembly. Moreover, it facilitates the identification of differentially expressed miRNAs and their target genes between the two willow cultivars, even in the absence of a genomic database. We generated both mRNA and siRNA libraries from the shoot tips of the two willow species to substantially increase the available data on willow mRNAs and sRNAs to obtain genomic resources for further investigating candidate genes especially the stress resistance genes in various metabolic pathways in willow species. Comparison of transcriptome sequence data of the two libraries revealed 292 differentially expressed genes with over 67% involved in stress response. Differential expression of miRNAs also showed many miRNAs and their target genes were stress related. The real-time PCR analysis of the differentially expressed miRNAs and genes further verified the miRNA and the transcript expression level revealed by small RNA and transcriptome comparison from RNA-seq data.

Transcription factors related to stress tolerance in S. matsudana
Transcription factors play a key role in modulating the acclimation response of plants to various internal or external cues. They potentially control downstream gene expression in stress signal transduction pathways through activation and repression of the genes after certain stress treatments. Plant genomes contain a large number of transcription factors (TFs). It has been estimated that the Arabidopsis genome codes for at least 1,533 transcription factors and over 5.9% of its total estimated genes. Among these, about 800 genes encode AP2/EREBP, bZIP, zinc finger proteins, Myb proteins, and AUX/IAA types of transcription factors. Sixteen transcription factors have higher transcription levels in SM compared to SMT, which include Crepeat binding factors (CBF), Myb proteins, TINY proteins, AP2/ EREBP, WRKY type, and other types of zinc finger proteins.
CBF proteins bind to the cold-responsive cis-acting elements that contain the same motif (CCGAC), are encoded by AP2/ EREBP multigene families, and mediate the transcription of genes in response to stress. Overexpression of two Arabidopsis AP2/ EREBP genes, CBF1 and DREB1A, resulted in stronger tolerance to drought, salt, and freezing [24][25][26]. Transgenic Arabidopsis plants overexpressing CBF3 under the control of cauliflower mosaic virus (CaMV) 35S promoter also showed a high tolerance to drought, high-salinity, and freezing stress [25,27,28]. Five CBF and AP2/EREBP genes were highly expressed in SM as compared to SMT. The sequenced CBF and AP2/EREBP gene reads in SM are 2 to 30 times over that in SMT indicating that the SM may have a high tolerance to abiotic stress. WRKY family is defined by a DNA binding domain that contains the strictly conserved amino acid sequence WRKY. Members of WRKY family transcription factors have been implicated in the response to biotic and abiotic stresses. Overexpression of the GhWRKY39-1 gene enhances the resistance to pathogen infection and tolerance to high salt and oxidative stress in transgenic Nicotiana benthamiana [29]. The same results were obtained from the research of Tamarix hispida, overexpression of ThWRKY4 conferred tolerance to salt, oxidative and ABA treatment in transgenic plants [30]. The overexpression of a subgroup IIe WRKY transcription factor CaWRKY27 of Capsicum annuum also positively regulates the resistance of tobacco to the pathogen Ralstonia solanacearum infection [31]. In the present study of the two transcriptome libraries, three WRKY showed higher expression levels in SM as compared to SMT that may enhance the tolerance of SM to biotic and abiotic stress. Plant zinc finger transcription factors play an important role in tolerance to various environmental stresses such as drought, cold, osmotic stress, wounding and mechanical loading. The zinc finger proteins chraracterized by zinc finger domains enable protein interaction with DNA. A large number of studies have shown that zinc finger proteins may not be specific to one particular stress, but may regulate responses to several stresses. In A. thaliana, the zinc finger proteins have been shown to be involved in cold, dehydration, reactive oxygen, and salt stress [32,33]. In Populus euphratica, the overexpression of a zinc finger protein gene PSTZ showed enhanced salt tolerance [34]. In rice, seven zinc finger transcription factors were involved in the response to different abiotic stresses such as cold, drought, and mechanical stress [35]. In Capsicum annuum, a zinc finger transcription factor gene functions as a pathogen induced early-defense gene [36]. Six zinc finger transcription factors have a higher transcription level in SM than SMT indicating that the SM has a stronger ability to tolerate abiotic stresses.
Hormones related genes involved in stress tolerance in S. matsudana  level of various phytohormones to protect the cell from infection when plants were attacked by these diverse pathogens [37,38]. Plants hormones also have the ability to adapt changing environments by mediating a wide range of adaptive responses [39,40]. It is notable that xx hormone-related transcripts have significantly higher expression levels in SM than in SMT, which include eighteen ABA related genes, ten auxin response genes, four JA, and two ET related genes. Under abiotic stresses such as drought, cold, and salinity, ABA levels increased in vegetative tissues, which are adaptive responses essential for the plant survival and productivity [41]. In Arabidopsis, increasing the transcription of a key ABA biosynthesis gene, NCED3, resulted in the increased production of ABA and enhanced tolerance to drought-shock and osmotic stress [42]. Overexpression of AtTRE1 gene that encodes the Arabidopsis trehalase conferred the ability to tolerate drought stress through increased sensitivity toward ABA-dependent stomatal closure [43]. A maize calcium-dependent protein kinase gene, ZmCPK4, positively regulates ABA signaling and enhanced drought stress tolerance in transgenic Arabidopsis [44]. Apart from the vital function for abiotic stress tolerance, ABA plays a modulating role in plant pathogen infections [45][46][47]. Many key components of ABA signal transduction have been identified including protein kinases, RNA processing factors, phosphatases, proteasome components, and chromatin remodeling proteins [48,49]. In S. matsudana, thirteen genes involved in the response to ABA transduction had a higher transcription level than SMT. Moreover, five miRNAs, which target five ABA transduction genes were down regulated as compared to SMT. These genes could increase the ABA levels and enhance the ability of S. matsudana to tolerate abiotic and biotic stresses.
Auxin responsive genes also regulate plant stress response. Overexpression of an auxin responsive gene GH3-8 resulted in stronger pathogen resistance in rice and this resistance was shown to be independent of SA and JA signaling [50]. A rice auxin responsive gene OsGH3.13, encoding IAA-amido synthetase was shown to enhance the expression of genes that confer drought tolerance [51]. Ten genes related to auxin response had higher transcription levels in SM than SMT, and one miRNA had a lower expression level in SM than SMT. These genes could change the auxin levels and confer increased stress tolerance ability upon S. matsudana. Auxin plays a major role of in integrating the activities of multiple phytohormones to control plant growth in response to environmental signals [52]. Except ABA and auxin, JA and ET related genes were also detected that had different transcription in SM and SMT, which in turn could result in variable stress response activity. Four JA related genes and two ET related genes had higher transcript levels in SM than SMT, which could enhance the ability of S. matsudana to tolerance stress.
Antioxidants and detoxification genes related to stress tolerance in S. matsudana Biotic and abiotic stresses such as pathogen infection, salt, drought, oxidation and heat stress are accompanied by the formation of reactive oxygen species (ROS). The ROS such as H 2 O 2 , O 2 , and OHare the main damages of membranes and macromolecules [53]. Plants have developed many antioxidation strategies to scavenge ROS. Enzymes such as thioredoxin, ascorbate peroxidase (APX), superoxide dismutase (SOD), glutathione peroxidase, glutathione transferases, phospholipid-hydroperoxide glutathione peroxidase, galactose oxidase, and glutathione reductase form a part of the main ROS scavenging strategy. The increased transcript levels of these enzyme genes can result in enhanced tolerance to oxidative stress. In potato, overexpression an APX gene resulted in enhanced tolerance to methyl viologeninduced oxidative stress and high temperature in transgenic potato plants [54]. Overexpression of a tobacco glutathione S-transferase with glutathione peroxidase activity in transgenic tobacco enhanced seedling growth under a variety of stress conditions [55]. In Solanum tuberosum, substantial accumulation of a plant thioredoxin named CDSP 32 mRNA and protein were revealed upon oxidative treatments [56]. In S. matsudana, eleven genes involved in response to oxidative stress had higher transcription levels than SMT. These genes encode the abovementioned enzymes including one APX, three glutathione transferases, four thioredoxin, one galactose oxidase, and two NADP reductase. The higher transcription levels of these oxidative stress response genes in SM compared to SMT suggestion that SM has stronger stress tolerance.