Genomic Responses Associated with Drought Stress in Walnut as Revealed Transcriptome Sequencing

Walnut production is challenged by abiotic stresses. We investigated the leaf transcriptome responses of walnut under control and drought stress in 9 and 18 days. We identied 921, 1035 differentially expressed genes (DEGs) between control and drought stress groups in 9 and 18-day, respectively. In control and drought stress conditions DEGs were signicantly enriched into the abscisic acid biosynthesis, regulation of stomata closure, leaf morphogenesis, carbohydrate metabolism, oxidative stress, cell wall macromolecule catabolism, and secondary metabolite biosynthesis pathways. We conrmed our RNA-Seq data using quantitative real-time PCR (qPCR) of six candidate genes. Our results indicated that more complicated transcript regulation of drought responses following prolong exposure to drought stress. In general, walnut activated more tolerance mechanisms 18 days after drought stress. Findings of this research would be useful for future studies on breeding for drought tolerance of Persian walnut and related species.


Introduction
Persian or English walnut (Juglans regia L.) is native to the central Asia, but is now distributed in worldwide in temperate areas for nut production 1 . The highest walnut production in the world belongs to the United States, China and Iran, respectively. From a nutritional point of view, walnuts provide appreciable amounts of proteins, carbohydrates and are rich in unsaturated fatty acids including, linoleic, linolenic acids and oleic 1 . In a large part of the areas where walnuts are produced, researchers have discovered that drought stress severely reduced yield [2][3][4] . In the current changing global climate scenario, water availability is one of the most limiting factors that determine plant distribution and productivity. Responses of plants to drought stress may occur from cell membrane to photosynthetic organs, to the whole plants 5 . Drought stress in plants is a complex phenomenon that triggers diverse modi cations from gene to morphology 6 . To survive, plant provides physiological (e.g. stomatal movements, ABA signaling, membrane transport, osmolyte regulation), biochemical (enzyme and non-enzyme activity), morphological (changes in leaf and root growth, wilting of leaves) and molecular (changing the expression of a number of genes, transcription factors and micro-RNAs) responses 7,8 . Hence, comprehensive surveys of mechanisms by which walnuts are affected by drought stress are needed, especially in arid and semi-arid regions. Thanks to next-generation sequencing (NGS)-based technologies, gene expression pro ling using RNA sequencing method has become a highly sensitive and accurate method, which can be applied to study the plant responses to biotic and abiotic stresses 9 . Recently, many studies have been conducted to understand the molecular mechanisms of drought tolerance using comparative transcriptome by RNA-Seq in many plant species such as tomato 10 , sorghum 9 , cotton 11 , camellia 12 , pear 13 , peach 14 , eucalyptus 15 , California oak 6 , and Jhar Ber 16 . In spite of drought stress importance, most of the studies considered morphological and physiological indices and no studies have been investigated the molecular mechanisms (including gene expression pro ling) in response of walnut to drought stress. A study reported that seed germination of walnut in drought stress conditions signi cantly reduced and soluble sugars and proline accumulate in shoots and roots in tolerant genotypes 17 . It has been also reported that petioles loss 87% of their maximum hydraulic conductivity in drought-exposed walnuts, causing leaf abscission 18 . In spite of these studies, the molecular mechanisms underlying drought stress on walnut remains largely unknown and needs to be further investigated. Therefore, in this study, RNA-Seq technology was used to understand the genome-wide transcriptome changes of walnut following mid and long-term exposure to drought stress condition and to characterize differential physiological response of walnut under normal and drought conditions. Results of this study will expand the knowledge of lncRNAs participating drought response in walnut and may be applied to improve productivity of walnuts under drought stress.

Results
Effect of mid-and long-term drought stress on physiological parameters of walnut plants

Rna Sequencing And Assembly
The plants were kept under control and drought stress conditions for 18th d. Leaf sampling was performed at two time points, 9th (C9 and S9) and 18th (C18 and S18) d after the start of the experiment, and three biological replications were conducted for each group. In total, 240 million paired-end reads were generated of 11 libraries with an average of 21.8 million reads per library. After removing the adapters, low-quality reads and correction of errors, 231 million clean reads were remained for further analysis (Table S1). Our de novo assembly pipeline, based on various quality metrics of assembly, showed that the transcriptome assembly produced by EvidentialGene can be prioritized over other assemblies in walnut leaves. The reference transcriptome assembly generated by EvidentialGene included 183,191 transcripts with an N50 length of 1,831 bp, among which 64,702 transcripts were longer than 1 kb (Table S2) ). All the annotated transcripts were classi ed to 56 functional categories using WEGO plot (Fig. 3). In the biological process category, most of the transcripts were related to 'cell', 'cell part', 'organelle' 'membrane part' and 'protein-containing complex'. Also, 'catalytic activity', 'binding', 'transcription regulator activity' and transporter activity' had highest percentage in the molecular function category. Finally, in the biological process category 'cellular process', 'metabolic process', 'biological regulation', 'response to stimulus' and 'regulation of biological process' were the most representative categories. A number of organelle plants including cell, cell part and membrane part are affected at the beginning of drought stress.
Identi ed DEGs between control (C9) and drought stress (S9) at 9th d The analysis showed that 921 genes were DEGs, of which 580 genes were up-regulated and 341 genes were down-regulated in S9 than C9 group. Results of the functional enrichment analysis of up-regulated genes revealed 27 signi cant GO terms in biological process category including 'secondary metabolite biosynthesis', 'response to oxidative stress', 'carbohydrate metabolism', 'polysaccharide digestion' and 'cell wall macromolecule catabolism' pathways ( Fig. 4a). Also, 11 biological process terms were signi cantly enriched in the down-regulated genes including 'leaf morphogenesis', 'polysaccharide biosynthesis', 'cell wall organization' and 'acetyl-coA biosynthesis' pathways were signi cant (Fig. 4b). On the other hand, KEGG pathway enrichment analysis revealed four and three signi cant enriched pathways in up-regulated and down-regulated genes, respectively ( Table 2). Identi ed DEGs between control (C18) and drought stress (S18) at 18th d A total of 1035 DEGs were identi ed between the control and drought stress group in 18 day after the start of experiment (P< 0.05). Of these, 850 genes were up-regulated, while 180 genes were downregulated in S18 than C18 group. GO enrichment analysis of up-regulated DEGs showed that the 'response to drought stress', 'abscisic acid', 'secondary metabolite biosynthesis', 'response to salicylic acid', 'response to oxidative stress', 'response to salt stress', 'response to osmotic stress', 'polysaccharide digestion', 'regulation of stomatal opening', 'malate transport' and 'carbohydrate metabolism' pathways were enriched (Fig. 5a). Likewise, in the down-regulated DEGs 'photosystem stabilization', 'ATP biosynthesis', 'negative regulation of response to water deprivation', 'vacuole organization' and 'starch catabolism' pathways were signi cant terms (Fig. 5b). KEGG pathway analysis predicted the involvement of up-regulated and down-regulated DEGs in seven and three pathways, respectively (Table 3). Table 3 Important genes for each pathway with fold change represented in walnut under drought stress in 18th d.

Pathways
Gene name (Fold change)

Upregulated
Response to drought stress Since the more pathways were activated 18 days than 9 day, up-regulated pathway enrichment conducted using the KEGG database. There were 20 KEGG pathway recognized as enriched signi cant with different p value. According to the Table 4. the important pathways that signi cantly in drought stress conditions including 'Transcription Factors', 'Plant hormone signal transduction', 'Starch and sucrose metabolism', 'Galactose metabolism' and 'Mitogen-activated protein kinase (MAPK) cascade'. Also, the pathway enrichment for down-regulated genes were showed in the Table 5.  Identi ed Degs Between Drought Stress (S9) And (S18) Groups A comparison was conducted to identify DEGs in drought stress groups at 9th and 18th d after the start of experiment. Overall, 841 DEGs between S9 and S18 groups were identi ed, of which 769 genes were up-regulated and 76 genes were down-regulated in S18 than S9 groups. On the basis of GO analysis 33 biological terms including 'carbohydrate metabolism', 'activation of MAPKK activity', 'response to osmotic stress', 'abscisic acid metabolism', 'proline catabolism' and 'secondary metabolite biosynthesis' were enriched in up-regulated DEGs (Fig. 6a). As displayed in Fig. 6b, down-regulated DEGs were signi cantly enriched 10 biological terms including 'regulation of photosynthesis', 'stomatal complex morphogenesis' and 'positive regulation of response to oxidative stress'. C18-S18 and S9-S18), respectively.

Transcription Factors (Tfs) In Degs
Because of the importance of TFs in gene expression regulations, especially in response to drought stress, TFs in the DEGs were identi ed. A total of 28 TFs belonged to 10 families were identi ed including MYB, NAC, WRKY and ERF families. Most of the identi ed TFs were belonged to MYB group, as MYB78 and MYB3R1 were up-regulated and showed the highest fold change (10.77 and 9.38 fold change, respectively) ( Table 6). Validation Of Degs Using Rt-pcr Validation of DEGs using RT-PCR To validate the RNA-seq results, six genes including four up-regulated (RCD1, ALDH, PGI1 and MPK3) and two down-regulated (ERL2 and LECRK1) were randomly picked from DEGs. The ratio of the Log2 fold change from the RNA-seq data was compared to the Log2 fold change obtained with RT-PCR. As shown in Fig. 8, all the six genes displayed similar expression patterns in both RNA-Seq and RT-PCR results, con rming the reliability of the RNA-Seq data.

Discussion
Drought stress is a multidimensional stress causes diverse modi cations including metabolic adjustment, transcriptomic changes, osmotic adjustment, leaf turgor pressure alterations, proteomic changes and retardation or deceleration of plant growth [19][20][21] . In the present study, RWC, SPAD index and MSI were down-regulated by 9 th and 18 th d exposure to drought stress condition. These ndings are in line with the previous reports that identi ed down-regulation of physiological parameters by drought stress in Persian walnut 2,17,22,23 . In addition, parameters related to the photosynthesis functionality such as Fv/Fm, PI ABS, DIo/RC and Fv/Fo related to photosynthesis, were negatively affected by drought stress condition at 18 th d (Fig. 1). These ndings con rms occurrence of photoinhibition in water de citexposed walnut plants. Imbalance between products of light-dependent reactions (ATP and NADPH) and consumption of these products in Calvin-Benson cycle is usually occurred as a result of long-term stomatal limitation for provision of enough CO 2 as the substrate for RUBISCO enzyme in Calvin-Benson cycle 24,25 . This problem emerges because there are many enzymatic reactions in Calvin-Benson cycle and they decelerate when they are working under low mesophyll CO 2 concentration 26 . Although at the beginning of drought stress exposure there would be no limitation on light-dependent reactions and they are still running and transfer the electrons to NADP + to produce NADPH; after a while the NADP + pool would be limited due to down-regulation of enzymatic reactions in Calvin-Benson cycle. This imposes an imbalance between two sets of the reactions directing the excited electrons toward oxygen and produces reactive oxygen species (ROS) 27 . Reduction in MSI as an indicator for oxidative stress can be related to the damages caused by ROS accumulation. It seems rational to postulate that these results are due to that drought stress induced much more adverse photoinhibitory impairment in drought stress (S18) than control (C18) conditions. Changes in the above indicators indicated that plants were continuously damaged by drought stress in 9 th and 18 th d.
KEGG pathway enrichment analysis showed that most of transcripts were related to metabolic pathway than the other pathways (Fig. 2). In this regard, carbohydrate metabolism and amino acid metabolism included most of transcripts (2028 and 1170 transcripts, respectively). Previous studies reported that carbohydrate metabolism was signi cantly enriched under drought stress 28-30 . Moreover, a large number of the transcripts (3,847) were annotated to be participated in signal transduction pathway. Many reports are available in the literature that described rapid responses of plants to drought stress (molecular changes, TFs activation, hormones transmission of hormones and stomatal closing) often mediated by signal transduction pathway 29,30 . According to the GO enrichment analysis, 'metabolic process, 'response to stimulus', and 'signaling' biological process terms had the most transcripts (Fig.   3). The ndings of the current study are consistent with other studies such as in apple 31 , sweet orange 32 and apricot 33 which reported that under drought stress 'metabolic process' and 'response to stimulus' are the highest transcripts in biological process category.
Response to oxidative stress was one of the enriched pathways in DEGs in S9 to C9 comparison. Drought condition leads to an imbalance between the production of ROS and the plant's antioxidant defenses, which causes oxidative stress in plants. For 37 . In the current study, secondary metabolite biosynthesis pathway activated in drought stress groups. The important gene belonged to secondary metabolite biosynthesis which up-regulated were CYP71, MO3, CYP89 and LKR/SDH. CYP71 and CYP89 genes belonged to cytochrome P450, previously reported to the highly expression by environmental stress 31,38,39 . Pathways related to leaf morphogenesis down-regulated signi cantly in drought stress at 9 day after start of experiment (Table. 2). In the current study TCP4 and TCP2 down-regulated, this indicated that walnut under drought conditions changes of leaf morphology. Previous research report that under drought stress leaf morphogenesis regulated by miR319 targets to the gene family of TCP transcription factors 40 . Overall in the present study, more drought tolerance pathways were activated at 18 days than at 9 days after the start of the experiment. The important pathways including up and down-regulated showed in Fig. 6. Biosynthesis and catabolism of ABA is essential for drought tolerance in plants, including stress-responsive gene expression, metabolic changes and stomatal closure 41 . ABA biosynthesis pathway enriched signi cantly in S18 than C18, list of gene expression in this pathway given in Table. 3. NCED6 encodes 9-cis-epoxycarotenoid dioxygenase, a key enzyme in the biosynthesis of abscisic acid 42 . ABA4 is essential for neoxanthin biosynthesis that intermediary step in abscisic acid biosynthesis 43 . Other genes play a key role in the synthesis and transport of ABA including ARIA 44 , AHG1 45 , AFP2 46 and XERICO 47 . Other pathway that enriched signi cant in S18 is regulation of stomatal opening, actually, the role of this pathway is to close the stomata opening in collaboration with the ABA biosynthesis pathway in order to prevent massive water loss from the plants 8 . One of the genes that encodes a calmodulin-binding protein involved plays key roles in stomatal movement is IQM1 48 . Response to oxidative stress, carbohydrate metabolism and secondary metabolite biosynthesis pathways activated in S9 and S18, and gene name with full change given in Table 2 and 3. These ndings suggest that walnut activate certain pathways in early drought, but at long droughts, different gene expressed (compare Table 2 and 3). For example, compared to its transcript levels in C18, HSPRO2 gene up-regulated in S18 group, but not expression in S9. In Arabidopsis transgenic overexpressing of HSPRO2 gene showed that high tolerance to oxidative stress than the wild type 49 . In the other study researcher reported that Arabidopsis transgenic with overexpression of GATL10 gene have increased galactionol and ra nose levels than the wild type plants 50 . Increasing the concentration of these two substances protects salicylic acid against the free oxygen radicals 51 . Interestingly, 'response to salt stress' pathway enriched signi cantly in the S18 than C18 (Fig.  5a). Several paper reported that there are similarities responses between salinity and drought in plants 51,52 . As shown in Fig. 5b ATP biosynthesis pathway down-regulated under drought stress in S18 group. This nding is in agreement with 53 report that the tricarboxylic acid (TCA) cycle and ATP biosynthesis are negatively affected under drought stress. Transcription factors such as WRKY, NAC, BZIP, MYB, ERF and ERF family members played vital roles in drought stress tolerance by regulation of hormone metabolism and hormone signalling transduction pathways 32,54 . Several groups of transcription factors have been identi ed under drought stress in S18 than C18 ( play a positive role in the growth and tolerance to drought stress in rice 66 . Also, reported that ERF1 transcription factor was directly involved in tolerance to drought stress in arabidopsis. This transcription factor is involved in the transfer of the ethylene, abscisic acid, and jasmonic acid 67 .

Conclusion
The present study was designed to determine the transcriptomic response of leaf against drought stress in walnut. We identi ed 921 DEGs between C9-S9, 1035 DEGs between C18-S18 and 841 DEGs between S9-S18 groups. Our results for GO enrichment showed that in day 9th of drought stress exposure, Oxidative stress, Carbohydrate metabolism, Cell wall macromolecule catabolism and Secondary metabolite biosynthesis pathways were up-regulated and Leaf morphogenesis, Polysaccharide biosynthesis and Acetyl-coA biosynthesis pathways were down-regulated in drought stress conditions.
Also, day 18th of drought stress exposure Response to drought stress, Abscisic acid, Secondary metabolite biosynthesis, Response to oxidative stress and Regulation of stomatal opening pathways were up-regulated and ATP biosynthesis, vacuole organization and starch catabolism pathways were down-regulated in drought stress groups. We also identi ed several groups of transcription factors associated with drought tolerance in walnut. A schematic model is represented for the transcriptional responses associated with the drought tolerance in walnut is shown in Fig. 9. This research extends our knowledge regarding the mechanisms that are involved in drought stress responses of walnut and provide useful genomic information for future studies.

Plant material, growth conditions and drought stress
Plant material and growth conditions in the current study described in the recent publication 67 . Brie y.
Two-year old clonally propagated walnut cv. 'Chandler' were used in present study. These plants were grown in two-liter pots containing peat and perlite (1/1, v/v) and kept in the greenhouse under control conditions (with 16 h light/8 h dark photoperiod). To apply drought stress to the chandler cultivar, eighteen pots were selected and divided in two groups including control and drought-exposed plants (i.e. nine pots per each group) for a total of three replicates, with three pots each. Control plants were watered once every two days with distilled water and drought stress was applied by water withholding. Finally, leaves samples were harvested from control and stress treatment at 9 and 18 days after starting treatments. The leaf samples were frozen in liquid N 2 and stored at -80 ˚C for subsequent uses.

Physiological Measurements
Relative water content (RWC) Fully developed young leaves were taken at 9 th and 18 th d after water holding from both control and water de cit-exposed plants between 9-10 am and RWC was calculated according the formula; (FW -DW)/(TW -DW) × 100 where FW is fresh weight, DW is dry weight (after drying disks in 70 ˚C for 48 h) and TW is turgid weight (adding approximately 3 ml of distilled water to disks for 4 h in room temperature).

Leaf chlorophyll content
Each treatment had ve replicates, and ve disks (10 mm 2 ) were used for each replicate. SPAD index (leaf chlorophyll content) was estimated nondestructively for ve leaves per treatment, using a portable SPAD-502 chlorophyll meter (Minolta Corp., Ramsey, NJ, USA).

Membrane stability index
Membrane stability index percentage was calculated in 20 leaf discs from young leaves using equation MSI= (1-C 1 /C 2 ) ×100 where C 1 and C 2 are electrical conductivity were measured before and after leaf discs incubated at 40°C according to the method described by 52 . (http://arthropods.eugenes.org/EvidentialGene/) and the Transfuse v0.5.0 (https://github.com/cboursnell/transfuse) into a single assembly. After assembly quality assessment the best assembler was chosen 67 . The clean reads were assembled with EvidentialGene version 2013.07.27;

Chlorophyll Fluorescence Measurements
(http://arthropods.eugenes.org/EvidentialGene/) software using the default setting. EvidentialGene pipeline enables us to reduce the complexity of the de novo transcriptome assemblies using remove the highly similar transcripts, sequence fragments and transcripts with low coding potential. In the following, quality of assembly was evaluated by N50 length, reads mapping back to transcriptome (RMBT) and Finally, edgeR package (version 3.26.4) of R software was used to perform count-based differential expression analysis. An false discovery rate (FDR) <0.001 and a normalized fold change ≥ 2 were set as the thresholds to identify DEGs. Overall, three pair-wise comparisons were considered including control (C9) and drought stress (S9) in 9-day, control (C18) and drought stress (S18) in 18-day and drought stress in 9 (S9) and 18 (S18) day after beginning the experiment.

Validation Of Rna-seq Results By Rt-pcr Analysis
A total of six DEGs (four up-regulated and two down-regulated) were randomly selected to validate the RNA-Seq results using RT-PCR. In brief, total RNA was extracted using TRIzol reagent (Invitrogen, USA).
Then, 2 µg of total RNA was used for rst-strand cDNA synthesis using the Thermo Scienti c kit (K1621, USA) according to the manufacturer's protocol. Primer sequences were designed using the Oligocalculator and primer3 programs (Table S3). The walnut Actin2 was used as a housekeeping gene to verify the expression results of sequencing. RT-PCR was conducted using a HIFI SYBR Green Master Mix (Mabnateb, Iran). The cycling instructions was 15 min at 95 ºC, followed by 40 cycle of 15 s at 95 ºC for denaturation. 60 s, at 60 ºC for annealing and 120 s at 60 ºC for extension. All reactions were done with triple biological replicates and two technical replicates. Relative expression of each gene was calculated using 2 ΔΔCT method.   enriched GO terms with log10 p-value < 0.05. The color and the size of bubbles show the p-value (legend in upper right-hand corner) and the frequency of the GO term in the underlying GOA database in REVIGO analysis, respectively (bubbles of more general terms are larger).  Comparison of the expression pattern achieved by RT-PCR and RNA-Seq.