Potato (Solanum tuberosum L.), which is the third most important food crop after rice and wheat in terms of human consumption, is cultivated throughout the world in temperate and continental areas including North China, North America, and the highlands of South America (https://cipotato.org/crops/potato). This plant is generally considered to be rather sensitive to drought stress due to its sparse and shallow root system of about 50–100 cm depth (Dalla Costa et al., 1997; Vos and Groenwold, 1986). Varieties with diverse traits have been developed through traditional breeding, most of which incorporate traits such as, disease resistance, increasing productivity and quality, and the control of maturity. Although developing new cultivars exhibiting drought tolerance is a key aim of potato breeders, the development of drought-tolerant cultivars by conventional breeding strategies is a very complicated process which takes an extended period of time. Traditionally, potato breeding for crop improvement has relied on phenotypic selection with a limited emphasis on genotypic isolation (Hirsch et al., 2014).

The production of mutants, either naturally or artificially, has been considered a very useful method to develop cultivars with new desired traits within defined germplasm pools (Pathirana, 2012; Shu et al., 2012). In potatoes, which are difficult to breed by conventional methods, such mutant production is an indispensable part of gene discovery for molecular breeding as well as breeding of new cultivars (Fischer et al., 2008; Muth et al., 2008). Treatment of ethyl methanesulfonate (EMS), an alkylating chemical mutagen, has successfully assisted the development of new cultivars in both seed and vegetatively propagated plants (Parry et al., 2009). This mutagen is easy to use without specialized equipment and can provide a very high random point mutation frequency (Sikora et al., 2011). Many cases of the development of varieties by EMS mutagenesis have been reported, but potatoes are not relatively common (Behera et al., 2012; Jabeen and Mirza, 2004; Luan et al., 2007; Taheri et al., 2017).

Next-generation sequencing (NGS) platforms that are relatively rapid and cost effective, such as Roche/454, ABI SOLiD, Illumina/Solexa, and the Helicos Genetic Analysis System, have made it possible to complete functional genomics studies to improve crop genetics at the whole-genome level (Cloonan et al., 2008; Vera et al., 2008; Wang et al., 2009). This NGS technology has been used for RNA sequencing and de novo transcriptome assembly of non-model organisms with or without a reference genome (Haas and Zody, 2010; Li and Dewey, 2011), particularly RNA sequencing of minor crops such as sorghum (Sorghum bicolor) (Mizuno et al., 2012), sunflower (Helianthus annuus) (Livaja et al., 2013), Jerusalem artichoke (Helianthus tuberosus) (Jung et al., 2014), red clover (Trifolium pretense) (Yates et al., 2014), and Camelina sativa (Mudalkar et al., 2014) have been reported. In potatoes, the whole-genome has been sequenced via the Potato Genome Sequence Consortium (PGSC) using the S. tuberosum group Phureja clone DM1-3 516 R 44 (doubled monoploid) and the S. tuberosum group Tuberosum RH89-039-16 (heterozygous diploid), released in 2011 (Xu et al., 2011). This information is provided to biologists, breeders, and geneticists as reference genomes (www.solgenomics.org). Massa et al (2011; 2013) has reported a reference for the potato transcriptome using the reference accession (Massa et al., 2011; 2013), and since then, many transcriptomes corresponding to specific conditions have been published (Gálvez et al., 2016; Gong et al., 2015). Nevertheless, since most potato cultivars are autotetraploid, there is a limit to utilize the completed genome sequencing for breeding.

In the current study, we developed EMS-mutagenized potatoes that showed significant tolerance to drought stress compared to the wild-type (WT) ‘Desiree’ cultivar. The RNA-seq approach was applied to generate time-course transcript expression profiles 48 h after PEG treatment. The sequences are analyzed by both de novo assembly of transcripts and alignment to the published diploid potato genome to identify differentially expressed genes (DEGs) related to drought tolerance at the initial stage of stress. In total, 304,976,315 high-quality reads were retained and 61,100 transcripts and 5,118 DEGs showing up- or down-regulation under drought conditions were identified. Those DEGs were then annotated and assigned putative functions, classifications or pathways by alignment with public databases. We considered this approach to be very meaningful for gene expression studies on the initial response of potato plants under drought stress, using this tolerant mutant. The comprehensive transcriptome data obtained from this study may be useful for further characterizing drought stress responses in potatoes and may enable the breeding of new drought-tolerant potato cultivars.

Generation of potato mutants

Nodal cuttings with one leaf node were harvested from 4-week-old in vitro ‘Desiree’ plantlets, a representative tetraploid commercial potato cultivar, and cultured on a Murashige and Skoog (MS; (Murashige and Skoog, 1962)) basal medium containing 9% sucrose and 8 g l⁻¹ agar in a growth room (17 ± 1°C, 24 h dark) to produce microtubers, which were harvested after 3 months. Media used in this study were autoclaved at 121°C for 20 min and the pH was adjusted to 6.0 before autoclaving. All macro and microelements for tissue culture media, sucrose, and agar were purchased at Duchefa (Haaelem, Netherlands), and other chemicals were purchased at Sigma-Aldrich (USA). To induce mutations, 1,000 microtubers were submerged in a freshly prepared 0.8% (v/v) EMS for 18 h at room temperature under a vacuum hood (following the safety instructions of the manufacturer’s Material Safety Data Sheet), then rinsed in water and kept at room temperature for germination.

Characterization of potato mutants

Mutagenized microtubers were grown to plants in pots (30 × 30 × 21 cm³) containing bed soil in the greenhouse. The water content of each pot was measured three times a week, and the water lost was supplemented to keep equivalent levels according to treatment requirements. After six weeks, plant height and phenotypes were measured before water was withheld from the plants for four weeks to screen for drought tolerance. Plants were re-watered, and tubers were harvested at 100 days total cultivation. WT plants and tubers were cultivated and harvested at normal conditions as a negative control. Germinated sprouts from M1 tubers were surface sterilized with 10% NaOCl and transferred into in vitro culture systems on a MS basal medium. One-week-old in vitro shoots were treated with −1.8 Mpa polyethylene glycol (PEG)-8000 in MS medium for two days. All experiments were performed in three replicates.

Leaf and root samples were fixed in 2.5% paraformaldehyde-glutaraldehyde mixture buffered with 0.1 M phosphate (pH 7.2) for 2 h, post-fixed in 1% osmium tetroxide in the same buffer for 1 h, dehydrated in graded ethanol, and substituted by isoamyl acetate. They were then dried at the critical point in CO₂. Finally, the samples were sputtered with gold in a sputter coater (SC502, Polaron) and observed using a scanning electron microscope, FEI Quanta 250 FEG installed in KRIBB.

Total RNA isolation

PEG-treated in vitro DR and WT plants were sampled at 0, 6, 12, 24, and 48 h post treatment for RNA sequencing. Total RNA was extracted from samples using the TRIzol Reagent (Invitrogen, USA), and were then treated with DNase I (Fermentas, USA) according to the manufacturer’s instructions. Three biological replicates for each treatment were used to prepare pooled total RNA samples. The RNA quality was assessed according to optical density ratios (260 nm:280 nm and 260 nm:230 nm) determined using a NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific, USA).

cDNA library construction and transcriptome sequencing

We prepared RNA-Seq paired-end libraries for the pooled total RNA of each treatment using the Illumina TruSeq RNA Sample Preparation Kit v2 (Illumina, USA). The mRNA was purified using poly (A) selection. The RNA was chemically fragmented and synthesized into single-stranded cDNA using random hexamer primers. The second-strand cDNA was synthesized to create double-stranded blunt-ended cDNA fragments. The cDNA fragments were extended by adenylating the 3′ blunt-end to enable the ligation of sequencing adapters. The adapter-containing cDNA fragments (approximately 200 bp) underwent agarose gel electrophoresis, and fragments were isolated based on size. The cDNA fragments were amplified by PCR using adapter-specific primers. The cDNA library was quantified using the KAPA Library Quantification Kit (Kapa Biosystems, USA) according to the manufacturer’s instructions. The library was used for high-throughput sequencing with the Illumina NextSeq platform, as well as a paired-end sequencing system, to generate raw sequence data.

De novo assembly and sequence analysis

Raw read data were acquired by image analysis and base calling using the Illumina Pipeline software. Raw reads were quality checked with the quality assessment software FastQC (v0.11.5) (Andrews, 2010). Raw reads were trimmed, and clean reads were obtained by removing low quality reads (Phred quality score: Q ≥ 30 for all bases), and adapter sequences were eliminated using Skewer (v0.2.2) (Jiang et al., 2014). De novo sequences were assembled using Velvet (v1.2.07) (Zerbino and Birney, 2008) and Oases (v0.2.08) (Schulz et al., 2012) based on the de Bruijn graph algorithm, as described in Kim et al (Kim et al., 2016). Velvet-Oases were used to assemble short read sequences, and they obtained a higher quality assembly using reads with long N50 lengths. The optimization of program parameters to the N50 length and the trimming of low-quality reads could improve the assembly output significantly (Garg et al., 2011). Velvet was used to assemble contigs from the clean reads with various k-mer lengths (i.e., 51, 53, 55, 57, 59, 61, 63, 65, 69, 71, 75, 79, 81, 85, and 89). The contigs were assembled into transcripts for each k-mer using Oases. After transcripts were assembled for the 57 and 59 k-mer lengths (i.e., optimal k-mer lengths), a final assembly was completed for the hash length (BioProject Accession number: PRJNA476484). Putative transcripts assembled from the total reads of each mRNA sample were confirmed with sequences available in the Phytozome database (http://www.phytozome.net/) using the BLASTx algorithm (E-value ≤ 1E-10).

Functional annotations of differentially expressed genes

Paired-end clean reads were aligned with the assembled loci and expression was quantified by counting the number of mapped clean reads using Kallisto (v0.43.1), which is based on the novel idea of pseudoalignment for rapidly determining the compatibility of reads with targets (Bray et al., 2016). Genes whose expression levels were affected during the PEG-treatment period were identified using NOISeq (v2.16.0) for normalization (TMM normalization of the TPM from Kallisto output) and differentially expressed gene (DEG) analysis (Tarazona et al., 2015). The probability ≥ 0.95 and log2 values ≥ 2 were used to calculate the significance of the altered expression levels with NOIseq.

Functional annotations were completed using the BLASTP algorithm (E-value ≤ 1.0E-10). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Moriya et al., 2007) were analyzed using DAVID (http://david.abcc.ncifcrf.gov/tools.jsp) (Huang et al., 2008). The enrichment p-values were derived from an EASE score, a modified Fisher’s exact p-value, for gene-enrichment analysis. The group enrichment score was calculated from the geometric mean (in −log scale) of all the enrichment p-values. Visualization of the analysis results was shown using NOIseq (hierarchical clustering of samples and various plots of gene expression profile) and In-house R script (heatmap of clustering).

Validation of differentially expressed genes by quantitative real-time PCR

The relative mRNA expression levels of DEGs in drought conditions were analyzed by qRT-PCR using the SYBR Green Master Mix (Enzynomics Co., Korea) and the CFX Connect Real-Time PCR System (Bio-Rad, USA) according to the manufacturer’s instructions. Total RNA (1 μg) was used to synthesize cDNA with the PrimeScript RT Reagent Kit (Takara, Shiga Japan). The qRT-PCR was completed using 1 μl cDNA, gene specific primers (Supplementary Table S1), and the following program: 95°C for 15 min; 45 cycles of 95°C for 10 s, 55°C for 15 s, and 72°C for 30 s. Relative expression levels were determined by normalizing the data for the target transcripts against the Actin transcript data (GenBank accession number: XM_006345899) according to a normalized expression method. Three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed.

Isolation and characterization of the drought-tolerant potato mutants

The 344^th M1 plants of the 1,000 EMS-mutagenized population survived under the no watering condition for 4 wks. This surviving plant was healthy with no drought-related symptoms (Fig. 1A) and tubers were harvested as normal. To verify the plants response to drought stress in detail, in vitro plants of WT and the 344^th plant were cultured in the MS basal medium containing PEG as an osmotic pressure agent. After culturing for 10 days, growth rates including the stems and roots of both types of plant were not different under the normal condition (Fig. 1B). The growth rate of the WT plants began to decrease in the drought stress condition. Additionally, root growth was not observed, while chlorotic leaf lesions developed and eventually became necrotic (Fig. 1C). In contrast, the 344^th plant exhibited normal leaf and shoot growth, and root growth was especially vigorous. This plant was regarded as drought tolerant, and we referred to it as the DR plant. Interestingly, more stomata were closed in the DR plant than in the wild-type plant under the same culture conditions, otherwise no significant differences in the root surface in both types of plant were observed (Figs. 1D–1E).

RNA-Seq analysis and de novo transcriptome assembly

The expression profiling of WT and DR plants in response to a PEG-induced drought was performed. In total, 306,308,385 raw reads were produced from the cDNA libraries of 10 samples of these two types of potato at five time points. After trimming adapters and low-quality reads, we obtained 304,976,315 clean reads with a mean quality score of Q34.73 (Phred quality score: Q ≥ 30). Of the 10 samples, WT reads were typically mapped to the potato reference sequence using Bowtie2 software for reference CDSs (PGSC_DM_v4.03) (Xu et al., 2011). The total mapping rate was only 51.94%, with 32.73% of the reads aligned once, and 19.21% of the reads aligned more than once (Supplementary Table S2). Furthermore, the reference transcript data has 38,982 predicted gene models. In this study, we identified expression of 61,100 transcripts that were difficult to cover with reference transcripts data. Thus, the remaining 48% transcripts were uncharacterized owing to their uniqueness or incompatibility with the current reference database. Clean reads were de novo assembled into the transcriptome containing 61,100 transcripts with over 200-bp transcript lengths and N50 length of 1,667 bp transcripts. The size distribution of the assembled transcripts are shown in Supplementary Fig. S1.

Analysis of differentially expressed genes in response to drought stress

To identify the differentially expressed gene levels between WT and DR under drought stress, clean reads were mapped to the assembled potato transcriptome sequences using Kallisto. The total average mapping rate was 85.09% with 258,649,240 read pairs in the 306,308,385 read pairs of the WT and DR plants. Transcriptome profiling revealed a total number of 909, 977, 1181, 1225, and 826 DEGs (absolute log2 fold change ≥ 2) between WT and DR plants at 0, 6, 12, 24, and 48 h, respectively. Venn diagram data is expressed using DEGs between wild-type and mutant corresponding to all conditions including normal condition (Fig. 2A). To understand the mutant-specific expression profiles, differences in transcript expression were determined by pairwise comparisons between WT and DR plants without drought stress (at 0 h), a total of 909 DEGs were found to be significantly (absolute log2 fold change ≥ 2, probability ≥ 0.95) differentially expressed with 299 up-regulated and 610 down-regulated DEGs (Fig. 2B). An increase in the total number of DEGs between WT and DR plants was observed within one day of the drought stress condition (at 6–24 h). A total of 977 DEGs (518 up-regulated and 459 down-regulated) were detected after 6 h of PEG treatment in WT and DR plants (WT6 vs. DR6), and 1181 DEGs (848 up-regulated and 333 down-regulated) were detected after 12 h (WT12 vs. DR12). Interestingly, the largest number of up-regulated DEGs and the least number of down-regulated DEGs were observed at 12 h of drought stress. To compare gene expression of each pairwise comparison, log10 (FPKM) values of DEGs are shown as scatter plots (Supplementary Fig. S2A). For these plots, ‘significant’ (marked as red color spots) represents a probability ≥ 0.95. The total number of DEGs in WT0 vs. DR0 was similar to other pairwise comparisons, but there were a few genes with significantly high expression levels. In addition, the highest expressed gene of WT0 vs. DR0 was found to be unknown, expressed significantly higher (7.958 log2 fold change) than other DEGs. The gene expression levels of WT and DR plants at 6 and 24 h after PEG treatment were significantly higher than other conditions.

To understand the transcriptional changes of DR plants in response to drought stress, tolerant-specific expression genes were determined by comparing the PEG-treated to non-treated DR plant (at 0 h). Transcriptome profiling revealed a total number of 3213, 2513, 2710, and 3158 DEGs in the drought tolerant DR plants at 6, 12, 24, and 48 h, respectively (Fig. 2C). We found that the number of DEGs was the highest in DR plants at 6 h (3213 DEGs; 931 up-regulated and 2282 down-regulated), compared to other time course DEGs. In addition, the highest number of up-regulated DEGs occurred at 6 h and gradually decreased until 48 h. A large amount of down-regulated DEGs were observed at 6 h, but the number of DEGs significantly dropped at 12 h (2513 DEGs; 866 up-regulated and 1647 down-regulated). Conversely, a gradual increase in DEGs was observed after 12 h and the highest number of DEGs was observed at 48 h (3158 DEGs; 711 up-regulated and 2447 down-regulated). Overall, the number of down-regulated genes in DR plants was greater than the number of up-regulated genes during drought stress, but the level of each gene expression was significantly higher in up-regulated genes (Supplementary Fig. S2B). The combined analysis were carried out to show the genes involved in the DEGs of DR plants under each drought conditions against identified DEGs between wild-type and DR plant in the same conditions (Fig. 2D). According to the drought conditions, 28, 19, 84, and 31 DEGs were matched between DR_control (DR0) vs. DR and wild-type vs. mutant, respectively (Supplementary Table S3).

Gene Ontology classification analysis of differentially expressed genes

To investigate the potential DEGs activated in response to drought stress in the tolerant plant compared to the wild-type, GO classification was performed on DEGs identified in each pairwise comparison (Fig. 3 and Supplementary Fig. S3). GO information in the assembled potato transcripts was generated by DAVID with a cut-off probability ≥ 0.95 and absolute log2 fold change ≥ 2. In the GO enrichment analysis of total DEGs including up-regulated and down-regulated, three GO categories (Defense response; GO:0006952 in Biological Process, Integral component of membrane; GO: 0016021 in Cellular Component, ADP binding; GO:0043531 in Molecular Function) were found to be related to drought tolerance during all drought conditions (Fig. 3A). Interestingly, these three GO categories included over 75% of up-regulated DEGs in the early hours (0–12 h), and over 50% down-regulated DEGs in the late hours (24–48 h) (Supplementary Fig. S4). A total of 9 GO categories with up-regulated and down-regulated DEGs were identified in two plant types in the normal condition (WT0 vs. DR0). Of these, the four GO categories (ATP binding, ADP binding, integral component of membrane, and defense response) were found to be rich not only among up-regulated DEGs, but also among down-regulated DEGs. However, there are no significant enriched categories (p-value ≤ 0.05) in these 9 GO categories between WT and DR plants under normal conditions (Supplementary Fig. S3A). On the other hand, significantly up-regulated DEGs were observed in the two GO categories (‘Terpene synthase activity; GO:0010333’ and ADP binding in Molecular Function) in DR plants at 6 h in the drought condition. The biggest change in DR plants was observed in GO enriched analysis at 24 h (WT24 vs. DR24) with the highest amount of DEGs. In the WT24 vs. DR24 pairwise comparison, a total of 45 GO categories with 15 significant enriched GO categories (p-value ≤ 0.05) were identified compared to other pairwise comparisons (Supplementary Fig. S3D). The water response (GO:0009415) of the GO biological process containing dehydrin (LOC102599708), a subfamily of LEA protein known to be involved in drought stress, showed a significant increase (5.387 log2 fold change, 1.00 probability) in DR24 compared to WT24 (Fig. 3B). Regarding plant photosynthesis, there were significant differences between WT and DR plants after 24 h of drought stress. A category of Molecular Function (Chlorophyll binding; GO:0016168), three categories of Cellular Component (Chloroplast thylakoid membrane; GO:0009535, Photosys-tem II; GO:0009523, and Photosystem I; GO:0009522), and two categories of Biological Process (Protein chromophore linkage; GO:0018298, and Photosynthesis, light harvesting; GO:0009765) had a highly significant association (−log10 (p-value) ≥ 10) in the down-regulated DEGs (Fig. 3B). Chlorophyll a–b binding proteins (LOC102603980, LOC102586284, LOC102583646, and LOC102589865 were expressed at −3.532, −2.618, −2.504, and −2.357 log2 fold change, respectively) enriched in the photosynthetic GO categories including photosystem I, photosystem II, photosynthesis, light harvesting, protein~chromophore linkage, and chlorophyll binding were found to have low expression in the DR plant at 24 h.

We also performed GO classification on DEG identified in each combined analysis “DR_control (DR0) vs. DR and WT vs. DR” (Fig. 3C). WT6 vs. DR6/DR0 vs. DR6 containing 28 DEGs was enriched in “terpene synthase activity” and “protein serine”, and only one category “response to biotic stimulus” was enriched in WT48 vs. DR48/DR0 vs. DR48 containing 31 DEGs. In WT24 vs. DR24/DR0 vs. DR24 containing 84 DEGs, it was enriched in 10 categories. In particular, as in the comparison of WT24 vs. DR24, the photosynthesis-related GO categories showed considerably higher enrichment.

Metabolic pathway analysis by KEGG

When plants are affected by drought stress, many genes are expressed or inhibited to protect plants through associated pathways. To identify the biological pathway of DEG associated with drought, all DEGs were compared against the KEGG using DAVID with a cut-off probability ≥ 0.95 and absolute log2 fold change ≥ 2 (Fig. 4 and Supplementary Fig. S5). In total, all DEGs were enriched in 33 different functional pathway categories with 20 categories of up-regulation and 18 categories of down-regulation in all pairwise comparisons. Of these, 12 different functional pathway categories are significantly enriched categories (p-value ≤ 0.05). In particular, the most important biochemical metabolic pathway was identified by the KEGG pathway. KEGG pathway analysis indicated that ‘Photosynthesis-antenna protein, sot00196’ was enriched most significantly (−log10 (p-value): 17.82) in WT24 vs. DR24 (Fig. 4B). Next, ‘protein processing in endoplasmic reticulum, sot04141’ showed a significant association (−log10 (p-value): 8.86) in WT12 vs. DR12 (Fig. 4A). These two different functional pathway categories were the most relevant compared to the other categories of WT and DR plants under PEG treatment, but WT0 vs. DR0 had no significant enrichment category, as with the GO analysis. Interestingly, two different functional pathways (‘Biosynthesis of secondary metabolites’ and ‘Metabolic pathways’) were found to be enriched not only among up-regulated DEGs, but also among down-regulated DEGs of all PEG treatment groups. Results of KEGG enrichment in this study indicate that the mechanism for drought tolerance in the DR plant can be mainly explained in two aspects, photosynthesis-antenna protein and protein processing in the endoplasmic reticulum.

In the combined analysis, KEGG pathway analysis of the DEGs of DR_control (DR0) vs. DR and WT vs. DR indicated that only “Photosynthesis–antenna proteins” and “Metabolic pathways” were enriched in WT24 vs. DR24/DR0 vs. DR24 (Fig. 4C). In particular, the “Photosynthesis-antenna proteins” pathway were significantly related in the drought condition. This result was similar to the down-regulation result of the KEGG analysis of WT24 vs. DR24, indicating that the analysis of DEGs between WT and DR in each drought conditions involves the results of a combined analysis with DR_control vs. DR under each drought conditions.

To investigate the functional significance of 84 DEGs of WT24 vs. DR24/DR0 vs. DR24, we focused on playing a role in the well-characterized pathway of photosynthesis using the MapMan software (Fig. 5A) (Thimm et al., 2004). Among these DEGs, photosynthesis was affected by drought stress, with down-regulated 10 DEGs involved in the light reaction, down-regulated 2 DEGs in the Calvin cycle, and down-regulated 1 DEG in the photorespiration under drought condition after 24 h. In the KEGG map of the “photosynthesis-antenna protein” pathway, down-regulated 10 DEGs corresponded to the Light-harvesting chlorophyll protein complex (LHC) (Fig. 5B). Among these 10 DEGs, 3 DEGs corresponded to light-harvesting complex I chlorophyll a/b binding protein 2-4 (Lhca2-4) of light-harvesting complex I, and 7 DEGs corresponded to light-harvesting complex II chlorophyll a/b binding protein 1-3 (Lhcb1-3) of the light-harvesting complex II. The KEGG analyses of WT24 vs. DR24 and combined analysis (WT24 vs. DR24/DR0 vs. DR24) showed the same results that expression of these LHC genes decreased.

DEG clustering analysis of DR plants and validation by qRT-PCR

To identify different expression patterns of DR plants under drought conditions, four pairwise comparisons (DR0 vs DR6, DR12, DR24, DR48) were performed using five libraries. Then, 4,395 DEGs identified in these pairwise comparisons were grouped using a clustering algorithm with Pearson’s correlation, and a heatmap and eight clusters were identified (Fig. 6). Cluster 1 and 2 represent up-regulated clusters, which contain 519 and 566 genes, respectively. Cluster 3 included 247 genes which showed up-regulation until 24 h but was not up-regulated in DR48. Cluster 4 included 81 genes, which were only down-regulated in DR6, whereas Cluster 7 included 60 genes that were up-regulated in DR6. Cluster 5 and 6 represented down-regulated clusters, which contained 1830 and 980 genes, respectively. These clusters contained over 60% of the total DEGs. Representative genes of each cluster are shown in Table 1. In this result, the families of LEA proteins, including dehydrin related to the water response mentioned in the GO analysis, showed an up-regulated expression pattern of Cluster 1. In addition, genes involved in the photosystem (e.g. Photosystem I subunit), which had a large impact on the DR plant under drought stress, showed a down-regulated expression pattern of Cluster 5. These results suggest that the genes involved in the water response and photosystem play a major role in the DR plant in the early stages of the drought environment.

To confirm the reproducibility of the RNA-seq by qRT-PCR validations, eight DEGs were randomly selected based on the time-course clustering analysis data in the DR plant (Fig. 7). CAP160 and LEA proteins belong to Cluster 1 of the RNA-seq expression pattern, which showed the highest expression at 6 h, and also showed the highest expression in the qRT-PCR result at 6 h. In addition, expansin-like B1 and laccase 14 proteins belong to Cluster 2, which showed the highest expression pattern at 48 h, and also showed the same expression pattern in qRT-PCR. As a result, the selected candidate genes determined by qRT-PCR were closely matched with the RNA-Seq results. This suggests that the data obtained from the DEGs analysis were credible.

Many studies to increase drought tolerance in potatoes have been performed (Kikuchi et al., 2015), and most of the work has been undertaken using transgenic techniques. However, the degree of tolerance to abiotic stress is not easily controlled due to the complexity of the genetic regulatory mechanism. Traditionally, mutants are not only very useful for the development of new varieties but are also a suitable method of research widely used in gene function studies. Genotypic variations of salt-tolerant potato mutants developed by irradiation were characterized using a Random Amplified Polymorphic DNA-PCR analysis (Yaycili and Alikamanoglu, 2012), and a network-based transcriptome analysis related to biological pathway modifications in irradiated rice mutants was performed based on microarray transcriptional profiling (Hwang et al., 2015). We found that DR plants obtained through EMS treatment had obvious drought tolerance under severe drought conditions, although there was no significant difference in plant growth and tuber yield compared to WT under the normal condition. In addition, transcriptome analysis was performed based on RNA-seq profiling to identify changes in the expression of drought-related genes using DR plants. Our potato (WT0) sequencing data was mapped to the published reference potato genome (DM1-3 516 R44) prior to de novo assembly, resulting in a slightly lower mapping rate. The published reference potato genome was analyzed against the doubled-monoploid potato genome, and many researchers have used this potato genome sequence since (Anithakumari et al., 2012; Gong et al., 2015; Uitdewilligen et al., 2013; Zhang et al., 2014). Using Illumina sequencing technology, we generated sequence data (61,100 transcripts) by de novo transcriptome assembly and identified drought response-related genes in DR plants. Based on the assembled transcriptome data, 909, 977, 1181, 1225, and 826 DEGs were identified in the comparison of WT and DR plants by time conditions (0, 6, 12, 24, 48 h after PEG treatment), respectively.

In the GO annotation in WT and DR plants, DEGs enrichment was observed in the photosynthesis-related category. Photosynthesis is one of the main processes affected by drought stress (Chaves, 1991). Among them, it has been reported that chlorophyll concentration, which is one of the evaluation indexes of environments, is decreased in drought stress conditions and it can be considered as a nonstomatal limiting factor (Khayatnezhad and Gholamin, 2012). Moreover, the decrease in expression of chlorophyll a–b binding proteins was consistent with the results reported previously in tomatoes under drought stress (Iovieno et al., 2016). The antenna protein present in light-harvesting chlorophyll protein complexes of photosynthetic plant act as a peripheral antenna system, allowing more efficient absorption of light energy. MapMan analysis and KEGG map analysis showed that the DR plant had decreased expression of the LHC subunits (Lhca2-4, Lhcb1-3) in drought stress than in the WT plant (Figs. 5A and 5B). In order words, the expression of chlorophyll a–b binding proteins was decreased under drought stress in DR plants, meaning that the synthesis of factors related to photosynthesis is reduced, and this reaction affects the photosynthesis system. Next, the reduction of photosynthesis limits water consumption through plant growth inhibition, thereby maintaining the water status of the plant and assisting the plant in carbon assimilation to survive under stress. Overall, the DEGs of DR plants after 24 h of PEG treatment showed significant down-regulation patterns in six enriched GO categories associated with photosynthesis-related pathways and the same down-regulation pattern in the ‘photosynthesis-antenna proteins’ functional pathway of KEGG pathway analysis. This suggests that the mutation of the photosynthetic genes mentioned in the GO analysis led to changes in the photosynthetic metabolism of the DR plant under drought stress. These results are again consistent with previous reports that photosynthesis related genes are down-regulated under drought stress conditions. These results also suggest that drought stress affects plant photosynthesis, as it occurs in other plant species including S. moorcroftiana (Li et al., 2015), C. morifolium (Xu et al., 2013), B. nivea (Liu et al., 2013b), and A. mongolicus (Liu et al., 2013a).

Inhibition of photosynthesis in a drought environment leads to a decrease in stomatal conductance, thereby improving water use efficiency. Especially, inhibition of CO₂ supply to rubisco causes down-regulation of photosynthesis when plants are exposed to high light and temperature. It leads to a decrease in stomatal conductance as shown in our results (Figs. 1D–1E), which can protect plants against stress by improving the efficiency of water consumption utilization (Chaves et al., 2009). Our result which showed that the rubisco activase gene induced significant down-regulation (−2.408 log2 fold change) in DR plants after 48 h of PEG treatment again supports the above-mentioned results. On the other hand, up-regulated DEGs in DR plants were observed in the water response (GO:0009415) category of the GO biological process (Fig. 3B). The LEA proteins, which are protection factors for macromolecules showing protein kinase and phosphatase activity under osmotic stress, are related to tolerance against water deficiency, as well as osmotic, salt, and cold stresses (Sivamani et al., 2000; Zhang et al., 2014). In our results, dehydrin (LOC102599708), a subfamily of LEA protein, showed a significant increase (5.387 log2 fold change) in DR compared to WT, again providing further evidence that the LEA protein is related to drought tolerance in potatoes.

Potato multicystatin, a member of the cystatin family and a subfamily of phytocystatin, is known to act as an inhibitor of Cys-proteases involved in the degradation of storage proteins, senescence, programmed cell death, and stress signaling. Thus, cystatin plays an important role in endogenous process regulation and plant protection from various environmental stresses by regulating the activity of Cys-proteases, but the underlying mechanisms are limited (Tan et al., 2017). In addition, overexpression of cystatins in transformed Arabidopsis (A. thaliana) showed enhanced tolerance to oxidative stresses, high salt, and drought stress (Zhang et al., 2008). In the DR plant, the multicystatin (Gi number: 6671196) gene was significantly higher (4.163 log2 fold change) than in the WT plant after 48 h under drought stress.

Conversely, several drought-related genes in DR plants showed different expression patterns than previously known gene information. The polygalacturonase-1 non-catalytic subunit beta (Gi Number: 565399375) and GDSL esterase/lipase 2 (Gi Number: 565344594) showed the same expression patterns in WT and DR plants for drought stress and induced down-regulation in all time conditions of PEG treatment (data not shown). PG-1 beta is a subunit of polygalacturonase (PG), one of the hydrolytic enzymes involved in cell wall pectin degradation. Overexpression of PG-1 beta induced pectin degradation and affected the integrity and proliferation of the cell wall, thereby reducing tolerance to abiotic stress (Liu et al., 2014). The GDSL esterase/lipase 2 gene involved in plant development, morphogenesis, synthesis of secondary metabolites, and defense responses have been reported to play an important role in abiotic stress including drought in rice (Chepyshko et al., 2012). These results demonstrate that these genes had no effect on the drought tolerance of DR plants, although the PG1 beta gene may have benefited both plants. Another case is that the Defense response (GO: 0006952) category of Biological Process in the DR plant was significantly down-regulated at 12 h after drought stress. Many defense related genes including NB-LRR were significantly increased after drought and salt treatments (Li et al., 2017), which is in contradiction to our results. This inconsistency needs to be confirmed by conducting the disease tolerance test in the future, whether the drought-tolerant DR plants we obtained have weaker disease tolerance than WT plants.

In comparative analysis of gene expression in DR plants by time course under drought stress, increased expression of the LEA family protein, expansin-like B1, and laccase 14 affected the enhanced drought tolerance in the DR plant (Fig. 7 and Table 1). The LEA family protein has already been mentioned above. Expansins are a family of proteins involved in cell wall loosening and cell enlargement during plant growth and development (Sampedro and Cosgrove, 2005; Zhao et al., 2012). In Jatropha subjected to drought stress, two kinds of expansins were highly expressed during the recovery phase (Cartagena et al., 2014). Although the process of cell wall loosening and expansion involves other genes or gene complexes, the up-regulation of expansin transcripts during recovery from stress is understandable. Laccase transcript levels are enhanced in non-lignifying cell types under abiotic stresses, such as drought (Cho et al., 2014). Arabidopsis plants containing three laccase gene mutations showed no root growth under PEG-induced drought conditions (Cai et al., 2006). Here, the expansin-like B1 and laccase 14 genes were continuously highly induced during 48 h of PEG treatment. Based on the up- and down-regulated gene analyses, some genes are consistent with previous reports, but others showed distinctly different expression patterns. For example, UDP-glucuronosyl-transferase (UGT), involved in flavonol glycoside biosynthesis and defense roles against stress, was strongly up-regulated under drought stress in potato leaves (Zhang et al., 2014), but in our results, UGT was the most down-regulated gene. Germin proteins, which have enzymatic activities, including oxalate oxidase and superoxide dismutase, that are involved in drought tolerance mechanisms (Lane et al., 1993), were found in a gene expression analysis of germination in wheat (Lane et al., 1991). The germin 3 protein shows a high sequence homology to auxin binding protein 19a, and it has a role in changing the cell wall characteristics associated with normal growth. Germin 3 is repressed by drought stress in Arabidopsis leaves (Bray, 2004). However, germin3/auxin-binding protein ABP 19a-like was down-regulated here, in contrast to our expectations. Because drought tolerance is a complex trait that involves intricate molecular mechanisms related to many genes, more detailed studies focused on these inconsistencies are needed in the near future.