Effect of mid- and long-term drought stress on physiological parameters of walnut plants
Overall, drought stress had significant effects (P< 0.01) on RWC, SPAD index and MSI at both 9th and 18th days of drought stress exposure. RWC in control plants was 81.2% and 80.23% at 9th and 18th d respectively, but under drought stress this index reduced by 61.1% and 48.7% at 9th and 18th d. Also, SPAD index reduced in water deficit-exposed plants from 38 to 30.3 at day 9th and 37.2 to 23.4 at day 18th of drought stress exposure. Finally, MSI was reduced by drought stress conditions (35.4% and 31.1% at 9 and 18 d) than control (57.3% and 53.6% at 9th and 18th d) (Table 1).
Table 1
change in RWC, SPAD and MSI in control and drought stress of walnut in 9th and 18th d after the start of experiment
|
RWC (%)
|
SPAD
|
MSI (%)
|
|
Control
|
Drought stress
|
Control
|
Drought stress
|
Control
|
Drought stress
|
9d
|
81.2
|
61.1
|
38
|
30.3
|
57.3
|
35.4
|
18d
|
80.23
|
48.7
|
37.2
|
23.4
|
53.6
|
31.1
|
Significant
|
**
|
**
|
**
|
At 18th d of drought stress exposure parameters related to fluorescence of chlorophyll were measured to evaluate the healthiness of photosynthesis apparatus of plants. Some of indexes such as Fv/Fm, ABS/RC, PI ABS, DIo/RC and Fv/Fo affected under drought conditions (Fig. 1). The parameter related to electron transport and functionality of the photosynthesis apparatus including Fv/Fm, PI ABS and Fv/Fo decreased in the plants under drought stress, simultaneous ABS/RC and DIo/RC parameter that associated to dissipation of energy increased under drought stress.
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).
Functional Annotation
All the assembled transcripts were analyzed and compared against public databases (NCBI nr database, UniProtKB database, Pfam database, Rfam database, miRbase database, PlnTFDB database, SignalP software, Rnammer software and tmHMM software) by different tools including BLAST, HMMER, TransDecoder, RNAMMER v1.2, TMHMM v2.0, SignalP v4.1, and KASS. The results revealed that 61% (111,451 transcripts) of 183,191 transcripts were annotated using one or more databases. Briefly, a total of 104,926 (57.27%) ORFs, 79,185 (43.22%) transcripts with at least one predicted domain, 3,931 (2.14%) transcription factor, 6,591 (3.59%) signal peptides, and 92,704 (50.6%) transmembrane domains were identified. This analysis identified 109,413 (59.72%) and 84,384 (46.06%) hits from the NCBI nr and the UniProtKB databases, respectively. In order to better comprehend the biological pathways in the transcripts, the annotated sequences were analyzed by KEGG database. Results of pathway analysis showed that all the transcripts were classified to five specific pathways including Organism Systems, Cellular Processes, Environmental Information Processing, Genetic Information Processing and Metabolism (Fig. 2). The largest number of the transcripts were belonged to Metabolic pathway (8,050, 32.58%) followed by Organism Systems (5,630, 22.79%), Environmental Information Processing (3,951, 15.99%), Genetic Information Processing (3,786, 15.32%) and Cellular Processes (3,284, 13.29%). Moreover, Signal transduction (3,857, 15.68%) and Carbohydrate metabolism (2,028, 8.21%) were the top two representative pathways (Fig. 2).
Gene Ontology (GO) classification of all the annotated transcripts were performed based on the three functional groups; biological process (BP), molecular function (MF) and cellular component (CC). Out of 183,191 transcripts, 81,494 transcripts were contained at least one GO term (cellular compounds (69,717 transcripts), molecular functions (70,689 transcripts) and biological processes (64,992 transcripts)). All the annotated transcripts were classified 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.
Identified 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 significant 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 significantly enriched in the down-regulated genes including ‘leaf morphogenesis’, ‘polysaccharide biosynthesis’, ‘cell wall organization’ and ‘acetyl-coA biosynthesis’ pathways were significant (Fig. 4b). On the other hand, KEGG pathway enrichment analysis revealed four and three significant enriched pathways in up-regulated and down-regulated genes, respectively (Table 2).
Table 2
Important genes for each pathway with fold change represented in walnut under drought stress in 9th d.
|
Pathways
|
Gene name (Fold change)
|
Up-regulated
|
Oxidative stress
|
BCB (14.2), SIP2 (10.49), CYC6 (9.34), APXS (9.16), PXD (5.3) and OSA1 (2.19)
|
Carbohydrate metabolism
|
MDI (9.9), BGAL1 (8.3), R1 (4.03), Lichenase (2.4) and HGN1 (2.4)
|
Cell wall macromolecule catabolism
|
HSP7P (10.77), KTI2 (3.96) and CALS3 (11.41)
|
Secondary metabolite biosynthesis
|
CYP71 (6.46), MO3 (4.54), CYP89 (3.66) and LKR/SDH (2.25)
|
Down-regulated
|
Leaf morphogenesis
|
TCP4 (1.64), GAE6 (1.95), LNG1 (2.09), TCP2 (3.31) and LUH (7.5)
|
Polysaccharide biosynthesis
|
GATL9 (9.32), AXS2 (1.5), GAE6 (1.95) and GATL9 (2.65)
|
Acetyl-coA biosynthesis
|
ACS (4.4)
|
Identified DEGs between control (C18) and drought stress (S18) at 18th d
A total of 1035 DEGs were identified 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 down-regulated 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 significant 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)
|
Up-regulated
|
Response to drought stress
|
MPK3 (10.8), MAPKKK17 (6.69), Map3k1 (7.71), ABR1 (4.95), RMA1H1 (4.6), EDL3 (4.71), XERICO (3.5),
|
Abscisic acid
|
AHG1 (9.58), ARIA (9.26), CCA1 (8.81), NCED6 (6.83), SAG21 (5.14), ABA4 (4.6), XERICO (3.5), AFP2 (4.35), NCED1 (3.15), MARD1 (2.79),
|
Secondary metabolite biosynthesis
|
CYP71B37 (5.91), CYP82A3 (5.42), MO3 (5.8), CYP71B34 (4.62)
|
Carbohydrate metabolism
|
rhiE (12.1), Chi3l1 (9.26), rfs5 (8.12), GNS1 (5.48), HGN1 (5.13), GN1 (4.81), NANA (3.7)
|
Response to oxidative stress
|
MAPK3 (10.83), GATL10 (5.5), SAG21 (5.14), CSE (3.97), HSPRO2 (3.36)
|
Response to salt stress
|
CCA1 (8.81), GLP9 (8.35), SRO5 (7.8), GLYI (4.44), MKK9 (4.11), SRM1 (3.9)
|
Regulation of stomatal opening
|
ALIS1 (10.9), IQM1 (8.87), MYB15 (4.92), RZPF34 (4.89), CKX6 (4.6)
|
Down-regulated
|
ATP biosynthesis
|
WAK2 (2.79), AHA10 (2.7), BAM3 (2.66), ABCC3 (2.73), AB3C (3.04)
|
vacuole organization
|
DTX41 (2.33), ABCC3 (2.73), YPQ1 (5.86), MON1 (8.16)
|
starch catabolism
|
WAXY (2.77), LECRK1 (2.65), BAM3 (2.65)
|
Functional KEGG pathways for DEGs in 18 day after the start of experiment
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 significant with different p value. According to the Table 4. the important pathways that significantly 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.
Table 4
Functional analysis up-regulated DEGs based on the KEGG pathways in walnut 18th d after the beginning of experiment.
Enriched pathway term
|
ID
|
Number of transcripts
|
P-Value
|
Transcription Factors
|
ko03000
|
20
|
0.02
|
Plant hormone signal transduction
|
ko04075
|
18
|
0.00025
|
Starch and sucrose metabolism
|
ko00500
|
14
|
0.00056
|
Galactose metabolism
|
ko00052
|
12
|
4.18×\({10}^{-7}\)
|
Mitogen-activated protein kinase (MAPK) cascade
|
ko04016
|
11
|
0.015
|
Glycolysis
|
ko00010
|
9
|
0.032
|
Tryptophan metabolism
|
ko00380
|
6
|
0.00012
|
Isoquinoline alkaloids biosynthesis
|
ko00950
|
5
|
0.00027
|
Lysine degradation
|
ko00310
|
5
|
0.0081
|
Tyrosine metabolism
|
ko00350
|
5
|
0.009
|
Cytochrome P450
|
ko00982
|
5
|
0.014
|
Ascorbate and aldarate metabolism
|
ko00053
|
5
|
0.015
|
Alpha-Linolenic acid metabolism
|
ko00592
|
5
|
0.016
|
Amino acid metabolism
|
ko00250
|
5
|
0.018
|
Glycine, serine and threonine metabolism
|
ko00260
|
5
|
0.02
|
Beta-Alanine metabolism
|
ko00410
|
5
|
0.025
|
Metabolism; Amino acid metabolism
|
ko00280
|
5
|
0.025
|
Histidine metabolism
|
ko00340
|
3
|
0.036
|
Betalain biosynthesis
|
ko00965
|
2
|
0.00018
|
Table 5
KEGG pathway analysis of down-regulated DEGs between the control and drought stress groups in 18th d after the start of experiment.
Enriched pathway term
|
ID
|
Number of transcripts
|
P-Value
|
Transporters
|
ko02000
|
14
|
0.0013
|
Flavonoid biosynthesis
|
ko00941
|
13
|
3×\({10}^{-17}\)
|
Ubiquinone and other terpenoid-quinone biosynthesis
|
ko00130
|
6
|
3.42×\({10}^{-6}\)
|
Phenylpropanoid biosynthesis
|
ko00940
|
4
|
0.0031
|
Glutathione metabolism
|
ko00480
|
4
|
0.0016
|
Circadian rhythm
|
ko04712
|
3
|
0.001
|
Identified 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 identified, 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 significantly enriched 10 biological terms including ‘regulation of photosynthesis’, ‘stomatal complex morphogenesis’ and ‘positive regulation of response to oxidative stress’. Fig. 7a and Fig. 7b depicted the number of common and unique up-regulated and down-regulated DEGs in the three performed comparisons (C9-S9, 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 identified. A total of 28 TFs belonged to 10 families were identified including MYB, NAC, WRKY and ERF families. Most of the identified 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).
Table 6
Transcription factors identified in walnut transcriptome under drought stress at 18th d after the beginning of the experiment.
Transcription factor
|
Up-regulated (Fold change)
|
Transcription factor
|
Up-regulated
(Fold change)
|
MYB78
|
10.77
|
NAC029
|
6.01
|
MYB3R1
|
9.38
|
NAC073
|
5.07
|
MYB2
|
6.03
|
NAC072
|
4.83
|
MYB73
|
5.31
|
NAC100
|
3.8
|
MYB108
|
5.12
|
NAC002
|
3.5
|
MYB15
|
4.9
|
ERF1
|
7.99
|
MYB69
|
4.8
|
ERF1B
|
5.3
|
MYB3
|
3.1
|
ERF1A
|
3.6
|
WRKY40
|
10.11
|
ERF4
|
3.56
|
WRKY75
|
6.57
|
SIB1
|
8.3
|
WRKY14
|
5.8
|
RAP2
|
2.9
|
WRKY6
|
4.08
|
EIN3
|
2.26
|
ARF5
|
8.8
|
EDL3
|
3.27
|
ARF6
|
3.3
|
ASIL2
|
9.19
|
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, confirming the reliability of the RNA-Seq data.