Insights into grapevine defense response against drought as revealed by biochemical, physiological and RNA-Seq analysis

Grapevine is an important and extensively grown fruit crop, which is severely hampered by drought worldwide. So, comprehending the impact of drought on grapevine genetic resources is necessary. In the present study, RNA-sequencing was executed using cDNA libraries constructed from both drought-stress and control plants. Results generated 12,451 differentially expressed genes (DEGs), out of which 8,021 genes were up-regulated, and 4,430 were down-regulated. Further physiological and biochemical investigations were also performed to validate the biological processes associated with the development of grapevine in response to drought stress. Results also revealed that decline in the rate of stomatal conductance, in turn, decrease the photosynthetic activity and CO2 assimilation in the grapevine leaves. Reactive oxygen species, including stress enzymes and their related proteins, and secondary metabolites were also activated in the present study. Likewise, various hormones also induced in response to drought stress. Overall, the present study concludes that these DEGs play both positive and negative roles in drought tolerance by regulating various biological pathways of grapevine. Nevertheless, our findings have provided valuable gene information for future studies of abiotic stress in grapevine and various other fruit crops.

Grapevine (Vitis vinifera L.) is an important crop, having 7.8 million hectares of cultivated land with an annual production of 67.6 million tons worldwide 1 . The climate change has noticeable effects on the survival and productivity of fruit plants 2 . Hence, the growth of grapevine is consequently, affected by abiotic stress, such as drought and salinity. Among these, drought has deleterious impacts on viticulture around the world 3,4 . Globally, 45% of the agricultural terrains are under constant/periodic water scarcities 5 , causing nearly 50% of yield losses. Plants as sessile organisms can make versatile vicissitudes in physiology and morphology that allow them to endure environmental stress. However, these adaptations are derisory to recover physiological water potential in the cell 6 . Plant response to these minimal water conditions is mediated by expression of numerous genes encoding stress-related proteins, enzymes and metabolites functioning in the various pathways of cell metabolism 7 . The genes induced by osmotic stress in plants are categorized into two groups, such as functional proteins and regulatory proteins 8,9 . In previous studies, several salient genes were identified in grapevine genome in response to biotic and abiotic stresses, 59 similar genes encoding putative WRKY proteins were identified from the available database 10 . Similarly, plant pathogenesis-related proteins were believed to be involved in plant defense and are usually induced during biotic and abiotic stresses 11 . In addition, over-expression of AtHSP70 has demonstrated its contribution in enhancing tolerance to abiotic stress 12 . Transcriptomic analysis of grapevine during drought stress is of vital importance whose results could provide a defense-related gene information, which offers a foundation for further development of drought-resistant grapevine cultivars.
Water scarcity is not the only threat to viticulture productivity, but also for wine quality 13,14 . Schultz proposed that an increase in environmental temperature due to rise in atmospheric CO 2 is a primary cause of water shortages for viticulture 15 . Grapevine possesses distinct molecular machinery which adjusts the circulation of water to leaf and then to the atmosphere by vessel anatomy 16 , stomatal conductance 17 and aquaporin 18 . Consequently, the sluggish leaf and shoot growth, elongation of tendrils, restrained internodes extension, leaf augmentation,

Results
The sequence data obtained from the Illumina deep-sequencing was submitted to Short Read Archive (SRA) database at NCBI under accession number SAMN04914490. After filtering, raw data yielded 42.47 and 53.05 million clean reads in control and drought-stressed leaf samples, respectively. The sequence alignment (soap2/ SOAPaligner; http://soap.genomics.org.cn) to the grapevine reference genome, allowed two base mismatch. The total mapped reads (73.44%) were corresponding to unique (72.01%) and multiple (1.44%) genomic positions (Supplementary Table S1).

Gene ontology (GO) and KEGG analysis of differentially-expressed genes. A total of 12,451
DEGs were subjected to GO and KEGG for annotation. All the DEGs were assigned to different groups based on their functional properties, such as biological process (21), cellular component (16) and molecular functions (14). Under the broad category of "Biological process", 4,537 transcripts were assigned to 'metabolic processes' (GO: 0008152), followed by 'cellular process' which consists of 3,632 transcripts (GO: 0009987) and 3,185 transcripts were involved in 'single organism process' (GO: 0044699). Similarly, under "Cellular component" category, 3,247 transcripts were assigned to 'cell' and 'cell part' (GO: 0005623, and GO: 0044464), followed by 2,358 transcripts in 'organelle' (GO: 0043226). Further, under "Molecular functions" category, 4,291 transcripts were in 'catalytic activities' (GO: 0003824), and 3,497 transcripts were assigned to 'binding activities' (GO: 0005488; Supplementary Table S3; Figure S1).
Several DEGs from the current study were subjected to KEGG annotation for further characterization of transcripts, where 12,451 transcripts were allotted to 306 KEGG pathways. Our results revealed that the highest numbers of transcripts (1,425) were involved in the "Metabolic pathway" (1,257 up-regulated, 168 down-regulated), followed by "Biosynthesis of secondary metabolites", which consists of 1,160 transcripts (681 up-regulated, 479 down-regulated), then 756 transcripts were recorded in "Plant-pathogen interaction pathway" (487 up-regulated, 269 down-regulated), while lowest transcripts (21) were recorded in "Sesquiterpenoid and triterpenoid biosynthesis pathway" in which 14 transcripts were up-regulated and 7 transcripts were down-regulated (Supplementary  Table S4).
Plant hormone signal transduction pathway under drought stress. The hormonal level, including auxin, was increased in drought treatment (1.626 ± 0.03 ng g −1 FW) compared to control (1.373 ± 0.02 ng g −1 FW). A similar trend was observed in abscisic acid that is 0.908 ± 0.01, and 0.257 ± 0.01 ng g −1 FW for drought and control treatments, respectively. In the same way, jasmonic acid in drought treatment sample was 1.67 ± 0.05 ng g −1 FW, whereas, in control, it was found to be 1.451 ± 0.03 ng g −1 FW. Further gibberellic acid (GA) in treated and control sample was recorded to be 1.671 ± 0.02, and 1.53 ± 0.02 ng g −1 FW, respectively. Alike brassinosteroid also showed 1.091 ± 0.01, and 1.073 ± 0.01 ng g −1 FW for drought and control treatment samples, respectively (Fig. 3). In grapevine transcriptome, several DEGs related to AUX, GA, ABA, JA, ET (ethylene), and BR were found in signal transduction pathways in drought stressed grapevine leaves. Under AUX signaling, three genes (down-regulated) related to auxin transport, eleven auxin-response factors (7 up-regulated and 4 down-regulated) involved in the transcriptional repressors were detected. Moreover, fifteen     Table S8).

Discussion
Drought stress can subdue the plant growth by hindering many physiological processes of plants. Chlorophylls (chls) are the principal light-absorbing pigments and key components of photosynthesis in plants. The physiological and transcriptomic studies of grapevine leaves responding to drought stress have revealed that chl contents were incredibly decreased, which in turn, inhibited the photosynthetic activity. Similarly, the reduction in chl contents has been reported in corn and chickpea in response to drought stress 21,22 . Moreover, transcriptomic data demonstrated that drought stress suppressed the chl biosynthesis process by inhibiting the activity of key enzymes, such as HemA (Glutamyl-tRNA reductase 1) and CHLH (Magnesium chelatase H subunit), which play an essential role in chla synthesis process 23 . Furthermore, in the chl cycle, the oxygenation reactions of chlorophyll(ide) a to chlorophyll(ide) b were catalyzed by chlorophyllide a oxygenase (CAO) 24 , whose activity likewise was decreased under drought stress, suggesting the inhibition of chl cycle. In contrast, the chlorophyll(ide) b to a conversion is catalyzed by chlorophyll(ide) b reductase NYC1 (CBR) and its activity was up-regulated, suggesting that chl cycle process was also suppressed by the drought treatment 25 . Furthermore, PAO (pheophorbide a oxygenase) is regarded as a vital chl catabolic enzyme 26,27 , and contributed well in chl deprivation process as its  Table S5). Buchert 28 and Du 29 have demonstrated the key role of PAO as an important chl degradation enzyme investigated during senescence of broccoli and banana, respectively. Meantime, the photosynthetic activity, stomatal conductance, and CO 2 assimilation rate were significantly decreased in drought treated grapevine leaves as compared to control. Similar findings have also been reported in grapevine under Cu and drought stresses 30,31 . Moreover, the photosynthesis-related genes, involved in PSII, PSI, cytochrome b6-f complex, photosynthetic electron transport, F-type ATPase and photosynthesis-antenna proteins were significantly down-regulated by drought stress, yet the extent of light-harvesting proteins (CP47, CP43), which binds the chla molecules was found down-regulated by drought stress (Supplementary  Table S6). Furthermore, PsaB is anticipated as the heart of PSI that binds P700 special chlorophyll pair 32 was also down-regulated by drought stress in our findings. Finally, drought stress gradually decreased the activities of PSII electron transport and light-harvesting complex (photosynthesis-antenna proteins). Stomatal closure is known to reduce the CO 2 absorption which limits the photosynthetic activity in plants under drought stress environment [33][34][35] . In the present study, the results of physiology and transcriptome analysis revealed that the decline in photosynthetic process was primarily connected with the chl degradation, which indicates that drought stress has direct influence on the photosynthesis process in plants.
ROS is a universal response of the plants to counter environmental stress to prevent oxidative damage. Several studies have already been performed to investigate the importance of MDA under oxidative stress in different crops, such as wheat (Triticum aestivum) and oilseed rape (Brassica napus), proposed that MDA was accumulated by drought stress [36][37][38] . On the contrary, plants have the ability to accrue different antioxidative enzymes to confer drought severity, while similar investigations in olive 39 and wheat 37 support our findings of enhanced activity of ROS enzymes against elevated levels of MDA. The transcriptomic investigation revealed that one NADPH oxidase and five amine oxidases were remarkably up-regulated, while both have a crucial role in the ROS synthesis and accumulation under various stress environments 40 . SODs are regarded as the first line of defense against ROS which has two isozymes Fe-SOD and Cu/Zn-SOD in plant chloroplast 40 . It is worth mentioning that Fe-SODs was up-regulated, but Cu/Zn-SODs was down-regulated, which is in agreement with our previous findings in grapevine under Cu stress conditions 30 . Other enzymes, including CAT, POD, GSH-AsA cycle, PPO, GST, AO, MDHAR, DHAR and GR also possess drought-responsive antioxidative defense system in grapevine 41 . Perhaps, non-enzymatic antioxidants, such as glutathione and proline also assisted in enhancing the ROS scavenging system in drought treated grapevine, which is consistent with the ROS scavenging system investigated in the V. vinifera and S. lycopersicum under drought stress 42,43 . The observed higher antioxidant capacity, increased activity of ROS enzymes and enhanced expression of genes-related to ROS system seems to be a mechanism operative in plant tolerance to drought stress.
Drought stress causes dehydration in plant cells. Plant hormones, such as abscisic acid, auxin, gibberellin, ethylene, jasmonic acid and brassinosteroid accumulate under dehydration conditions and play a major role in stress tolerance in plants 44 . In Arabidopsis, ABA activates the subclass III protein kinases of SnRK2 family, which further facilitate the regulation of stomatal conductance to optimize plant water status through guard cells 45,46 , favor our findings of increased activity of SRK2I protein kinase in drought-treated grapevine. The activation of PP2C genes in grapevine responding to drought stress proposed that PP2C has its primary role in stress tolerance, especially in regulating stomatal responses to cope transpiration losses 47 . The AUX gene family includes early response AUX genes, Aux/IAA, GH3 and SAUR and the regulators of AUX genes, ARF, while their activities were down-regulated in our findings. ARFs are also the important abiotic stress-responsive genes and have their crucial role in physiological processes of fruit plants. Wang et al. 44 investigated the AUX gene family in sorghum (Sorghum bicolor) and proline related protein 6 6 12 Table 5. List of differentially-expressed genes related to heat-shock proteins (HSPs) and pathogens resistance (PRs) proteins in grapevine perceived during drought stress.
specified that most of these genes were induced by the exogenous application of IAA under drought stress conditions. Moreover, GA activity and the accumulation of DELLA proteins were up-regulated in drought treatment, while similar findings in Arabidopsis have suggested that DELLA proteins restrain the plant growth and promote survival under drought stress 48 . JA biosynthesis and signaling together with ABA and other hormones have been extensively studied in many crops. In current investigations, JA amino acid conjugate (JAR1) was up-regulated, while JAR1 are enduringly present in the plant leaves and together with ABA induce the stomatal closure under osmotic stress, have been extensively studied in Arabidopsis 49 . Interestingly, jasmonate-zim domain proteins (JAZ) were significantly down-regulated, which was observed up-regulated in another study in rice 50 , may be because of prolonged duration of drought stress. Moreover, the activity of AOS and LOX were increased, which is similar to the findings of Leng et al. 36 in V. vinifera. Ethylene is an important stress hormone because its synthesis is induced under different oxidative environments. Under drought stress, the synthesis of ethylene precursor 1-aminocyclopropane-1 -carboxylate oxidase was up-regulated in grapevine, which can stimulate plant development and functioning by prompting the diffusion possibility of ABA to its active site 51,52 . Furthermore, the expressions of the ethylene-related regulatory genes (ETR1 and CTR1) were increased in our findings, suggesting their functional role in ethylene biosynthesis as described by Schachtman and Goodger 53 . BRs are the only plant steroids, which elicit the expression of many genes, especially during stress environments. Brassinosteroid Insensitive 1 (BRI1) was up-regulated in our findings, which is known to play the key role in plant growth, morphogenesis and response to drought stress. Feng, et al. 54 created RNAi mutants for bdBRI1 in Brachypodium distachyon and suggested that this gene produces a dwarf phenotype with enhanced tolerance towards drought stress. BR signal transduction, from cell surface perception to activation of specific nuclear genes will be interesting to investigate in the future. Plants cope with environmental stress by the accumulation of certain compatible osmolytes, such as proline, which is known to induce drought-tolerance in plants 55 and up-regulation of all the genes related to proline metabolism is the clear evidence of induced grapevine tolerance in our study. Proline biosynthesis commenced with the phosphorylation of glutamate, which then converted into gulatamic-ϒ-semialdehyde by Pyroline-5-carboxylate synthetase (up-regulated). Similarly, arginine is converted into orthinine by arginase (up-regulated) and then into GSA by the ornithine-δ-aminotransferase (not-detected). GSA is then converted into pyrroline 5-carboxylate (P5C) by spontaneous cyclization. Finally, proline is synthesized from the P5C by P5C reductase (P5CR) enzyme 55,56 . In proline degradation pathway, proline is re-converted into P5C by Proline dehydrogenase (PDH; up-regulated) and then into glutamate by Pyrroline-5-carboxylate dehydrogenase (P5CDH; up-regulated). Thus PDH and P5CDH are believed to be most important enzymes in proline degradation to glutamate 57,58 . Hence, proline metabolism may regulate the gene expression during drought stress.
In higher plants, accumulation of various secondary metabolites, such as amino acids, carbohydrates and lipids occur when a plant is subjected to environmental stress 59 . Shikimate pathway not only acts as a bridge between central and secondary metabolism but also serves as a precursor for other secondary metabolites 60 . Additionally, Tyr is a precursor of IAA and initiate the synthesis of indole alkaloids and isoquinoline alkaloids, which prevent plants from oxidative stress 61 . Phe is considered as the precursor of secondary metabolites family, and PAL participates in phenylpropanoid biosynthesis; a key step towards biosynthesis of stilbenes, flavonoids, lignins and various other compounds 62 . STS (stilbene synthase) catalyzes the initial step of flavonoid biosynthesis pathway, which has the protective function during drought stress 63 . Overall, 4 PAL and 6 STS were significantly up-regulated in our findings, proposing the innate link with drought stress. The respective, up and down-regulation of 1-deoxy-D-xylulose 5-phosphate reductoisomerase and 1-deoxy-D-xylulose-5-phosphate synthase can act as a rate limiting enzymes in MEP pathway, also found in Cu-stressed grapevine leaves 30 . Dimethylallyl diphosphate and isopentenyl diphosphate are the universal 5 carbon precursors found in terpenoid synthesis. It has been reported that one isopentenyl-diphosphate isomerase II can catalyze isopentenyl diphosphate to form dimethylallyl diphosphate and one terpene synthase 64,65 , while both were up-regulated in our findings. The down-regulation of most of the genes related to anthocyanin, lignin and terpenoid biosynthesis have elucidated the negative role of drought stress on an accumulation of secondary metabolites in grapevine leaves.
HSPs are ubiquitous stress-related proteins that act as molecular chaperone, HSP members participate in the protein synthesis, folding, aggregation and transportation from the cytoplasm to different intracellular compartments 66,67 . In the present study, majority of HSP genes were up-regulated (35), however genes related to sHSP along with heat stress transcription factors and other HSP were found to be predominant in their expression, in contrast, it was noted that drought has also significantly suppressed few genes (14) related to HSP. Heat-shock protein 70 (Hsp70) proteins are one of the large families of highly conserved molecular chaperones and are extensively found in almost all organisms. HSP70 is one of the highly conserved protein, known to induce during environmental stress, for example, in Erianthus arundinaceus, HSP70 was expressed 7-fold higher under drought-stressed condition 68 . sHSPs are plant tolerant proteins generally induced upon abiotic stress 69 . In the same way, Vasquez reported that expression of HSP70 and sHSP was higher during stress conditions in pine tree, indicating their major role in mitigating the drought stress 70 . Therefore, transcriptional analysis of the present study revealed that HSP also play a fundamental role in defense mechanism of grapevine during drought stress. In addition, pathogenesis-related (PR) proteins are derived from plant allergens and act as defense-responsive proteins by increasing their expression under pathogen attack and variable stress environments. Depending on the functions and properties, PR-proteins are classified into 17 families, such as beta-1,3-glucanases, chitinases, thaumatin-like proteins, peroxidases, small proteins (defensins and thionins) and lipid transfer proteins (LTPs) 11,71 . Most of the PR-proteins were down-regulated in our study, suggesting that drought stress posed a negative effect on PR-proteins defense response. Contrarily, most of the genes related to dirigent-proteins (DIR), play a role in lignin formation and proline-related proteins were up-regulated, suggesting their possible defensive role in grapevine in response to drought stress.

Conclusion
Our results have provided substantial shreds of evidence to demonstrate that grapevine adaptation to drought stress is a multi-step component system consisting of several genes that regulate various pathways. Out of 12,451 DEGs, 8021 transcripts were up-regulated and 4,430 transcripts were down-regulated. Nearly 2 fold up-regulations of DEGs have clearly indicated their defense-related role in grapevine responding to drought stress. The significant increase in the activity of ROS enzymes and hormones level revealed the defensive function of these enzymes and hormones during drought stress in grapevine leaves. Transcription analysis unveils that drought has affected the overall physiology of the grapevine; however, the regulatory mechanism of certain key genes like chlorophyll synthesis, ROS system, HSPs and other defense-related pathways assist grapevine in mitigating drought severity.

Materials and Methods
Plant material and drought treatments. Two-year old grapevine (V. vinifera cv. 'Summer Black') pot grown plants were selected as experimental material which were grown in standard greenhouse condition (25 ± 5 °C) under 16-h light/8-h dark photoperiod and 65% relative humidity (RH) at the Nanjing Agricultural University-Nanjing, China. Grapevine plants were subjected to drought with an interval of 20 days against control, each with three biological replicates. The fourth unfolded leaf from the shoot apex was collected from each replicates of both control and drought treatment with the interval of 0 and 20th day, respectively, and the three samples were mixed to make one composite sample. After harvesting, the samples were immediately put in liquid nitrogen and then stored at −80 °C until analysis.
Determination of important biochemistry and physiology-related traits. The chlorophyll a and b contents was determined using spectrophotometer at 663 and 645 nm, respectively as briefly explained by Leng, et al. 26 . Photosynthesis activity, stomatal conductance and CO 2 assimilation rate were carried out on mature leaf between 4 th to 7 th nodes from the shoot base for both control and drought treatment; between 9:00-11:00 AM measured using LI-COR (LI-6400XT, Germany) meter as described by Tombesi et al. (2015). Malondialdehyde (MDA) contents were quantified by using thiobartiburic acid. The activities of antioxidant enzymes (SOD, POD and CAT) were measured using the method briefly described by Haider, et al. 72 . The activities of indole-acetic acid (IAA), abscisic acid (ABA), jasmonic acid (JA), gibberellic acid (GA) and brassinosteroid (BR) were measured following the method of Tombesi, et al. 31 . Three technical repeats were generated for all the quantifications. Data was subjected to one-way analysis of variance (ANOVA) at p < 0.05, using MINITAB (ver. 16) and represented as mean ± standard deviation (SD).
RNA extraction, cDNA library construction and Illumina deep sequencing. Total RNA from leaf samples of both control and drought-stressed were extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA) (1% agarose gel buffered by Tris-acetate-EDTA was run to indicate the integrity of the RNA.) and subsequently used for mRNA purification and library construction with the Ultra ™ RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's instructions. The samples were sequenced on an Illumina HiseqTM2500 for 48 h.

Analysis of gene expression level, gene ontology (GO) and Kyoto encyclopedia of genes and genomics (KEGG).
After adaptor trimming and quality trimming, the clean reads were mapped to the V. vinifera transcriptome using Bowtie 1.1.2. Then, Sam tools and BamIndexStats.jar were used to calculate the gene expression level, and reads per kilobase per million (RPKM) value was computed from SAM files 73 . Gene expression differences between log 2 and early stationary phase were obtained by MARS (MA-plot-based method with Random Sampling model), a package from DEGseq. 3.3 (Leng et al., 2015). We simply defined genes with at least 2-fold change between two samples and FDR (false discovery rate) less than 0.001 as differential expressed genes. Transcripts with |log2FC| < 1 were assumed to have no change in their expression levels. The gene ontology (GO) enrichment (p-value < 0.05) was investigated by subjecting all DEGs to GO database (http://www.geneontology. org/) in order to further classify genes or their products into terms (molecular function, biological process and cellular component) helpful in understanding genes biological functions. Kyoto encyclopedia of genes and genomics (KEGG; the major public pathway-related database) was used to perform pathway enrichment analysis of DEGs 74 .
Illumina RNA-seq results validation by qRT-PCR. In order to validate the Illumina RNA-seq results, drought-stressed grapevine leaf samples of each collection were applied to qRT-PCR analysis. Total RNA of the collected samples was extracted following the above mentioned method, and then was reverse-transcribed using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China), following the manufacturers' protocol. Gene specific qRT-PCR primers were designed using Primer3 software (http://primer3.ut.ee/), for 20 selected genes with the sequence data in 3′UTR (Table S12). qRT-PCR was carried out using an ABI PRISM 7500 real-time PCR system (Applied Biosystems, USA). Each reaction contains 10 µl 2 × SYBR Green Master Mix Reagent (Applied Biosystems, USA), 2.0 µl cDNA sample, and 400 nM of gene-specific primer in a final volume of 20 µl. PCR conditions were 2 min at 95 °C, followed by 40 cycles of heating at 95 °C for 10 s and annealing at 60 °C for 40 s. A template-free control for each primer pair was set for each cycle. All PCR reactions were normalized using the Ct value corresponding to the Grapevine UBI gene. Three biological replicates were generated and three measurements were performed on each replicate.