Introduction {#Sec1}
============

Purple-fleshed sweet potato (PFSP) (*Ipomoea batatas* (L.) Lam*.*) is a type of sweet potato, which has purple skin and flesh color of the storage roots due to the large accumulation of anthocyanin pigments^[@CR1]^. It is a root crop widely grown in some Asian and African countries (i.e. China, India, and Kenya). Sweet potato is regarded as the seventh most important food crop. It is among the crops which generate large amounts of food per unit area per unit time^[@CR2]^. Furthermore, it is considered to be a cheap source of energy, carbohydrates, vitamin, potassium, iron, fiber, and protein^[@CR3]^. Purple-fleshed sweet potato (PFSP) is also an important source of anthocyanin with strong antioxidant properties. Therefore, it can be a quick solution for nutritional problems in many developing countries^[@CR4]^. Due to high anthocyanin content, sweet potato shows excellent antioxidant properties.

Fluctuations in climatic conditions and unprecedented rainfall have resulted in increases in the exposure of crop plants to different types of abiotic stresses, such as drought, salt, high temperature, and cold stress. Among these abiotic stresses, drought is the major yield-limiting factor in a variety of crop plants, including sweet potato. Drought stress causes plant water deficits in sweet potato, which results in decreased cell turgor and cell enlargement, closed stomata, decreased photosynthetic rate, and shortened vegetative period^[@CR5]^. Moreover, it adversely affects anatomical traits of plants, physiological processes, biochemical pathways and ultimately inhibits plant growth and incurs severe yield losses^[@CR6]^. Drought tolerance in sweet potato remains a challenge for plant breeders due to its large homogeneity of the chromosome number (2n = 6X = 90) with the estimated genome size of more than 2.4 GB. Furthermore, it expresses self-incompatibility with greater levels of cross-incompatibility^[@CR7]^.

In non-model plants without a reference genome such as sweet potato, transcriptome sequencing approaches adopt a rapid method to improve drought tolerance through genetic manipulation^[@CR8]^. Transcriptome sequencing is considered an effective strategy for studying the expression of a large number of genes under drought stress conditions^[@CR9]^. Furthermore, RNA sequencing using high throughput next-generation and third-generation sequencing technologies can be used for the discovery of genes involved in drought tolerance and their related signaling pathways^[@CR10]^. Transcriptome studies of a few cultivars of sweet potato have been conducted by de novo transcriptome assembly using second and third-generation sequencing technologies with Illumina platforms^[@CR11]--[@CR18]^. However, transcriptome sequencing was mostly performed from Ipomoea trifida (2n = 2x = 30), the closest wild relative of Ipomoea batatas, thus making it the sweet potato progenitor^[@CR19]^.

In the present study, second and third-generation sequencing technologies were used to establish a useful database for transcriptome sequencing and differentially expressed genes in purple-fleshed sweet potato leaves under different durations of drought stress. A total of 184,259,679 high-quality reads were assembled into 18,209 transcripts and 17,508 unigenes. Our results provide novel insights into sweet potato response dynamics of changes to drought stress and revealed potential defense mechanisms of specific genes involved in drought tolerance under different time points. Their role in drought tolerance can help guide future efforts to breed for drought-tolerance in sweet potato.

Materials and methods {#Sec2}
=====================

Plant materials and growth conditions {#Sec3}
-------------------------------------

In this study, Xuzi-8, a high yielding drought tolerant and early maturing sweet potato cultivar with storage roots having purple flesh was used. Its storage roots contain about 6% soluble sugar content and more than 80 mg anthocyanin content/100 g fresh weight. This cultivar was obtained from Xuzhou Institute of Agricultural Sciences in Jiangsu Xuhuai District, China.

This sweet potato cultivar was propagated by apical stem cuttings of 15--20 cm length. The cuttings were grown hydroponically in the Hoagland nutrient medium for four weeks. The cultures were incubated at 25 ± 3 °C under 16 h light/8 h dark photoperiod. Seedlings were then exposed to drought stress by the addition of 30% PEG-6000 to Hoagland solution. Leaf samples (0.1 g) were collected from the upper third of the plant at 0 (control), 1, 6, 12 and 48 h after drought treatment and were used for next-generation sequencing. However, leaf samples taken from the upper third of sweet potato seedlings after 1, 6, 12 and 48 h of drought exposure were mixed and used for third-generation sequencing; untreated control samples were collected and examined in parallel. All samples were frozen in liquid nitrogen immediately upon collection and stored at − 80 °C until processing for RNA extraction. Each sample was analyzed in triplicate.

Total RNA isolation {#Sec4}
-------------------

Total RNA was extracted from both drought treated and control plant leaves using the TRIZOL reagent (Life Technologies, Carlsbad, CA) according to the manufacturer's protocol. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using the Nanophotometer spectrophotometer (IMPLEN, CA, USA). The Qubit RNA Assay Kits were used with the Qubit 2.0 Fluorometer (Life Technologies, CA, USA) to determine RNA concentration. The RNA integrity was assessed using the RNA 6,000 Nano Kit of the Agilent 2,100 Bio-analyzer system (Agilent Technologies, CA, USA).

cDNA library preparation for transcriptome analysis {#Sec5}
---------------------------------------------------

For second generation sequencing, a total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations. NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) was used for generating sequencing libraries following the manufacturer's protocol and index codes were added to attribute sequences to each sample^[@CR20]^. For third generation sequencing RNA was extracted from different time points and equal amount of each sample were mixed to produce one library beside to the control library. Briefly, poly-T oligo-attached magnetic beads was used to purify mRNA from total RNA. NEBNext First Strand Synthesis Reaction Buffer (5X) under elevated temperature was used for fragmentation through divalent cations. First-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-)^[@CR21]^. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After that, NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. The library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA) to selected the fragments more than 250--200 bp. After that, 3 µL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Universal PCR primers and Index (X) Primer. Finally, PCR products were purified and library quality was assessed on the Agilent Bioanalyzer 2,100 system.

Sequence clustering and analysis {#Sec6}
--------------------------------

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-H (Illumina) following the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina platform and paired-end reads were generated.

Quality control data analysis {#Sec7}
-----------------------------

At first, raw data (raw reads) of fastq format were processed through in-house Perl scripts. In this step, reads containing adapter, ploy-N, and low-quality reads were removed from raw data to obtain clean data (clean reads)^[@CR22]^. The Q20, Q30, GC-content and sequence duplication level of the clean data were simultaneously calculated. All downstream analyses were based on clean, high-quality data^[@CR23]^.

Transcriptome assembly and gene functional annotation {#Sec8}
-----------------------------------------------------

Transcriptome sequencing was accomplished based on both NGS and 3rd GS, and TPM, FPKM, RPKM and fold change (FC) were recorded for each replicate of each library separately. Obtained sequence from NGS and 3rd GS was aligned and similar sequence data from all libraries/samples were pooled^[@CR24]^ and ultimately used for further analysis using Trinity^[@CR25]^ with min_kmer_cov set to 2 by default and all other parameters set default.

The functional annotation of genes was performed based on different source databases such as NCBI non-redundant protein sequences (Nr) (<https://www.ncbi.nlm.nih.gov/>) with an e-value cut-off of 1E-5; Nt (NCBI non-redundant nucleotide sequences) (<https://www.ncbi.nlm.nih.gov/>) with an e-value cut-off of 1E-5; Pfam (Protein families), (<https://pfam.sanger.ac.uk/>) with an e-value cut-off of 1E--2; KOG/COG (Clusters of Orthologous Groups of proteins) (<https://www.ncbi.nlm.nih.gov/cog/>) with an e-value cut-off of 1E--3 and Swiss-Prot (reviewed) (A manually annotated and non-redundant protein sequence database) (<https://www.ebi.ac.uk/uniprot/>) with an e-value cut-off of 1E--5.

Based on the annotation result of NR and Pfam, Blast2GO (v2.5) was used to determine GO (Gene Ontology) annotations (<https://www.geneontlology.org>) according to the three GO ontologies (molecular function, biological process, and cellular component)^[@CR26]^. GO enrichment analysis of the differentially expressed genes (DEGs) was performed by the GO-seq R packages based on Wallenius non-central hypergeometric distribution^[@CR27]^, which can adjust to gene length bias in DEGs.

Kyoto Encyclopedia of Genes and Genomes (KEGG)^[@CR28]^ is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (<https://www.genome.jp/kegg/>). The KOBAS software was used^[@CR29]^ to test the statistical enrichment of DEGs in KEGG pathways.

Quantification of gene expression levels and differential expression analysis {#Sec9}
-----------------------------------------------------------------------------

Gene expression levels were estimated using RSEM^[@CR30]^ for each sample. Clean data were mapped back onto the assembled transcriptome. Read counts for each gene were obtained from the mapping results.

Differential expression analysis of two conditions/groups was performed using the DESeq2R differential expression package for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The obtained *P* values were adjusted using the Benjamini and Hochberg method for controlling the false discovery rate (FDR). Genes with the FDR-adjusted *P* value of ≤ 0.05 were defined as differentially expressed genes (DEGs)^[@CR31]^. FDR ≤ 0.05 and log~2~ of the fold change (\|log~2~ FC\|) ≥ 1 were set as the thresholds for significant differential expression.

Results {#Sec10}
=======

Transcriptome sequencing, assembly, and mapping {#Sec11}
-----------------------------------------------

In order to conduct gene expression profiling of hexaploid sweetpotato exposed to PEG-induced drought stress at different time points, ten cDNA libraries were constructed for NGS from leaves of treated (with 30% PEG-6000) and untreated seedlings after 0, 1, 6, 12, and 48 h. For 3rd GS, four time points RNA samples including1, 6, 12 and 48 h were mixed to produce one library beside to the control library. These libraries were separately sequenced using Illumina high-throughput third generation sequencing platform. TPM, FPKM, RPKM and fold change (FC) were recorded for each replicate of each library separately on both NGS and 3rd GS. Obtained sequence from NGS and 3rd GS were aligned and similar sequence data from all libraries/samples were pooled. Due to the lack of a reference genome, the clean reads resulted in from the transcriptome sequences were aligned and assembled using Trinity software. In total, approximately 184,259,679 clean reads were obtained from the deep sequencing with 34--38 million reads per library for all samples and mapped reads ratios of 59.38, 61.08, 61.52, 57.73, and 57.24% at 0, 1, 6, 12, and 48 h after treatment, respectively (Table [1](#Tab1){ref-type="table"}). The GC content ranged between 45--46%, Q20 ranged between 96--97% and the proportion of unknown nucleotides in clean reads (N %) was 0%. A total of 17,508 unigenes and 18,209 transcripts were obtained from all libraries with average lengths of 1,783 and 1785 base pairs (bp), respectively (Table [2](#Tab2){ref-type="table"}). Furthermore, a total of 1827 (10.4%) unigenes and 1931 (10.6%) transcripts were longer than 3,000 bp and 84% of unigenes and transcripts were longer than 1,000 bp (Fig. [1](#Fig1){ref-type="fig"}A).Table 1Next generation sequencing statistical summary of sequenced and assembled results.0 h1 h6 h12 h48 hTotal reads38,990,28837,595,49437,441,88134,736,08535,495,931Mapped reads23,178,49022,970,72423,021,70420,039,97620,308,985Mapped reads (%)59.3861.0861.5257.7357.24Nt11,697,086,54011,278,648,20011,232,564,30010,420,825,70010,648,779,500GC (%)46.6046.5146.7345.5345.59Q20 (%)96.0196.7696.6097.0396.81Q30 (%)92.0192.0091.5292.4692.07N%00.0000.0000.0000.0000.00*Nt*, total number of clean nucleotides; *The GC%* is the proportion of guanidine and cytosine nucleotides among total nucleotides; *The Q20 and Q30%* is the proportion of nucleotides with a quality value \> 20 and 30, respectively; *The N%* is the proportion of unknown nucleotides in clean reads. Table 2Third generation sequencing statistical summary of sequenced and assembled results.UnigenesTranscriptsTotal number of sequences17,50818,209Total sequences length31,222,20932,513,465Maximum length8,4558,455Minimum length227227Average length1,7831,785GC (%)43.7243.74N402,1382,147N501,8531,857N601,6501,651N701,4541,455N801,2501,250N901,0521,052*N50*, represents sorting the assembled transcripts from long to short by length, accumulating the length of the transcript to 50% of the total length, corresponding to the length of the transcript, and so on. Figure 1(**A**) Assembly result sequence length distribution map of transcripts and unigenes in Xuzi-8 sweetpotato cultivar under drought stress conditions. The horizontal axis represents the length intervals of the transcripts and unigenes, and the vertical axis represents the number of transcripts and unigenes. (**B**) Species distribution of the top BlastX matches of the transcriptome unigenes in the non-redundant protein database (Nr) data base under drought stress conditions.

The species which provided the best BLASTX matches was *Nicotiana tomentosiformis* (19.86%), with more than 8,000 genes showing high homology, followed by *Nicotiana sylvestris* (19.83%), *Solanum tuberosum* (8.72%), *Sesamum indicum* (7.78%), *Solanum lycopersicum* (4.93%), and *Solanum pennellii* (4.51%). The reference sequence genome for Ipomoea batatas was ranked at the 8th position with only 3.45% sequence similarity. These data represent a sequence resource for genetic and genomic studies of sweet potato that will broaden the understanding of drought resistance mechanisms (Fig. [1](#Fig1){ref-type="fig"}B).

The assembled unigenes were validated and annotated using BLASTX search based on the sequence similarity search against the NCBI non-redundant (NR) protein database, Clusters of Orthologous Groups (COG) database, Gene Ontology (GO) database, Kyoto encyclopedia of genes and genomes (KEGG) database, Eukaryotic orthologous groups (KOG) database, protein family (Pfam) database and Swiss-Prot database. In total, 16,891 unigenes were detected which had good comparability with known gene sequences in at least one database, corresponding to 96.47% of the total unigenes. A total of 7,603 (45.01%), 13,860 (82.06%), 6,482 (38.38%), 10,740 (63.58%), 14,646 (86.71%), 14,620 (86.55%) and 16,866 (99.85%) unigenes were annotated in COG, GO, KEGG, KOG, Pfam, Swiss-Prot and Nr, respectively (Table [3](#Tab3){ref-type="table"}).Table 3Statistics of unigenes annotated in public database.Annotated databaseAnnotated number value (%)300 ≤ length \< 1,000 value (%)Length ≥ 1,000 value (%)COG7,603 (45.01)984 (05.83)6,619 (39.19)GO13,860 (82.06)2,040 (12.08)11,817 (69.96)KEGG6,482 (38.38)949 (05.62)5,533 (32.76)KOG10,740 (63.58)1,327 (07.86)9,413 (55.73)Pfam14,646 (86.71)2,009 (11.89)12,637 (74.81)Swiss-Prot14,620 (86.55)1,986 (11.76)12,631 (74.78)Nr16,866 (99.85)2,474 (14.65)14,389 (85.19)All16,8912,49314,395

Comparison of differentially expressed genes among four libraries between treated versus control samples {#Sec12}
--------------------------------------------------------------------------------------------------------

In order to quantify gene expression, the expression of each unigene was calculated using FPKM values (obtained fragments per kilo base of transcript per million mapped reads). Differentially expressed genes (DEGs) were identified by performing pairwise comparisons among the samples. Thresholds of false discovery rate (FDR) ≤ 0.05 and \|log~2~ ratio\|≥ 1 were set to identify significantly different expression between treated and untreated plants at each time point^[@CR32]^.

A rigorous analysis of the algorithm was performed at FDR 0.01 and 0.05 to determine the differences between the two levels. At 1, 6, 12, and 48 h of stress exposure, compared to the no-treatment control, 2, 724, 1,381 and 2,177 significantly up-regulated unigenes and 3, 545, 1811 and 3,132 significantly down-regulated unigenes were detected, respectively at the transcriptional level, with values of FDR ≤ 0.001. Overall, 414 up-regulated unigenes and 293 down-regulated unigenes were common at all three time points (6, 12, and 48 h) (Fig. [2](#Fig2){ref-type="fig"}). At 0.05 false discovery rate (FDR), the number of 2, 1,041, 1,720 and 2,583 significantly up-regulated and 3, 1,064, 2,366 and 2,366 significantly down-regulated were observed at 1, 6, 12 and 48 h after exposure to drought stress, respectively. Furthermore, when exposed to drought for 6, 12 and 48 h, 624 unigenes were found to be up-regulated, while 626 unigenes were down-regulated at 0.05 FDR (Fig. [2](#Fig2){ref-type="fig"}A). Moreover, the comparison of FPKM values from the different drought treatment time-point libraries versus control samples for each time point was more accurate than comparison with control values at 0 time point. Ultimately, in the present study, a pairwise comparison algorithm of FPKM was used for comparison of control libraries at each sampling point (FDR at 0.05 and \|Log~2~ FC ≥ 1\|).Figure 2Statistical chart of differentially expressed genes (DEGs) transcriptome in response to drought stress. Number of DEGs (up- and down- regulated) of four libraries including 1, 6, 12 and 48 h of salt stress treatment as compared to control at two different levels of false discovery rate (FDR) including 0.001 and 0.05 (**A**). Venn diagram analysis of all induced unigenes comparing each treatment with the control at the same time point including C1 versus D1, C4 versus D4, C5 versus D5 and C6 versus D6 at 1, 6, 12 and 48 h of drought stress (**B**). Venn diagram analysis of all induced unigenes comparing all treatments with the zero time point one control against all treatments C0 versus D1, C0 versus D4, C0 versus D0 and C0 versus D6 at 1, 6, 12 and 48 h of drought stress (**C**).

There were only 2 up-regulated unigenes TFs with significant expression levels. On the other hand, there were three significantly down-regulated unigenes, including thylakoid lumenal protein, Thermo-spermine synthase, and uncharacterized protein sequence after 1 h of drought stress. Upon exposure to drought stress for 6, 12 and 48 h, the number of commonly up-regulated unigenes was 624. Among these unigenes, 13 unigenes were annotated to an uncharacterized protein sequence. On the other hand, the total number of down-regulated unigenes was 626 at three different time points. Among 293 unigenes, 15 unigenes were annotated to an uncharacterized protein sequence (Fig. [2](#Fig2){ref-type="fig"}A, B).

After 6 h of drought stress, 1,041 up-regulated genes were differentially expressed in the leaf tissues. Among the up-regulated unigenes, ATP-dependent CLP protease was the most highly expressed gene (eightfold higher than control), followed by homeobox- leucine zipper protein (ATHB-7) and EID1-like F-box protein-3 which are expressed fivefold higher in PEG-treated samples than in control treatment and were. Furthermore, there were seven unigenes involved in drought tolerance, which were classified into different families, including GID1-like gibberellin receptor, LanC-like protein (GCL2), cinnamoyl-Co-A reductase-1, 2-hydroxy isoflavone dehydratase and carboxylesterase-1 (CXE-1). In addition, there were 10 unigenes involved in biological processes related to drought stress response, which belonged to four different families, including abscisic acid insensitive-5 like protein (bZIP-8), EID1-like F-box protein-3, 9-cis-epoxy carotenoid dehydrogenase, and cysteine protease inhibitor. On the other hand, there were 1,064 down-regulated unigenes in plants under PEG-induced drought stress compared to the control plants. The most pronounced change in gene expression occurred in down-regulated unigene was aligned to trypsin and protease inhibitor that exhibited a ninefold reduction in PEG-treated plants compared to control, followed by the unigene Methyl-CPG binding domain-11, MBD-11 which had a fivefold reduction in PEG-treated samples compared to control. In addition, two down-regulated unigenes were aligned to the defense mechanisms cluster, encoding p*rotein tyrosine* phosphatases and carboxylesterase-8. Furthermore, the *expression levels* of aspartic protease in guard cell and Jasmonate O-transferase were significantly down-regulated which contribute to processes involved in drought stress responses, including, respectively (Supplementary Table [S1](#MOESM1){ref-type="media"}).

After 12 h of drought stress, there were 1,720 up-regulated unigenes and 2,366 down-regulated unigenes in PEG-treated plants compared to control plants. Among the significantly up-regulated unigenes, 1, 5 and 3 unigenes were expressed eightfold, sixfold, and fivefold higher in treated samples than in control, while the others were one--twofold higher than those of control. According to the information provided by the NR database, the highly expressed unigenes were annotated with the EID1-like F-box protein-3 and home-box leucine zipper protein HB-12 like. In addition, there were 17 unigenes associated with tolerance to drought stress, which encodes the gibberellin receptor (GID1-B), carboxylesterase, tetraketide alpha-pyrone reductase 1, lanC-like protein GCL2, 2-hydroxyisoflavanone dehydratase-like and 2-hydroxyisoflavanone dehydratase. Furthermore, 17 unigenes were involved in drought stress responses, which belonged to seven different families (*i.e*., abscisic acid insensitive like protein-5 (bZIP-8), heat shock protein-70 (Hsp-70), gibberellin receptor (GID1B), S-adenosyl-L-methionine-dependent methyltransferase, Xanthoxin dehydrogenase-like, cysteine protease inhibitor and EID1-like F-box protein-3). Out of 2,366 unigenes, 9 down-regulated were related to the category of defense mechanisms in which [l]{.smallcaps}-Galactono-1,4-lactone dehydrogenase, GID-1C, cinnamoyl Co-A reductase and carboxylesterase-8 play an important role (Supplementary Table [S1](#MOESM1){ref-type="media"}).

The largest number of genes was recorded at 48 h of drought stress among all time points. Out of 2,583 up-regulated unigenes, 2, 9, 6, 16, 27, 38 and 76 unigenes showed 11, 10, 9, 8, 7, 6 and 5-fold higher expression, respectively in treated samples than in control. The superior unigenes with expression levels over tenfold higher than those of control belonged to various families including trypsin and protease inhibitor, FAD-dependent urate hydroxyl-like, YT-521-B like domain, protein phosphatase, arogenate dehydrogenase, heat shock protein-70, phosphor-enolpyrovate carboxykinase, STAY-GREEN protein, and EID1-like F-box protein-3. Moreover, 16 genes were assigned categories to defense mechanisms against drought stress, which were annotated to GID1-like gibberellin receptor, LanC like protein (GCL1), carboxylesterase, protein phosphatase, tetraketide alpha-pyrone reductase 1, phosphoglucan phosphatase LSF2 (chloroplastic), and 2-hydroxyisoflavanone dehydratase. In addition, 16 unigenes belonging to seven different families (WRKY, apolipoprotein D-like, ABA insensitive like protein-5 (bZIP-8), FAD binding domain, 9-cis-epoxycarotenoid dioxygenase-1, protein farnesyltransferase, and EID1 like f-box protein-3) were associated with the drought stress response. On the other hand, 3,481 unigenes were down-regulated at 48 h of stress exposure. Out of 3,132 unigenes, 16 unigenes were down-regulated to levels ninefold lower than those of control.

GO and KOG classifications {#Sec13}
--------------------------

The *Gene Ontology* (*GO*) *classification* system was used to classify the possible functions of the obtained unigenes into 46 elements, which were assigned to three main categories, including biological process, molecular function, and cellular component (Fig. [3](#Fig3){ref-type="fig"}). Each unigene was at least assigned to one element of a category. In total, 38,524 unigenes were assigned to the category of biological processes which was further divided into 17 sub-categories. Among these sub-categories, "metabolic process", "cellular process", and "single organism process" included the largest numbers of 9,228, 8,452, and 6,777 unigenes, respectively. Furthermore, the sub-categories "biological regulation" and "response to stimulus" were ranked second with 3,390 and 3,083 unigenes, respectively. The "molecular function" category contained 16,956 unigenes, which were distributed in 16 different functional groups. The greatest number of unigenes (7,120 and 7,037) was classified into the sub-categories "catalytic activity" and "binding", respectively. Moreover, 24,504 unigenes were assigned to the "cellular component" category and were further divided into 13 different sub-categories. For the "cellular component" category, the top largest categories were "cell part" (8,076 unigenes), "membrane part" (3,928), "organelle" (3,491), and "organelle part" (3,238). Additionally, a total of 286 unigenes were identified to be involved in drought stress tolerance response biological process according to Gene Ontology (GO) categories.Figure 3Gene ontology (GO) classifications in sweetpotato (Xuzi-8 cultivar), the percentage indicate the proportion of obtained unigenes with the GO annotations.

The COG database was used to classify the potential functions of unigenes (Fig. [4](#Fig4){ref-type="fig"}A). In total, 10,740 unigenes (63.58%) were annotated and classified into 25 functional classifications. The results of the functional annotation analysis showed that the majority of unigenes (1608, 13.63%) were associated with "posttranslational modification", "protein turnover", and "chaperones" categories, followed by the "signal transduction mechanisms" category which contained 1,000 unigenes (8.85%). Furthermore, 102 unigenes were found to be involved in mechanisms of drought tolerance in sweet potato. In addition, 1962 (16.63%) and 493 (4.18%) unigenes with uncharacterized function or poorly characterized belonged to the "general function prediction only" and "function unknown" categories, respectively. On the other hand, no genes were found to be involved in "cell motility" category.Figure 4(**A**) Clusters of orthologous groups (COG) classification in Xuzi-8 sweetpotato cultivar for unigenes obtained under drought stress conditions. Genes from the same Orthologous have the same function, so that direct functional annotations to other members of the same KOG cluster. (**B**) the most enriched KEGG clusters in Xuzi-8 sweetpotato cultivar under drought stress conditions. The most enriched 22 clusters out of 123 clusters were presented in this figure.

KEGG pathway mapping {#Sec14}
--------------------

The KEGG pathway analysis showed that a total of 16,891 unigenes and 6,482 sequences were annotated and assigned to 123 pathways. At the first hour of PEG-induced drought stress only Ascorbate and Aldarate metabolism KEGG pathway (Ko00053) contained one up-regulated unigene encoded myo-Inositol oxygenase as compared to control. A number of 82, 104 and 117 pathways included differentially expressed genes as compared to control. The two largest pathway groups were metabolic pathways (ko01100, 1829 unigenes) and biosynthesis of secondary metabolites (ko01110, 943 unigenes). Moreover, the pathways such as "Carbon metabolism (ko01200)", "Ribosome (ko03010)", "Biosynthesis of amino acids (ko01230)", "Plant hormone signal transduction (ko04075)", "Photosynthesis (ko00195)" and "Protein processing in endoplasmic reticulum (ko04141)" contained 390, 359, 327, 270, 236, and 263 unigenes, respectively (Fig. [4](#Fig4){ref-type="fig"}B). The KEGG pathway enrichment analysis provides insights into the mechanisms involved in the drought tolerance of Zi-8 sweet potato cultivar.

Differential expression of transcription factors (TFs) regulated by PEG treatment {#Sec15}
---------------------------------------------------------------------------------

Among the 16,891 DEGs differentially regulated by PEG treatment, 5, 1754, 3,709 and 5,755 encode TFs at 1, 6, 12 and 48 h of PEG treatment, respectively, belonging to 60 different families. The most abundant transcription factor family found during all time points under drought stress was zinc finger proteins (ZFPs) (212 unigenes), followed by MYB (147 unigenes), ERF (88 unigenes), TSA (73 unigenes), bZIP (63 unigenes), WRKY (33 unigenes), NAC (30 unigenes), and bHLH (5 unigenes) (Fig. [5](#Fig5){ref-type="fig"}A, Supplementary File [1](#MOESM2){ref-type="media"}).Figure 5Transcription factors differentially expressed (**A**) and stress related protein genes (**B**) induced under drought stress in Xuzi-8 sweetpotato cultivar.

After 1 h of PEG treatment, only the beta-amylase 1(BRI1-EMS-SUPPRESSOR1 (BES1) TFs family) and inositol oxygenase 2-like (Mediator of RNA polymerase II Transcription, Subunit 6 (MED6) TFs family) genes were up-regulated. BES1TFs family plant-specific transcription factors that regulate BR-responsive genes. MED6 TFs family among its related pathways are regulation of lipid metabolism and plant growth.

At 6 h of PEG-induced drought stress, there were 1754 unigenes encoded as TFs including 920 up-regulated and 834 down-regulated. These up-regulated TFs were classified into 61 different families including AP2/ERF (48 unigenes) followed by NAC (31), bZIP (27), DBP (26), MYB (25), HB-HD-ZIP (20), WRKY (19), B3-\>B3 (19), C2C2-\>C2C2-GATA (19), C3H (19) and bHLH (15). Furthermore, the most pronounced TFs families at this time point included DBP (3), HB-\>HB-HD-ZIP (2), NF-Y-\>NF-YC (1) and SNF2 (4) giving more than 4 folds of expression level higher than control. Among all TFs the highest one was expressed eight fold higher than control was aligned to ATP-dependent Clp protease ATP-binding subunit CLPT1 (chloroplastic). On the other hand, among the 834 down-regulated TFs, the lowest expressed as compared to control were aligned to Nitrate reductase 2 (C3H), Germin-like protein (SET), 14 kDa proline-rich protein (Trihelix), glucan endo-1,3-beta-glucosidase 12 (AP2/ERF-ERF), Sporamin B and Methyl-CpG-binding domain-containing protein 11 (MBD) (Fig. [5](#Fig5){ref-type="fig"}A, Supplementary File [1](#MOESM2){ref-type="media"}).

After 12 h of PEG-induced drought stress there were 3,709 unigenes encoded TFs including 1578 up-regulated and 2,131 down-regulated. The up-regulated TFs were classified into 81 families including same families present at 6 h beside 20 TF families more. The highest number of genes was included under MYB (63), AP2/ERF-ERF (53), NAC (52), WRKY (42), bZIP (33), HB-HD-ZIP (32), C3H (31), C2C2-GATA (30) and C2H2 (29). The highest expressed TFs were annotated to C2C2-GATA (14), DBP (1), HB-BELL (1), HB-HD-ZIP (3), LOB (2), MADS-MIKC (1) and SNF2 (1) families. Among all, the highest TFs was aligned to Protein phosphatase 2C (DBP TF family) which was expressed nine folds higher than control. On the other hand, among the down-regulated TFs the lowest expression were given by gibberellin induced protein (mTERF), cytochrome P450 (Trihelix), glutamine synthase 2,Tubulin beta-2 chain (SWI/SNF-BAF60b), bHLH63-like (bHLH), chlorophyll a-b binding protein of LHCII type I-chloroplastic (GNAT), magnesium-chelatase subunit ChlH-chloroplastic (Jumonji) and 21 kDa protein-like isoform X2 (Fig. [5](#Fig5){ref-type="fig"}A and Supplementary File [1](#MOESM2){ref-type="media"}).

At 48 h, the plant gave the maximum response to drought stress giving 5,755 significantly expressed TFs including 2,418 up-regulated and 3,337 down-regulated. The number of active TFs families at 48 h was similar to 12 h of drought. However the highest number of genes were included under C3H (70), MYB (68), AP2/ERF-ERF (63), NAC (61), WRKY (50), C2H2 (44), bZIP (41), B3-\>B3 (35), bHLH (35), PHD (35), SNF2 (34), ARID (33), HB-HD-ZIP (32) and FAR1 (31). Among all TFs families the highest unigenes were belonging to B3-\>B3-ARF, C2C2-GATA, C3H, DBP, LIM, MADS-MIKC and SNF2 reaching more than 7 folds higher than control. On the other hand, among all down-regulated TFs the lowest TFs as compared to control were classified as B3-ARF, C2C2-CO-like, GNAT, HSF, LUG, SET and Trihelix TF families (Fig. [5](#Fig5){ref-type="fig"}A and Supplementary File [1](#MOESM2){ref-type="media"}).

Protein kinases differentially expressed in response to drought stress {#Sec16}
----------------------------------------------------------------------

Out of 376 protein kinases (PKs) unigenes which were differentially expressed in leaf tissues under drought stress**,** 45 unigenes encoded receptor-like protein kinase (RLKs) in different libraries of control and drought stress (Table [4](#Tab4){ref-type="table"}). Cysteine-rich receptor-like protein kinase and receptor like protein kinase HAIKU-2 genes were significantly up-regulated after 6, 12, and 48 h of drought stress. Besides these two genes, G-type lectin S-receptor-like protein kinase gene was up-regulated after 6 h and putative receptor-like protein kinase gene was up-regulated after 12 and 48 h of drought stress. On the other hand, the proline-rich receptor-like protein kinase (PERK-15) gene was down-regulated at all three time points (6, 12, and 48 h) of drought stress. Moreover, after 48 h of PEG treatment, the number of down-regulated RLKs unigenes was 10, which were then aligned to glycerophosphodiester phosphodiesterase protein kinase domain-containing (GDPDL2), carbohydrate-binding protein of the ER, leucine-rich repeat receptor-like protein kinase, receptor-like protein kinase 5, and receptor-like protein kinase HSL1.Table 4Classification of drought-inducible regulated proteins into different categories according to sequence similarity with previously known genes.Regulatory protein categoryNo of total unigenesC1 versus D1C4 versus D4C5 versus D5C6 versus D6Up-regulatedDown-regulatedUp-regulatedDown-regulatedUp-regulatedDown-regulatedUp-regulatedDown-regulated**Protein kinases and phosphatases**Histidine Kinase1200040505Calcium-dependent protein kinase900001010Mitogen-activated protein kinase26000002010Mitogen-activated protein kinase kinase1500000412Mitogen-activated protein kinase kinase kinase600002000Serine/threonine-protein kinase2390018840173237**Protein phosphatase**SNF1-related protein kinase1400000010LRR receptor-like kinase2800112238Receptor-like kinase11200941591925Protein phosphatase167002713633517Serine/threonine-protein phosphatase4800003023Phospholipase30003473116Calmodulin1600406050Calcineurin B-like proteins1200301202CBL-interacting protein kinases (CIPK)370072170123**Plant hormones**Abscisic acid2180029748124529Aldehyde oxidase500000030Zeaxanthin epoxidase200000001Molybdenum cofactor sulfurase3000000019-Cis-epoxycarotenoid dioxygenase700506060Abscisic acid receptor2100919565**Ethylene**1-Aminocyclopropane-1-carboxylate oxidase11000010051-Aminocyclopropane-1-carboxylate synthase100000001Ethylene receptor200000000Ethylene insensitive3100103020Ethylene response factor2000101811EIN3-binding F-box protein1300202017**Jasmonic acid**Lipoxygenase (LOX)1300001014Allene oxide synthase100000010Jasmonate O-methyltransferase300010002**Salicylic acid**Salicylic acid-binding protein15800971948976**Auxin**Tryptophan aminotransferase-related protein200000102Auxin transporter-like protein400010321Auxin efflux carrier component100000001Auxin response factor500010101

The other PKs expressed under drought stress belonged to the SNF1-related protein kinase (11 unigenes, 1 up-regulated at 48 h), phosphatidylinositol 4-phosphate 5-kinase (8 unigenes, 1 up-regulated at 48 h), cyclin-dependent kinase (17 unigenes, 3 down-regulated at 48 h), CBL-interacting protein kinase 2-like (9 unigenes, 2 down-regulated at 6 h, 2 up-regulated at 12 h, 3 up- and 3 down-regulated at 48 h), histidine kinase (12 unigenes, 5 down-regulated at 6, 12, and 48 h) and mitogen-activated protein kinase kinase kinase YODA (47 unigenes, 6, and 12 down-regulated at 12, and 48 h) (Table [4](#Tab4){ref-type="table"}).

Stress-related protein genes induced by drought {#Sec17}
-----------------------------------------------

### Genes encoding redox regulation-related proteins {#Sec18}

A set of 189 unigenes encoding redox-regulated proteins was identified under drought stress in leaf tissues of Xuzi-8 sweet potato cultivar. These genes encoded ascorbate peroxidase (APX), peroxidase (POD), oxidoreductase, monodehydroascorbate-reductase (MAD), glutathione S-transferase, Acyl Co-A reductase and Acyl Co-A syntheses^[@CR33]^. All of these unigenes were commonly expressed in both PEG treated and untreated plants at 1 h after drought stress. The unigenes encoding ascorbate peroxidase and Acyl CO-A reductase were down-regulated after 6, 12, and 48 h of PEG treatment. The activities of bifunctional monodehydroascorbate reductase, carbonic anhydrase nectarin-3-like and monodehydroascorbate reductase were up-regulated after 12 and 48 h of PEG-induced drought stress. In addition, 2, 1, and 10 unigenes encoding Acyl Co-A syntheses (3-ketoacyl-CoA synthase 6, 3-ketoacyl-CoA synthase 1, 3-ketoacyl-CoA synthase 19, Long-chain acyl-CoA syntheses 8, 3-ketoacyl-CoA synthase 2 and Long-chain acyl-CoA syntheses 8) were up-regulated after 6, 12, and 48 h of drought stress, respectively, while the unigenes encoding 3-ketoacyl-CoA synthase 4 were down-regulated. Furthermore, 51 unigenes encoding glutathione S-transferase were down-regulated in leaf tissues under drought stress; however, 14 unigenes which encoded glutathione S-transferase zeta class like and T1-like were up-regulated (Fig. [5](#Fig5){ref-type="fig"}B and Supplementary File [1](#MOESM2){ref-type="media"}).

### Genes involved in carbohydrate metabolism and osmotic adjustment {#Sec19}

Thirty strongly up-regulated unigenes involved in carbohydrate metabolism and osmotic adjustment were constantly expressed in leaf tissues at 6, 12, and 48 h under drought stress conditions. These genes encoded beta-glucosidase-like SFR2, cinnamate beta-D-glucosyltransferase, UDP-glucuronic acid decarboxylase 2, 7-deoxyloganetin glucosyltransferase, fumarylacetoacetase, Alpha-amylase, glucose-6-phosphate 1-epimerase, Galactinol synthase 2, scopoletin glucosyltransferase, alpha, alpha-trehalose-phosphate synthase \[UDP-forming\], and histidine phosphatase superfamily. The expression level of the top up-regulated genes encoding galactinol synthase 2 and UDP-glycosyltransferase 85A5-like was sevenfold higher in treated than in control plants after 48 h of drought stress. On the other hand, 59 unigenes annotated as the "carbohydrate transport and metabolism" category were down-regulated at 6, 12, and 48 h of drought stress. These down-regulated unigenes encoded aquaporin (TIP1-1, TIP2-1, PIP2-7, TIP1-1, and NIP2-1), beta-galactosidase, glyceraldehyde-3-phosphate dehydrogenase B, (chloroplastic) phosphoglucomutase, chloroplastic, phosphoglycerate kinase (chloroplastic), protein plastid transcriptionally active 16 (chloroplastic), sedoheptulose-1, 7-bisphosphatase (chloroplastic), and UDP-D-apiose/UDP-D-xylose synthase 2. Moreover, the down-regulated unigenes with the lowest expression level (sevenfold lower in PEG treated than in control plant) encoded Sedoheptulose-1, 7-bisphosphatase, chloroplastic, and aquaporin PIP2-7 (Fig. [5](#Fig5){ref-type="fig"}B, Table [5](#Tab5){ref-type="table"} and Supplementary File [1](#MOESM2){ref-type="media"}).Table 5Classification of drought-inducible functional proteins into different categories according to sequence similarity with previously known genes.Protein functions categoryNo of total unigenesC1 versus D1C4 versus D4C5 versus D5C6 versus D6Up-regulatedDown-regulatedUp-regulatedDown-regulatedUp-regulatedDown-regulatedUp-regulatedDown-regulated**DEGs involved in osmotic adjustment**Arginine decarboxylase1100107010Choline monooxygenase100000001Pyrroline-5-carboxylate synthase200000010S-adenosylmethionine decarboxylase1700015002Trehalose-phosphate phosphatase100001000Galactinol synthase300303030**DEGs involved in detoxification**Respiratory burst oxidase100000001superoxide dismutase2200010205Catalase48000116148Peroxidase94001156192231Glutathione reductase300002020Dehydroascorbate reductase1500006383Glutathione peroxidase700002010Glutathione S-transferase620070882114Ascorbate peroxidase1500031408Metallothionein1100200040**DEGs involved in osmotic protection, transporters and channels**Aquaporin5800019229335ABC transporter family protein38200284361547284Multidrug resistance protein ABC-transporter protein2100235256Chloride channel900003101Betaine aldehyde dehydrogenase400000020Delta-1-pyrroline-5-carboxylate synthetase100000010Proline oxidase/dehydrogenase400400010**LEA proteins**Late embryogenesis abundant protein1200233313Dehydrin1800140160150**Others**Pathogenesis-related protein700005070Drought-related protein1200101000Dehydration-related protein200000022Heat shock protein60002018280

Validation of differential gene expression {#Sec220}
------------------------------------------

In order to validate differential gene expression results obtained by transcriptome sequencing, the expression levels of six randomly selected DEGs under PEG-induced drought stress were evaluated by qRT-PCR. Specific primer pairs for selected unigenes were designed as shown in Table [6](#Tab6){ref-type="table"}. The expression levels determined by qPCR followed the same trend with the changes in the results of transcriptome sequencing, which suggests that our transcriptome sequencing data are highly reliable (Fig. [6](#Fig6){ref-type="fig"}).Table 6Forward and reverse primers used for QPCR validation experiment.Unigene codeGene nameForward (5′-\>3′)Reverse (3′-\>5′)1g14139Protein disulphide-isomerase (PDI)TTTTCGTAATCCTGGCGGCGACCGTGTCGGAGAAGTTGGT2g2083Sugar carrier protein (STC)CCGCCCGAGAGGAAGTACAGGGGCCAGACTCGTCGTACTC3g16087Nucleoredoxin 2GGAGGCGTTTGCTGCTTACCTGCCCTCGACGTTGAACCTT4g2298Thioredoxin-likeTTCGACGCCTACTTCGGGACGTCCGGCCCGAGAATTACCA5g8096Thioredoxin superfamily protein isoform 1TCTAGCAATGGTGCTGCGGAAGCCCTACTGCTCCAACTACG6g6453Zinc finger protein CONSTANS-LIKE 5 (COL5)TATGGTGTCGGTACGGCGTCTGACTCAGTTCTAGCGGCCC Figure 6Quantitative real-time PCR analysis and fragments per kilo base million mapped reads (FPKM) of for 6 randomly chosen differentially expressed genes (DEGs). A, B, C, D, E and F; expression level determined by qPCR for the unigenes g14139, g2083, g16087, g2298, g8096 and g6453, respectively.

Discussion {#Sec21}
==========

Sweet potatoes are a rich source of nutrients. The production of sweet potato is limited by several abiotic stresses, resulting in tremendous yield losses^[@CR6]^. Among these abiotic stresses, drought is a major threat to the sustainable production of sweet potatoes. Due to the complex nature of the sweet potato genome, extensive genomic, transcriptomic, and proteomic approaches are required to better understand its genome^[@CR34],[@CR35]^. The lack of a reference genome is also a barrier to the accurate identification of candidate genes involved in drought tolerance. The present study was conducted to identify drought stress-responsive genes using next and third-generation sequencing techniques considering the severity of drought stress, the complex nature of the sweet potato genome, and the absence of a reference genome. Dual RNA-sequencing analysis provides molecular insights into defense mechanisms in plants against drought stress, which may further facilitate breeding for drought tolerant sweet potato cultivars^[@CR32]^.

In the present study, HiSeq Illumina sequencing was used to characterize the transcriptome profiles of purple-fleshed sweet potato cultivar Xuzi-8 in response to different durations of drought stress. The second and third-generation sequencing methods produced a total of 184,259,679 clean reads, which were assembled into 17,508 unigenes with an average length of 1,783 base pairs (bp). According to GO, the top cellular component involving the largest number of genes was "cell part" sub-category. The balance between the metabolic and catabolic activities estimates the degree of resistance of a plant to drought stress at each time point^[@CR36]^. Our results are in agreement with those of a previous study, reporting that the top three GO terms were metabolic process, cellular process, and catalytic activity in Ipomoea trifida under drought stress^[@CR10]^.

The results of COG function showed that the greatest numbers of genes were assigned to "posttranslational modification", "protein turnover", "chaperones" and "signal transduction mechanisms" categories^[@CR37]^. When plants are exposed to drought stress, the signal transduction mechanisms, which are responsible for detecting a physical or chemical signal transmitted through the cell as a series of molecular events, become activated, which ultimately results in a cellular response^[@CR38],[@CR39]^. After the detection of signals, the plant modulates the activities by enzymatic modification of proteins following protein biosynthesis and protein turnover through achieving the balance between protein synthesis and degradation. In addition, chaperone proteins, which assist in protein folding as a chemical reaction under stress conditions, act to prevent or correct damage caused by miss folding^[@CR40]^.

Drought inducible genes have recently been classified into two categories, namely, functional and regulatory proteins. Functional proteins are involved in abiotic stress tolerance, i.e., LEA, chaperons and dehydrin, detoxification, osmolytes biosynthetic (proline and sugar), protease, and water channel transporters proteins. Regulatory proteins are involved in regulating other genes that participate in drought tolerance, including TFs (bZIP, DREB2, MYC, AREB, MYB, and NAC), PKs, phosphatase, lipid metabolism, and ABA biosynthetic involved genes^[@CR40]^.

Differentially expressed genes among the four PEG treated libraries as compared to control library {#Sec22}
--------------------------------------------------------------------------------------------------

In order to establish an accurate base for DEGs comparison, different significant levels were used (0.05 and 0.01). Furthermore, the use of the same control library for comparing with the treated libraries did not give representative results. Therefore, finally, a pairwise comparison algorithm of FPKM was used for comparison of control libraries of each sampling point comparing with the treated libraries were chosen at *P* value 0.05.

In the current study, at the first hour of PEG-induced drought stress, there were only three down-regulated unigenes in the leaf tissues which were aligned to thylakoid luminal protein and thermosporamine synthesis proteins. These genes were found to be involved in cell tip growth, meristem initiation, elongation, and cell wall organization during the first hour of drought stress. This explains that in sweet potato at the first hour of drought stress, the plant starts hormonal signaling, and growth diminishes^[@CR41],[@CR42]^.

A set of gene families, including GID1-like gibberellin receptor (GID1), 2-hydroxy iso-flavanone dehydratase (HFD), Carboxylesterase-1 (CXE-1), Cinnamoyl-Co-A reductase-1 (CCo-A-R), and LanC-like protein (GCL2), were continued to be upregulated at high levels in Xuzi-8 leaf tissues at the different time points of PEG-induced drought stress (6, 12, 48 h). These genes play a vital role in plant defense mechanisms against drought stress^[@CR43]^. Besides their role in drought tolerance, the four genes, including GID1, HFD, CCo-A-R, and CXE-1, control the hydrolase activity, which is related to energy metabolism and storage. Moreover, hydrolase activity is a chemical process in which a molecule of water is added to a substance that decreases the water loss from plant^[@CR44]^. Thus, according to the current results of Xuzi-8, the plant resorted to reducing water loss through activating the genes of hydrolases activity in leaf tissues. LanC-like protein (GCL2) was reported previously to be one of the main drought-tolerant genes in Arabidopsis^[@CR45]^ wheat^[@CR46]^, and peach^[@CR47]^, which explains the current results.

According to the current results, at 6 h of PEG- induced drought stress, leaf tissues start to sense the stress through the early perception of water deficiency leading to activation of the stress-responsive genes and stimulation of stomatal closure to reduce water loss^[@CR48]^. The highest expressed gene in leaves at 6 h of drought stress was ATP-dependent Clp proteases (eightfolds higher than control). This gene minimizes the protein damage and maintains the quality of cellular proteins, thereby increases the tolerance for drought stress^[@CR49],[@CR50]^. The current results indicated that the EID1-like protein 3 was expressed five-fold higher in leaf tissues than that control and continued with a high level at 6, 12, and 48 h of PEG-induced drought stress. This gene is involved in regulating the responses controlled by plant hormones, which not only regulate many aspects of plant development but also integrate responses toward drought stress^[@CR51]^. On the other hand, the current results found that several genes were significantly down-regulated as compared to control. Trypsin-protease inhibitor showed a nine-fold reduction in Xuzi-8 leaf tissues as compared to control, which is in agreement with a previous study conducted on tomato drought resistant cultivar^[@CR52]^. Protein aspartic protease down-regulated targets in ABA-signaling increasing drought avoidance through guard cells closure^[@CR53]^.

In the current study, the highest number of drought stress-related genes was active at 12 h of PEG-drought stress. This means that Xuzi-8 sweet potato cultivar is at the peak of restoring to resist the stress after 12 h of its onset. In total, two gene families were acting as the superior genes at 12 h of drought stress giving an expression level that was more than eight folds than that in the control. These sets of genes were annotated to EID1-like F-box protein-3, which continued at the highest level from 6 h till 48 h of exposure to PEG-drought stress.

After 48 h, the expression of the fast action and non-specific drought stress genes fall back to the normal level, and only the drought specific genes were active. Remarkably, trypsin and protease inhibitor, FAD-dependent urate hydroxyl-like, YT-521-B like domain, protein phosphatase, arogenate dehydrogenase, heat shock protein-70, phosphor-enolpyruvate carboxykinase, STAY-GREEN protein, and EID1-like F-box protein-3 genes were expressed at ten-fold higher level than control at 48 h of PEG-induced drought stress. Trypsin and protease inhibitors were expressed at tenfold higher level than the control at 48 h of PEG-induced drought stress. However, they were expressed at nine-fold lower level than that of the control at 12 h. Trypsin protease inhibitor potentially increases drought stress tolerance and inhibits pathogen infection; therefore, it can be used as an excellent candidate as a lead compound for the development of a drought-tolerant character^[@CR54]^. FAD-dependent urate hydroxyl-like was found to be involved in the catalysis of purine, which plays a role in the recycle of nitrogen for remobilization to support new growth and reproduction^[@CR55]^. YT-521-B like domain is involved in cell growth, differentiation, and performs an important function in mediating a wide range of cellular responses under multiple abiotic stresses^[@CR56]^. Arogenate dehydrogenase is responsible for root elongation at low water potential^[@CR57]^. Heat shock protein-70 regulates the abscisic acid-induced antioxidant response under drought and heat stress^[@CR58]^. Phospho-enolpyruvate carboxykinase was found to be involved in acclimation of the plant to drought stress^[@CR59]^. STAY-GREEN protein encodes the regulatory proteins trigging chlorophyll catabolism and known to delay senescence, and improve plant survival and yield^[@CR60]^.

Transcription factors (TFs) regulated by PEG treatment {#Sec23}
------------------------------------------------------

In the current study, at the first hour of PEG-induced drought stress, only two TFs were up-regulated, including beta amylase-1, which is involved in the starch catabolic process and inositol oxygenase-2 like, which is involved in inositol catabolic process. These two TFs belong to BES1 and MED6 TF families. BES1 TFs family is a class of plant-specific transcription factors that cooperate with other transcription factors to regulate BR-induced genes. Moreover, it starts stress tolerance through the control of many basic cellular processes and stress responses, promotes vigor in plants, and prepares them to mount a dynamic response upon environmental challenges. MED6 TFs family among its related pathways is involved in the regulation of lipid metabolism and plant growth. The transcriptional changes in genes related to lipid metabolism indicate its function at the interface in response to drought stress^[@CR61]^. These two TFs, according to the current results, can be considered as early response genes in sweet potato leaf tissues for stress and may not be specific for drought stress. They were also found to be involved in the primary sensing of abiotic stress and in triggering and regulating cellular hormonal signaling^[@CR62]^.

The majority of drought-related unigenes in this study were classified under different families of transcription factors (TFs) such as AP2/ERF, NAC, bZIP, DBP, MYB, HB-HD-ZIP, WRK, B3-\>B3, C2C2-\>C2C2-GATA, C3H, and bHLH, which were identified at all studied time points except for the first time point of drought stress. It was reported that TFs are considered as early activated genes that are targeted by phosphatases and PKs^[@CR63]^. Among the active TFs, WRKY are largely up-regulated and known to be involved in drought stress defense pathways in plants^[@CR64]--[@CR66]^, and the members of AP2-EREBP TFs family are also shown to be highly up-regulated. They play an important role in drought stress signaling pathways induced by hormones^[@CR67]^. In addition, NAC factors were differentially expressed and played specific roles in response to drought stress^[@CR68]^.

The bHLH transcription factors (TFs) play a crucial role in controlling gene expression by regulating diverse biological processes such as growth, development, and drought stress response^[@CR69],[@CR70]^. In the current study, bZIP TFs were significantly up-regulated at all time points (6, 12, and 48 h of stress) and played crucial roles in response to drought stress and regulated numerous growth and developmental processes^[@CR71]^. These bZIP TFs included common regulatory factor (CPRF1), ABSCISIC ACID-INSENSITIVE 5-like protein 5, ABF2, bZIP28, and the transcription activator (TAF1). Interestingly, GO-annotation results showed that the aligned sequences of ABF2 are involved in response to water deprivation^[@CR72]^. C2H2-type zinc finger proteins belong to the major family of transcription factors that can be classified as osmotic stress-responsive genes^[@CR73]^ and bind to DNA, RNA, or protein^[@CR74]^. C3H transcription factor (TF) family included 8, 16, and 21 unigenes that were significantly up-regulated and assigned to "posttranslational modification", "protein turn over" and "chaperones" categories. These genes included C3HC4-type Ring-like zinc finger protein, XERICO^[@CR74]^, and SPAP32A, which have only been characterized in *Schizosaccharomyces pombe* species but till now have not been studied in sweet potato. MYB transcription factor family in the present study included two up-regulated genes *MYB44* (*MYB44*) and MYB1-R1, which are involved in posttranslational modification, protein turn over, and chaperone-like activities^[@CR75],[@CR76]^.

During the exposure to drought stress, the plant uses ABA as one of the key signaling molecules. In the current study, 9-cis-epoxycarotenoid dioxygenase was differentially expressed, which could be considered a key step in ABA synthesis^[@CR77]^. Further, many DEGs encoded protein phosphatase PP2C were differentially expressed and played an important role in ABA regulation^[@CR78]^. These results are in agreement with previous studies focused on drought stress^[@CR79]^. Ethylene responsive factors (ERFs) play a regulatory role in water logging responses, are also considered in ethylene sensor promotion, and reported to play a key role in cell senescence^[@CR10],[@CR80]^. They were down-regulated and contributed to the mitigation of drought stress in sweet potato leaves. In addition, the upregulated TF belongs to APETALA2/ethylene-responsive factor (AP2/ERF), which was observed to reduce transpiration and prevent a rapid decline in photosynthesis^[@CR18],[@CR81]^.

Protein kinases differentially expressed in response to drought stress {#Sec24}
----------------------------------------------------------------------

Remarkably, the DEGs encoded kinases, such as receptor-like kinases and LRR receptor-like kinase, were found to be largely up-regulated after 6, 12, and 48 h of PEG-induced drought stress. It is well established that LRR-RLKs play a valuable role in response to drought stress^[@CR82]^. On the other hand, a substantial number of DEGs were downregulated under drought stress as compared to the control at different time points. In the current study, another major class of differentially expressed regulatory proteins is PKs, which fincluded SNF1-related protein kinase, phosphatidylinositol 4-phosphate 5-kinase, cyclin-dependent kinase, CBL-interacting protein kinase 2-like, histidine kinase and mitogen-activated protein kinase kinase kinase YODA. PKs are considered drought stress sensor response genes that activate phosphorylation, which play important roles in response to drought stress^[@CR62],[@CR83]^.

Stress-related protein genes induced by drought stress in Xuzi-8 sweet potato cultivar {#Sec25}
--------------------------------------------------------------------------------------

In the current study, the second class are function proteins including osmolytes biosynthesis, water channel, sugar and proline transporters, detoxification enzymes, chaperones, and late embryogenesis abundant (LEA) proteins. This class of protein helps plants in controlling water efflux influx through various mechanisms during drought stress. Thus they protect the integrity of the cell membrane and provide an aqueous environment for plant organelles^[@CR84],[@CR85]^. In this study, many genes were identified to be involved in the metabolism of these osmolytes. Osmotically active solutes (osmolytes) are induced for the maintenance of cellular functions through biosynthesis of different substances as well as enzyme activity and proteins^[@CR86]^. In the current study, proline, trehalose, and glycine betaine were found to have a relation in combating drought stress. Another important osmolyte is proline, which accumulates under the first phase of stress. Proline is an essential component of drought tolerance, especially at the first phase. Its contribution is diverse as it participates in osmotic adjustment and stabilization of subcellular structure^[@CR87]^. In this study, two major proline components, including proline oxidase/dehydrogenase and Delta-1-pyrroline-5-carboxylate synthetase were found to be up-regulated at 6 h then decreased and upregulated again at 48 h of drought stress. It is required for maintaining growth at low water availability especially at the first phase of stress and to maintain growth at low water availability.

In the current study, it was found that a large number of DEGs encode proteins required for detoxification, such as protein detoxification-30, ascorbate peroxidase, glutathione peroxidase, glutathione S-transferase, ABC-2 transporter (Pleiotropic drug resistance, PDR), and superoxide dismutase. These genes have a major role in detoxification of ROS in the plant cell, which is extremely harmful as it causes lipid oxidation, DNA damage, and apoptosis. In agreement with our study, the role of protein detoxification has been documented in many plant species, i.e., *Ipomoea batatas*^[@CR18],[@CR88]^, Ipomoea pes-caprae^[@CR89]^*, Arabidopsis thaliana*^[@CR90]^*, Nicotiana tabacum*^[@CR91]^, *Oryza sativa* L.^[@CR92]^ and *Solanum lycopersicum*^[@CR93]^*.* In addition, TSA genes are involved in the detoxification of reactive oxygen species (ROS), reactive sulfur species (RSS), and reactive nitrogen species (RNS). TSA genes were significantly up-regulated, as they have a role in drought avoidance in plant^[@CR94]^.

Aquaporin encoded DEGs were identified in our study, including PIP2-1, SIP2-1 at 6, 12, 48 h of drought stress, which considered as a part of stress protecting protein families in sweet potato that can facilitate water uptake. Late embryogenesis proteins (LEAs), including late embryogenesis abundant protein and dehydrin, were highly up-regulated in the current experiment at 6, 12, and 48 h of drought stress. These LEAs proteins play a protective role under osmotic stress conditions, homeostasis of proteins and nucleic acids, stabilization of cell membranes, and maintenance of redox balance^[@CR95]^. In addition, a set of gene families including ABC transporter B, C, G family; LanC- like protein Gcl2 and MATE efflux family protein-5 play a vital role in a plant's defense mechanisms against drought stress.

Taken together, Xuzi-8 sweet potato plants start sensing drought stress after 1 h but the highest action taken against drought stress started after 6 h till and reached the maximum response at 48 h. There is a huge network of genes involved in drought stress tolerance mechanisms, which contribute in determining the level of tolerance of each cultivar. A substantial proportion of DEGs was identified in our study, which did not match with any database. Many of these DEGs should show a significant change in their expression levels as compared to control during the different time points of drought stress. Therefore, these DEGs are an unexplored resource that might be involved in responses to drought stress or may be enriched during drought tolerance in sweet potato. Thus our study of Xuzi-8 sweet potato transcriptome sequencing could be used for understanding drought tolerance mechanisms in sweet potato and further exemplify the function aspects of the sweet potato genome.

Supplementary information
=========================

 {#Sec26}

Supplementary file1 Supplementary file2

**Publisher\'s note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information
=========================

is available for this paper at 10.1038/s41598-020-69232-3.

This work was supported by the earmarked fund for China Agriculture Research System (CARS-10, Sweetpotato), the "333 talent project" fund of Jiangsu Province (BRA2016265), and China, Talented Young Scientist Program of Ministry of Science and Technology.

Li Qiang and Arisha designed the study. Arisha conducted the experiments, analyzed the results, performed the figures and wrote the manuscript. Li Qiang and Muhammad Qadir Ahmad revised the manuscript. Tang Wei, Kou Meng, Yan Hui, Liu Yaju, and Wang Xin, Zhang Yungang assisted in the experimental procedures.

The authors declare no competing interests.
