Abstract
In this study, high-throughput sequencing (RNA-Seq) was utilized to
evaluate differential expression of transcripts and their related genes
involved in response to terminal drought in root tissues of bread wheat
landrace (L-82) and drought-sensitive genotype (Marvdasht). Subsets of
460 differentially expressed genes (DEGs) in drought-tolerant genotype
and 236 in drought-sensitive genotype were distinguished and
functionally annotated with 105 gene ontology (GO) terms and 77
metabolic pathways. Transcriptome profiling of drought-resistant
genotype “L-82” showed up-regulation of genes mostly involved in
Oxidation-reduction process, secondary metabolite biosynthesis, abiotic
stress response, transferase activity and heat shock proteins. On the
other hand, down-regulated genes mostly involved in signaling,
oxidation-reduction process, secondary metabolite biosynthesis,
auxin-responsive protein and lipid metabolism. We hypothesized that the
drought tolerance in “L-82” was a result of avoidance strategies.
Up-regulation of genes related to the deeper root system and adequate
hydraulic characteristics to allow water uptake under water scarcity
confirms our hypothesis. The transcriptomic sequences generated in this
study provide information about mechanisms of acclimation to drought in
the selected bread wheat landrace, “L-82”, and will help us to unravel
the mechanisms underlying the ability of crops to reproduce and keep
its productivity even under drought stress.
Introduction
Wheat (Triticum aestivum L.) is one of the main food crops, being
consumed by humans for more than 5,000 years [[30]1]. It provides
nearly 55 percent of the carbohydrates requirements of the world
population [[31]2]. Water deficit is considered to be among the most
severe environmental stresses [[32]3] and adversely affects more than
50 percent of the wheat production area in the world [[33]4].
Maintenance of productivity during adverse environmental conditions is
one of the priority areas for plant science studies [[34]5, [35]6].
Root system traits are considered to be important in maintaining plant
productivity under drought stress [[36]5, [37]7, [38]8]. Extensive and
deep root systems make the plant maintain higher water potential and,
therefore, have a longer period of evaporation under drought stress
condition [[39]9]. The optimal root system for drought condition is
controlled by a wide range of the genes that are activated to enhance
tolerance of crops to water deficiency [[40]7, [41]10]. Converting this
physical stress into a biochemical response will happen after the
perception and recognition of the external change. Different signaling
pathways are activated, each of them stimulating a set of
stress-responsive genes. These signal cascades induced by a given
stress finally lead to stress tolerance. On the other hand, the roots,
as the first recipient of drought stress, stimulate molecular responses
of the shoots through signaling cascades, allowing the plant to protect
itself against water stresses [[42]11]. It seems that plant landraces
can provide genetic resources that meet forthcoming challenges for
struggling with water shortage.
Landraces are defined as “a mixture of genotypes that evolved, largely
by natural selection, under the environmental conditions in which they
were grown” [[43]12]. They generally have a tolerance to abiotic
stresses and can produce reasonable yield under low-input farming
conditions. Landraces provide a valuable gene pool for the improvement
of bread wheat to be adapted to drought condition. It is imperative to
study stress-responsive genes in those landraces which evolved in
drought conditions [[44]13]. Iran is considered to be a part of the
center of origin [[45]14] and genetic diversity of bread wheat
(Triticum aestivum L.) [[46]15]. Evaluation of these landraces will
increase adaptability to drought, through the identification of genes
involved in this stress [[47]13]. Several mechanisms have been
identified including drought avoidance, escape or tolerance that plants
can overcome drought stress. It has been shown that plants respond to
drought stress with a wide range of genes that regulate several
metabolic pathways such as carbohydrate metabolism, lipid metabolism,
protein folding, secondary metabolic process and Mitogen-activated
protein kinase (MAPK) signaling pathway [[48]16, [49]17]. A recently
published review provided an overview of the candidate genes estimated
that participating in drought response in plants [[50]18]. Drought
avoidance, in particular, is a strategy in which plants maintain their
tissue water content relatively high, despite exposure to drought
stress. In this strategy, the relatively high tissue water content is
created through a deeper root system and increased hydraulic
conductance, etc [[51]19]. Despite numerous studies on drought stress
in plants, regulatory networks of transcripts in response to this
stress have not been completely identified, especially in the root
tissue.
The availability of the next generation sequencing technologies
provided high-throughput tools for studying genes expression profiles
at the level of whole genome [[52]20] that widely used for studying the
plant gene response to abiotic and biotic stresses [[53]21, [54]22]. A
variety of studies have been performed on the transcriptional responses
of wheat crop under drought stress [[55]22–[56]24]. Most of these
studies focused on gene expression profiles of above-ground organs or
well-known varieties, but not on the roots of wheat landraces.
Therefore, in this research, our main objectives were to identify the
genes participating in response to drought stress and investigate gene
regulation networks related to drought adaptation mechanisms in the
bread wheat landraces.
Materials and methods
Plant materials and drought treatment
Based on our previous survey of 123 Iranian bread wheat landraces,
“L-82”, a spring wheat landrace was selected for this experiment due to
its tolerance to drought stress. A two-stage greenhouse experimentation
was conducted to evaluate these landraces for drought tolerance. The
first part of the screening was executed based on the root and shoot
characteristics, relative water content (RWC) and parameters derived
from polyphasic chlorophyll fluorescence under drought and normal
conditions. To evaluate the plant materials, a 123×2 factorial
experiment in randomized complete block design with 4 replications was
applied. The GGE biplot function allowed us to choose vertex genotypes
close to effective traits as an indicator of water stress tolerance.
Based on the results, four landraces including ‘L-118’, one from each
sector, along with “Marvdasht” a well-known drought susceptible spring
variety were considered for the second step of the experiment. Although
the first part of the screening was done during the seedling stage of
plant growth, the second stage of screening was executed at the
reproductive stage. A 5 (genotype) × 2 (well-watered and drought
stress) factorial experiment in randomized block design with 2
replications and five seeds per pot was used in this stage of
experimentation. Similar to the first experiment, at the second stage,
screening of landraces was done based on the root and shoot
characteristics, RWC and chlorophyll fluorescence parameters under
contrasting water regimes. Besides, seed weight, seed number and
thousand kernel weight (TKW) were also applied for this stage of
screening. The results of these experiments showed a great variability
for all morphological and physiological traits among wheat genotypes in
the evaluated germplasm. This enabled us to screen genotypes that are
close to the well-known traits as an indicator of water stress
tolerance. It has been shown that genotype ‘L-82’ is located close to
most of the well-known traits as an indicator of water stress tolerance
(i.e., Performance index for the photochemical activity, TKW, seed
weight and root length). Therefore, this genotype was used for further
evaluation in this experiment. “Marvdasht” (HD2172/Bloudan//Azadi), a
well-known Iranian drought-susceptible variety of bread wheat [[57]25],
was used as a sensitive control in the experiment. The experiment was
conducted in greenhouse conditions (16 h daylight at 25 ± 3 °C and an 8
h dark period at 17 ± 3 °C) at Hamedan, Center of Agriculture and
Natural Resources Research (Iran). Each of the selected landraces was
planted (five seedlings per pot) in 120-cm deep pipes with 25.4 cm
diameter, filled with the same weight of field soil (33.6% sand, 34%
silt and 32.4% clay). Two watering treatments were (a) control, which
the soil moisture of the pots where maintained near to field capacity
(FC) as calculated by [58]Eq 1, and (b) drought stress, that the soil
water content was kept at 45% FC. Each growth experiment was done with
two replications. The field capacity (FC) and permanent wilting point
(PWP) of the soil were measured by applying 0.3 and 15 bar pressure on
the soil sample in the pressure chamber device as proposed by Abbott,
1985 [[59]26]. The drought treatment was started by withholding water
at the heading stage (Feekes’ growth stage 10.5) onwards. The soil
moisture for the pots of the well-watered and drought-stressed
conditions was maintained using the required amounts of water as the
benefit of the handy TRIME-FM IMKO GmbH, Germany, following the
manufacturer’s instructions. Pots on well-watered irrigation received
water on each day, as calculated by [60]Eq 1.
[MATH:
Dn=∑i=1
m[(θFCi-θ1i)×Bdi
×DiBd
w×100] :MATH]
(1)
in which, θ[FCi] represent the percent of moisture content at field
capacity, θ[1i] is the percent of moisture content before irrigation,
Bd[i] shows soil bulk density (g cm^−3), D[i] express soil depth (cm),
Bd[w] exhibit water bulk density (g cm^−3) and D[n] is equal to
irrigation water depth (cm). Pots on drought-stressed treatment
received only the corresponding fraction of the water provided to
well-watered irrigation, as described above.
A total of eight samples (two genotypes × two biological replications ×
two irrigation regimes) were used for cDNA library construction and
sequencing.
RNA extraction, cDNA synthesis and sequencing
For RNA extraction, the root tissues of all plants in each pot were
pooled and frozen immediately in liquid nitrogen and stored at -80°C.
Frozen tissues were ground with a mortar and pestle and approximately
100 mg of powdered tissue was sampled. RNA was isolated using the RNA
extraction kit (ZR plant RNA MiniPrep, Zymo Research, USA) according to
the manufacturer’s instructions. RNA integrity and purity were assessed
using a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific,
USA), an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara,
CA, USA) and by 1% agarose gel electrophoresis. RNA samples were sent
to Beijing Genomics Institute (BGI), China ([61]https://www.bgi.com)
for library construction and sequencing. Libraries were sequenced on an
Illumina HiSeq 2500 generating 2x150 bp paired-end reads.
DEGs identification and functional annotation
The sequenced reads were processed applying CLC Genomics Workbench 10.1
(CLC-BIO. Aarhus, Denmark). The raw reads were filtered by removing
adapter sequences, ambiguous nucleotides and low-quality sequences.
Trimming of reads were performed based on a minimum quality of Q27 and
a maximum of two ambiguous nucleotides. All downstream analyses were
based on trimmed data. Pearson’s correlation coefficients were applied
to quantify the correlation between biological replicates.
The high-quality reads were assembled using the wheat reference genome
sequences TGACv1 (Ensemble plants release 36) as a guide. The de novo
assembly tool within the CLC Genomic Workbench software allowing the
software to optimize the kmer size. The assembly was further refined by
mapping all the reads back (Mismatch cost = 2, insertion cost = 3,
deletion cost = 3, length overlap = 0.5 and similarity fraction = 0.8).
Reads from each individual sample were mapped to the transcriptome
assembly using the RNA-Seq tool in CLC with the same similarity
parameters used in the assembly phase detailed above. Empirical
analysis of differentially expressed genes (DEGs) developed by Robinson
and Smyth [[62]27] were employed in this study to find DEGs. The test
uses the raw counts and implicitly carries out normalization and
transformation of these counts. This test was carried out to select
differentially expressed genes by calculating the fold change of each
transcript for each genotype under stress condition in comparison to
their respective control. The genes having a minimum of 50 reads in at
least one sample type were removed prior to differential expression
analysis. The criteria to identify the putative differentially
expressed transcript were false discovery rate (FDR) ≤ 0.01 and fold
change (FC) ≥ 2 or ≤ -2. Statistical test was done between each
genotype under water stress compared to its control.
Gene ontology enrichment was carried out using the singular enrichment
analysis (SEA) function of the web-based tool AgriGO [[63]28]. The
input list consisted of different sets of DEGs, and annotation of wheat
genome (Ensemble plants release 36) was used as background.
Overrepresented terms in the three main categories, “Biological
Process”, “Molecular Function” and “Cellular Component” were filtered
using Fisher’s exact test and the Benjamini-Hochberg multiple testing
corrections (Q-value < 0.05). The resulting GO terms were plotted with
the Web Gene Ontology Annotation Plot (WEGO) tool [[64]29]. Pathway
enrichment analysis of DEGs was performed using the Kyoto Encyclopedia
of Genes and Genomes (KEGG, [65]http://www.genome.jp/kegg/). Heat maps
were generated using the web-based tool Morpheus
([66]https://software.broadinstitute.org/morpheus/).
RT-qPCR validation
To validate the RNA-Seq analysis, RT-qPCR was performed on a set of
eight randomly selected genes. DNase I treated RNA samples were reverse
transcribed using First Strand cDNA Synthesis Kit (Thermo Scientific).
Primers for RT-qPCR were designed using the Primer3Plus online software
([67]www.bioinformatics.nl/primers3plus) and their sequences are
available in [68]S1 Table. The RT-qPCR reactions were performed with
the MiniOpticon Real-Time PCR System, Bio-Rad (Hercules, CA, USA) and
the SYBR Green 2x Master Mix (Ampliqon) were used to detect transcript
abundance. The reaction was performed using 1.5μl of first strand cDNA,
3μl of each primer and 7.5 μl SYBR Green Master Mix in a final volume
of 15 μl. Negative control was also considered for each run. Cycling
programs were incubation at 95 °C for 5 min, then 40 cycles of
denaturation at 92 °C for 45 s, annealing at 60 °C for 45 s and
extension at 72 °C for 45 s. The specificity of all products was
verified via melting curve analysis by increasing the temperature from
60 °C to 95 °C and read every 0.5 °C. Three replicates were considered
for each gene. Normalization of reads was done concerning the Actin 2
as the reference gene. The relative quantitative method (2^-ΔΔCt) was
used to estimate quantitative gene expression.
Results
Morphological, physiological and seed yield-related traits of the evaluated
genotypes
In this study, a landrace selected from our previous study, “L-82”, was
used because of its tolerance to drought stress. The selected landrace
had a deep rooting system, high root weight and high root/shoot weight
ratio related to “Marvdasht”, the drought susceptible spring wheat
variety ([69]Table 1). There was a slight decrease in leaf RWC (20%)
after applying drought stress treatment, while the amount of RWC
decreased by ≈ 44% under the same situation. TKW, as a promising trait
for increasing grain yield in wheat, was reduced by 19% after drought
treatment in “L-82”. The reduction of this trait at the same water
availability was much lower (40%) in susceptible genotype. Among
parameters derived from polyphasic chlorophyll fluorescence, PIabs
(performance index for the photochemical activity) was a very sensitive
parameter to environmental stress. It has been shown that this
parameter was decreased by 11% in “L-82” after drought treatment,
whilst this decline was much higher (≈ 40%) in sensitive genotypes. The
value of each parameter in resistant landrace relative to the
corresponding value of the susceptible genotype (which thus become
equal to 100%) is plotted in [70]Fig 1. It is clear that root weight,
RWC, seed weight, TKW and PIabs have been decreased in the susceptible
genotype relative to resistant landrace, under drought treatment.
Table 1. Mean of morphological, physiological and seed yield-related traits
of the evaluated genotypes under control and stress conditions.
Genotype Condition Root length (cm) Root weight
(gr) Root/shoot weight ratio RWC
(%) Seed number Seed weight
(gr) TKW
(gr)
L-82 Normal 99±4[71]^* 305±46 0.936±0.222 87.81±1.14 929±119 34.36±4.36
36.99±0.45
Stress 100.5±2.5 197±11 1.57±0.047 69.62±4.06 790±20 23.64±0.66
29.93±0.14
Marvdasht Normal 52.5±7.5 91.5±0.77 0.464±0.0013 89.2±3.47 759±153
28.55±5.14 37.78±0.85
Stress 54±9 49±3.3 0.803±0.0737 50.15±1.87 612±98 13.77±1.46 22.70±1.25
Genotype Condition Fo [72]^a Fm [73]^b Ft [74]^c Fv [75]^d Fv/Fm [76]^e
PIabs [77]^f
L-82 Normal 1.0401±0.0097667 5.547±0.109 283.1283333±2.475
4.5069±0.09905 0.81249±0.0020777 2.1100168±0.28893
Stress 1.118017±0.44983 5.291±0.201 293.895±4.3516667 4.1806±0.135933
0.7907±0.00464 1.77301±0.04226
Marvdasht Normal 1.018±0.0281133 5.218±0.216 253.26833±43.515
4.1996±0.1873967 0.8049±0.0024502 2.262912±0.0822166
Stress 1.1027±0.042933 4.917±0.143 253.595±45.93833 3.8163±0.079633
0.7761±0.006345 1.3665943±0.1025208
[78]Open in a new tab
* Mean of values ± standard errors of the respective traits
^aFo is the initial/minimal fluorescence;
^bFm is the maximal fluorescence value;
^cFt is the fluorescence intensity at time t;
^dFv is variable fluorescence, Fv = Fm–Fo;
^eFv/Fm represents the maximum quantum yield of Photosystem II (PSII);
and
^fPIabs is the performance index for the photochemical activity
[[79]30].
Fig 1. The value of each parameter in resistant landrace relative to the
corresponding value of the susceptible genotype.
[80]Fig 1
[81]Open in a new tab
Fo is the initial/minimal fluorescence; Fm is the maximal fluorescence
value; Ft is the fluorescence intensity at time t; Fv is variable
fluorescence, Fv = Fm–Fo; Fv/Fm represents the maximum quantum yield of
Photosystem II (PSII); and PIabs is the performance index for the
photochemical activity [[82]30].
Transcriptomics profiles and analysis of DEGs
As roots are the first organs to be exposed to drought, this tissue was
sampled from drought-tolerant and drought-sensitive genotypes under
well-watered and drought-stressed conditions. The aim was to achieve a
broad survey of genes associated with drought stress in bread wheat
landraces under water scarcity. Sequencing and trimming yielded 86 and
113.14 million trimmed reads from tolerant (ED, L-82 under drought and
EN, L-82 under normal conditions) libraries. Likewise, 121.83 and
106.55 million trimmed reads were generated from sensitive libraries
(MD, Marvdasht under drought and MN, Marvdasht under normal
conditions). The sequence reads of this study are available in the NCBI
Sequence Read Archive under accession number SRP140591.
These reads were assembled into 118 918 transcripts representing 97 918
gene loci. Nearly, 79 percent of the reported transcripts for these
libraries had a length from 500 to 1900 bp ([83]S1 Fig). It was found
that 106,681 transcripts expressed continuously in all samples. The
number of expressed transcripts varied from 69.1 percent to 72.8
percent in the samples ([84]S2 Fig). Studying of differential gene
expressions was done by calculating the fold change of each transcript
for each genotype under stress condition in comparison to their
respective control. The criteria to identify the putative
differentially expressed transcripts were false discovery rate (FDR)
p-value less than 0.01 and fold change (FC) ≥ 2 or ≤ -2. In total,
3,840 transcripts representing 1,690 unique gene loci were
differentially expressed under the various experimental conditions.
Among these differentially expressed gene loci, 460 (43 up-regulated
and 417 down-regulated) and 236 (106 up-regulated and 130
down-regulated) genes were differentially expressed in ED vs. EN and MD
vs. MN, respectively. The number of DEGs showing overlapping and
specific response under various experimental conditions is plotted in
[85]Fig 2. Overall, only a small number of genes were found to be
common in different comparisons between genotypes.
Fig 2. The number of DEGs (≥ 2-fold change; p ≤ 0.01) showing overlapping and
specific response under various experimental conditions in drought.
[86]Fig 2
[87]Open in a new tab
ED, L-82 under drought and EN, L-82 under normal conditions. MD,
Marvdasht under drought and MN, Marvdasht under normal conditions. ↑
and ↓, up and down-regulated DEGs.
Functional classification of DEGs
The singular enrichment analysis (SEA), carried out with AgriGO
software [[88]28] on the 419 and 195 unique DEGs against genome
reference in tolerant and susceptible genotypes, highlighted 115 and 12
GO terms, respectively ([89]S3 and [90]S4 Figs). GO term enrichment
analysis of 37 up and 382 down-regulated unique DEGs of tolerant
genotype resulted in 1 and 110 significantly enriched GO terms ([91]S5
Fig). “Respond to stress” (GO:0006950) in the biological process (BP)
category was the only GO term enriched in up-regulated DEGs in tolerant
genotype. On the other hand, “secondary metabolic process”
(GO:0019748), “cell volume homeostasis” (GO:0006884), “glycerol
transport” (GO:0015793), “cellular water homeostasis” (GO:0009992),
“polyol transport” (GO:0015791), “water transport” (GO:0006833),
“organic hydroxy compound transport” (GO:0015850), “fluid transport”
(GO:0042044), “water homeostasis” (GO:0030104) and “secondary
metabolite biosynthetic process” (GO:0044550) in the BP category,
“glycerol channel activity” (GO:0015254), “glycerol transmembrane
transporter activity” (GO:0015168) and “O-methyltransferase activity”
(GO:0008171) in the molecular function (MF) category and finally
“intrinsic component of plasma membrane” (GO:0031226), “integral
component of plasma membrane” (GO:0005887) and “mitochondrion”
(GO:0005739) in the cellular components (CC) category were the most
enriched GO terms among the 382 uniquely down-regulated DEGs in
tolerant genotype.
GO term enrichment analysis of 99 up-regulated DEGs of susceptible
genotype highlighted 13 significantly enriched GO terms but no GO term
was significantly enriched among 96 down-regulated unique DEGs of
susceptible genotype ([92]S6 Fig). Calcium ion binding (GO:0005509) and
hydrolase activity, acting on glycosyl bonds (GO:0016798) in the
molecular function category and plasma membrane (GO:0005886) and
cytoplasm(GO:0005737) in the cellular component category were some of
the significantly enriched GO terms in up-regulated DEGs of susceptible
genotype ([93]S6 Fig).
The entire set of DEGs was subjected to GO analysis to obtain deep
functional characterization. Altogether, 419 unique DEGs (37 up and 382
down-regulated) in tolerant and 195 unique DEGs (99 up and 96
down-regulated) in susceptible genotypes were divided into 39
subcategories within three main categories including biological process
(19), molecular function (9) and cellular components (11), in GO level
2. In total, 139 and 41 DEGs in the tolerant and susceptible genotypes
were associated with “cellular component” terms ([94]Fig 3). The number
of DEGs associated with “molecular function” terms were 251 for
tolerant and 60 for susceptible genotypes. For “biological process”,
the number of associated DEGs was 182 and 51 for tolerant and
susceptible genotypes ([95]Fig 3). In both genotypes, “metabolic
process”, “cellular process”, “biological regulation” and “respond to
stimulus” were the most represented “biological process” subcategories.
As for “molecular function”, the major subcategories were “binding” and
“catalytic activity”, followed by “transporter activity” and
“antioxidant activity”. “Cell” and “cell part” followed by “organelle”
and “membrane” were the dominant subcategories in “cellular components”
term ([96]Fig 3). In addition, it was clear that “immune system
process” (GO:0002376), “growth” (GO:0040007), “cell killing”
(GO:0001906), detoxification (GO:0098754) and “negative regulation of
biological process” (GO:0009892) in “biological process” term were
associated with tolerant genotype DEGs ([97]Fig 3). Besides,
“transcription regulator” (GO:0030528), “molecular function regulator”
(GO:0098772), “nutrient reservoir activity” (GO:0045735) and “signal
transducer activity” (GO:0004871) in MF term and “membrane-enclosed
lumen” (GO:0031974) in CC term were related to tolerant genotype DEGs
([98]Fig 3).
Fig 3. Histogram of GO terms assigned to DEGs in roots of “L-82” (n = 419)
and “Marvdasht” (n = 195) in black and purple column, respectively.
[99]Fig 3
[100]Open in a new tab
The DEGs are categorized into three main groups: cellular component
(CC), molecular function (MF), and biological processes (BPs). ED, L-82
under drought and EN, L-82 under normal conditions. MD, Marvdasht under
drought and MN, Marvdasht under normal conditions.
The most represented pathways
Having understood the functions of DEGs, we mapped all non-redundant
DEGs to terms in the KEGG database and found that 47 and 46 pathways
(77 KEGG pathway items) were enriched in ED vs. EN and MD vs. MN,
respectively ([101]Fig 4). These different pathways were assigned to 16
clades under four major KEGG categories, namely, “metabolism”, “genetic
information processing”, “cellular processes” and “environmental
information processing” ([102]S7 Fig). “Carbohydrate metabolism”,
“signal transduction” and “xenobiotics biodegradation and metabolism”
were the top three up-regulated pathways represented by the uniquely
expressed genes in ED vs. EN. On the other hand, the most enriched
down-regulated pathways were “biosynthesis of other secondary
metabolites”, “amino acid metabolism” and “xenobiotics biodegradation
and metabolism” ([103]S7 Fig).
Fig 4. The over-represented unigenes from differential gene expression data
of tolerant and sensitive libraries and their KEGG terms.
Fig 4
[104]Open in a new tab
ED, L-82 under drought and EN, L-82 under normal conditions. MD,
Marvdasht under drought and MN, Marvdasht under normal conditions.
The most enriched metabolic pathways in up-regulated DEGs of ED vs. EN
were “amino sugar and nucleotide sugar metabolism”, “glyoxylate and
dicarboxylate metabolism”, “MAPK signaling pathway–plant”,
“phenylpropanoid biosynthesis”, “protein processing in endoplasmic
reticulum” and “tryptophan metabolism”. In addition, down-regulated
DEGs were enriched mostly to “phenylpropanoid biosynthesis”,
“ubiquinone and other terpenoid-quinone biosynthesis”, “stilbenoid,
diarylheptanoid and gingerol biosynthesis” and “phenylalanine
metabolism” pathways. Additionally, most of the up-regulated DEGs of MD
vs. MN were sorted to “arginine and proline metabolism, “drug
metabolism-cytochrome P450”, “galactose metabolism”, “glutathione
metabolism” and “metabolism of xenobiotics by cytochrome P450”
pathways. On the other hand, top enriched pathways of down-regulated
DEGs were “stilbenoid, diarylheptanoid and gingerol biosynthesis”,
“oxidative phosphorylation” and “amino sugar and nucleotide sugar
metabolism”.
All unique DEGs across different comparisons were clustered into six
clades according to their expression patterns ([105]S8 Fig). The genes
in cluster I (382) down-regulated in drought-tolerant genotype while
having no differential expression in control plants. On the other hand,
cluster II had genes (96) that down-regulated in control plants but had
no differential expression in drought-tolerant genotypes. Cluster III
(35) and IV (6) comprised down and up-regulated genes with almost the
same expression profile across the comparisons. Cluster V (37)
exhibited DEGs which were up-regulated in drought-tolerant genotype
while they had no expression in control plants. Reversely, the last
cluster (VI) showed genes (99) that were up-regulated in control plants
did not have not significant expression level in drought-tolerant
genotype. Assessing each cluster for related GO terms showed that genes
involved in “methylation”, “phenylpropanoid biosynthetic process”,
“cellular detoxification” and “cinnamic acid biosynthetic process” in
the biological process category and “phenylalanine ammonia-lyase
activity”, “O-methyltransferase activity” and “protein dimerization
activity” in the molecular function category were enriched in cluster 1
genes, which were down-regulated in drought-tolerant genotype under
drought stress condition. However, genes involved in “transport”,
“localization” and “response to stress” from BP category and
“transporter activity” and “hydrolase activity” in MF category were of
enriched GO terms in cluster 2, showing genes that down-regulated in
control plants and had no DEGs in tolerant genotypes. No GO term was
enriched in cluster III, IV and VI, however, “response to stimulus”,
“response to stress”, “single-organism metabolic process” and
“single-organism process” GO terms were enriched in cluster V, in which
genes upregulated significantly in drought-tolerant genotype. The KEGG
pathway analysis of each cluster was performed to assign the related
biological pathways for each group of DEGs. Genes located in cluster
I-VI were assigned to 31, 5, 2, 5, 16 and 37 KEGG pathways,
respectively ([106]S9 Fig).
The most up and down-regulated genes in both genotypes
Analysis of DEGs in drought treated genotypes along with their controls
will assist our understanding of the molecular events involved in
drought stress response. Therefore, significantly up and down-regulated
DEGs in each set of comparisons were considered in this study. Out of
the 417 down-regulated genes in ED vs EN, only one gene was in common
with genes up-regulated in MD vs MN; however, 34 genes were in common
with genes down-regulated in MD vs MN ([107]Fig 2). Venn diagram also
showed that only 6 genes out of 43 genes were common between the
up-regulated genes of ED vs EN and MD vs MN comparisons ([108]Fig 2).
The top up-regulated DEGs in tolerant genotype under stress condition
in comparison to its respective control were annotated to genes
functional in secondary metabolite biosynthesis, oxidation-reduction
process, heat shock proteins and chaperones, cell wall modification,
transferase activity, co-factor and vitamin metabolism, DNA binding,
abiotic stress response, calcium-binding domain and transmembrane
transport signaling ([109]Table 2). The most significantly
down-regulated transcripts in tolerant genotype were annotated to
Auxin-responsive proteins, Secondary Metabolite Biosynthesis, Ion
binding, DNA binding, Protein degradation and signaling ([110]Table 3).
Table 2. List of specific genes which are significantly up-regulated in ED
vs. EN.
Gene ID Fold Change Annotation Function
TRIAE_CS42_7AL_TGACv1_556333_AA1761120 6.29 Cinnamyl alcohol
dehydrogenase Secondary Metabolite Biosynthesis
TRIAE_CS42_1AS_TGACv1_020323_AA0076650 3.25 Cytochrome P450
Oxidation-reduction process
TRIAE_CS42_1BS_TGACv1_049748_AA0160770 2.45 Cytochrome P450
Oxidation-reduction process
TRIAE_CS42_6BS_TGACv1_513206_AA1634480 2.31 Catalase isozyme 2
Oxidation-reduction process
TRIAE_CS42_6BS_TGACv1_513206_AA1634470 2.69 Catalase isozyme 2
Oxidation-reduction process
TRIAE_CS42_4DS_TGACv1_362595_AA1180930 2.31 Heat shock protein 70 Heat
shock proteins and chaperones
TRIAE_CS42_5DL_TGACv1_433530_AA1415640 3.45 UDP-glucose 4-epimerase
Cell wall modification
TRIAE_CS42_5BL_TGACv1_404728_AA1309330 3.16 Putrescine hydroxycinnamoyl
transferase-like Transferase activity
TRIAE_CS42_7DL_TGACv1_605357_AA2006220 12.15
Dolichyl-diphosphooligosaccharide—glycosyltransferase 48 kDa subunit
Transferase activity
TRIAE_CS42_3AS_TGACv1_210771_AA0678600 3.12 Fatty acyl-CoA reductase
Lipid biosynthesis
TRIAE_CS42_4DS_TGACv1_361374_AA1166440 2.78 Phosphomethylpyrimidine
chloroplastic Co-factor and vitamin metabolism
TRIAE_CS42_4AS_TGACv1_306585_AA1010680 2.91 B3 DNA binding domain DNA
binding
TRIAE_CS42_5BS_TGACv1_423362_AA1375120 2.43 Osmotin/thaumatin-like
protein Abiotic stress response
TRIAE_CS42_4AL_TGACv1_289983_AA0979890 2.13 Osmotin/thaumatin-like
protein Abiotic stress response
TRIAE_CS42_7BS_TGACv1_593723_AA1954050 3.03 Caleosin-related
Calcium-binding domain
TRIAE_CS42_6BS_TGACv1_513535_AA1644380 2.94 Bidirectional sugar
transporter SWEET14-like Probably SWEET sugar transporter
TRIAE_CS42_6AL_TGACv1_471216_AA1504800 2.14 Mechanosensitive ion
channel Transmembrane transport signaling
TRIAE_CS42_2DL_TGACv1_161880_AA0560970 5.24 Uncharacterized protein
Uncharacterized function
TRIAE_CS42_3AL_TGACv1_193635_AA0615920 4.08 Uncharacterized protein
Uncharacterized function
TRIAE_CS42_3DL_TGACv1_250209_AA0864280 2.75 Chitinase Cell wall
macromolecule catabolic process
[111]Open in a new tab
Table 3. List of specific genes which are significantly down-regulated in ED
vs. EN.
Gene ID Fold Change Annotation Function
TRIAE_CS42_7DL_TGACv1_603709_AA1987720 -3.97 AUX/IAA protein
Auxin-responsive protein
TRIAE_CS42_1AL_TGACv1_001992_AA0037540 -2.04 Small auxin-up RNA
Auxin-responsive protein
TRIAE_CS42_1AL_TGACv1_002495_AA0042400 -8.02 Probable methyltransferase
Methyltransferase activity
TRIAE_CS42_1AS_TGACv1_019041_AA0058780 -20.29 Histone -like DNA binding
TRIAE_CS42_2DL_TGACv1_158781_AA0526060 -4.73 Phenylalanine
ammonia-lyase Secondary Metabolite Biosynthesis
TRIAE_CS42_2AL_TGACv1_093330_AA0277780 -4.11 Phenylalanine
ammonia-lyase Secondary Metabolite Biosynthesis
TRIAE_CS42_2BL_TGACv1_129397_AA0381780 -3.17 Glycoside hydrolase
Secondary Metabolite Biosynthesis
TRIAE_CS42_4AL_TGACv1_288958_AA0962110 -3.49 Dirigent protein Secondary
Metabolite Biosynthesis
TRIAE_CS42_6BL_TGACv1_499669_AA1588820 -4.43 Serine threonine- kinase
Protein kinase activity
TRIAE_CS42_7DL_TGACv1_602629_AA1963150 -4.36 Plant peroxidase
Oxidation-reduction process
TRIAE_CS42_5AL_TGACv1_376349_AA1235590 4.04 Carotenoid hydroxylase
Oxidation-reduction process
TRIAE_CS42_2DL_TGACv1_159781_AA0542640 -5.11 Cellulose synthase
Cellulose biosynthetic process
TRIAE_CS42_1BL_TGACv1_032207_AA0126810 -4.77 Phosphoethanolamine
N-methyltransferase 3 Lipid metabolism
TRIAE_CS42_5DL_TGACv1_434151_AA1430490 -3.37 NAD-dependent
epimerase/dehydratase Amino acid metabolism
TRIAE_CS42_1DS_TGACv1_080223_AA0243550 -8.07 Cysteine protease Protein
degradation
TRIAE_CS42_2BS_TGACv1_147844_AA0488140 -4.47 Aspartic peptidase Protein
degradation
TRIAE_CS42_4DL_TGACv1_342374_AA1111780 -45.93 Serine/threonine-protein
kinase Signaling
TRIAE_CS42_1DS_TGACv1_080147_AA0241250 -4.03 E3 ubiquitin- ligase
RDUF1-like Signaling
TRIAE_CS42_5DS_TGACv1_456816_AA1478520 -4.66 Aquaporin Transport
TRIAE_CS42_5BL_TGACv1_405562_AA1330380 -3.28 NRT1 PTR FAMILY -like
Transport
TRIAE_CS42_1BL_TGACv1_030599_AA0095340 -4.07 Protein of unknown
function DUF538 Unknown function
TRIAE_CS42_2BL_TGACv1_129508_AA0386470 -4.07 WAT1-related protein
Presumably transport
TRIAE_CS42_7DL_TGACv1_603849_AA1990130 -2.82 Benzyl alcohol
O-benzoyltransferase-like Transferase activity
TRIAE_CS42_5AL_TGACv1_375844_AA1227900 -17.16 Unknown Not specified
TRIAE_CS42_3DL_TGACv1_251387_AA0881030 -2.09 Glutaredoxin-C2 Redox
activity
TRIAE_CS42_2AS_TGACv1_112139_AA0331210 -2.03 Heavy metal-associated
isoprenylated plant 23-like Metal ion binding
TRIAE_CS42_7DL_TGACv1_604406_AA1997870 -10.55 Germin Ion binding
TRIAE_CS42_2DL_TGACv1_159589_AA0540300 -2.79 Dehydrogenase/reductase
Oxidoreductases
[112]Open in a new tab
When sensitive genotypes under drought stress were compared to control,
the top upregulated DEGs were annotated to stress-responsive genes,
calcium ion binding, protein degradation and transcription ([113]Table
4). On the other hand, down-regulated DEGs in sensitive genotype were
annotated to genes functional in transcription, transport, ion binding,
protein degradation, lipid metabolic process and oxidation-reduction
process ([114]Table 5).
Table 4. List of specific genes which are significantly up-regulated in MD
vs. MN.
Gene ID Fold Change Annotation Function
TRIAE_CS42_7DS_TGACv1_625862_AA2065960 5.3 MYB transcription factor
Transcription
TRIAE_CS42_4BS_TGACv1_327885_AA1077740 4.06 Phosphatidylinositol
4-phosphate 5-kinase Protein degradation
TRIAE_CS42_2BL_TGACv1_130527_AA0412950 6.04 EF-hand domain Calcium ion
binding
TRIAE_CS42_3AL_TGACv1_195082_AA0644480 3.82 EF-hand domain Calcium ion
binding
TRIAE_CS42_1DS_TGACv1_080885_AA0255140 8.73 Glycosyltransferase
Transferase activity
TRIAE_CS42_1BS_TGACv1_049649_AA0158860 10.89 AWPM-19 Response to stress
TRIAE_CS42_2BL_TGACv1_130567_AA0413710 16.51 Late embryogenesis
abundant proteins Response to stress
TRIAE_CS42_4AS_TGACv1_308636_AA1028740 2.89 Glutathione S-transferase
Oxidoreductase activity
TRIAE_CS42_4BS_TGACv1_327885_AA1077740 4.06 Phosphatidylinositol
4-phosphate 5-kinase Signaling
TRIAE_CS42_1AS_TGACv1_019260_AA0064010 3.35 Lipase PLAT LH2 Protein
binding
TRIAE_CS42_1AL_TGACv1_000553_AA0014560 4.37 Choline-phosphate
cytidylyltransferase 2 Catalytic activity
TRIAE_CS42_3B_TGACv1_221059_AA0728710 5.84 Heat shock protein 70 Heat
shock proteins and chaperons
TRIAE_CS42_6DL_TGACv1_526901_AA1694680 2.37 Bidirectional sugar
transporter SWEET Transport
TRIAE_CS42_2BL_TGACv1_129668_AA0392100 2.07 Niemann-Pick C1 -like Lipid
transporter activity
TRIAE_CS42_3AS_TGACv1_210771_AA0678600 2.1 Fatty acyl-CoA reductase
Lipid biosynthesis
TRIAE_CS42_6AL_TGACv1_472258_AA1520010 2.37 Annexin Calcium ion binding
TRIAE_CS42_2AS_TGACv1_113656_AA0359020 2.37 Unknown Not specified
TRIAE_CS42_1BL_TGACv1_032777_AA0134000 3.55 Uncharacterised protein
family, basic secretory protein Are believed to be part of the plants
defense mechanism
TRIAE_CS42_6DS_TGACv1_544329_AA1747800 16.53 Not specified Not
specified
[115]Open in a new tab
Table 5. List of specific genes which are significantly down-regulated in MD
vs. MN.
Gene ID Fold Change Annotation Function
TRIAE_CS42_6BL_TGACv1_499354_AA1579060 -4.56 bHLH transcription factor
Transcription
TRIAE_CS42_2DL_TGACv1_158876_AA0527740 -3.38 NAC transcription factor
Transcription
TRIAE_CS42_4BS_TGACv1_328328_AA1086300 -9.15 Zinc finger Transcription
TRIAE_CS42_2AL_TGACv1_093527_AA0281800 -6.65 NAC transcription factor
Transcription
TRIAE_CS42_2BS_TGACv1_146827_AA0473530 -6.42 Protein kinase domain
Signaling
TRIAE_CS42_1DL_TGACv1_061111_AA0185720 -7.92 Subtilisin-like protease
Protein degradation
TRIAE_CS42_2DL_TGACv1_159703_AA0541750 -5.66 Serine-type peptidase
Protein degradation
TRIAE_CS42_2DS_TGACv1_177710_AA0582890 -2.65 UDP-glycosyltransferase
Transferase activity
TRIAE_CS42_6AL_TGACv1_470863_AA1497090 -19.18 Aquaporin-like Transport
TRIAE_CS42_7AS_TGACv1_570003_AA1828470 -15.8 Phosphate transporter
Transport
TRIAE_CS42_3B_TGACv1_221341_AA0737900 -3.7 NADP-dependent malic
chloroplastic Metabolism
TRIAE_CS42_4BL_TGACv1_320461_AA1040130 -3.21 phosphatase 2C 48 Cation
binding
TRIAE_CS42_3DL_TGACv1_251103_AA0877420 -9.14 Glycerophosphodiester
phosphodiesterase GDPD2 Lipid metabolic process
TRIAE_CS42_1AS_TGACv1_019480_AA0067160 -5.71 Peroxidase
Oxidation-reduction process
TRIAE_CS42_4BL_TGACv1_322934_AA1073440 -9 Bidirectional sugar
transporter SWEET12-like Transport
TRIAE_CS42_2BS_TGACv1_146508_AA0467110 -3.96 WAT1-related
At5g64700-like Presumably transport
TRIAE_CS42_2BL_TGACv1_134090_AA0443680 -11.62 Germin Ion binding
TRIAE_CS42_U_TGACv1_641430_AA2094940 -17.46 Phosphatase phospho1
Phosphatase activity
[116]Open in a new tab
Distribution of DEGs across the wheat genome
The distribution of DEGs across the wheat genome showed that more
up-regulated genes in drought-tolerant genotype were located in the B
genome ([117]Fig 5). However, it is obvious that all chromosomes except
2A, 4B and 6D had contributed in response to drought. In addition, most
of the DEGs mapped to chromosome D were down-regulated in tolerant
genotype in response to drought. On the other hand, in control plants,
the most up-regulated genes were mapped to B genome while
down-regulated genes in this genotype mostly occurred in A genome
([118]Fig 5).
Fig 5. Distribution of DEGs across the wheat genome.
[119]Fig 5
[120]Open in a new tab
Number of DEGs that are up or down-regulated in drought tolerance and
drought sensitive genotype roots are shown in genome A, B and D. ED,
L-82 under drought and EN, L-82 under normal conditions. MD, Marvdasht
under drought and MN, Marvdasht under normal conditions.
RT-qPCR validation of differential gene expression
To confirm the results of the RNA-Seq analysis, a total of eight
randomly selected genes (four genes for each set of comparisons),
phosphomethylpyrimidine synthase, cytochrome P450 family 76 subfamily
C, Niemann-Pick C1 protein, choline-phosphate cytidylyltransferase,
beta-carotene 3-hydroxylase, benzyl alcohol O-benzoyltransferase,
pyruvate dehydrogenase phosphatase and
phosphoethanolamine/phosphocholine phosphatase ([121]S1 Table), were
considered for RT-qPCR validation ([122]Fig 6a). The expression level
(up and down-regulated) of selected genes for qRT-PCR was the same as
that in RNA-Seq except for phosphoethanolamine/phosphocholine
phosphatase. However, a similar trend was apparent in the expression of
this gene in both assays, but the level of expression in qRT-PCR was ≈
¼ of that in RNA-Seq. Totally, the results illustrated that the
expression data of the selected DEGs in qRT-PCR correlate highly to
that in RNA-Seq analysis ([123]Fig 6b and 6c), indicating the
reliability of the transcriptomics profiling data.
Fig 6. Expression patterns of 8 selected genes in “L-82” and “Marvdasht” in
contrasting water regimes determined by RNA-Seq and qPCR.
[124]Fig 6
[125]Open in a new tab
RNA-Seq values represent the ratio of the expression level in drought
treated genotypes to the expression level of their related controls.
For qPCR, the data are relative expression ± SD from three independent
replicates (a). The qPCR primers for each contig are listed in
additional file 1. Linear regression containing the RNA-Seq log[2]
value and the RT-qPCR validation date, for drought tolerance (b) and
drought sensitive (b) genotypes. ED, L-82 under drought and EN, L-82
under normal conditions. MD, Marvdasht under drought and MN, Marvdasht
under normal conditions.
Discussion
Exposure to drought elicits a range of responses in plants that involve
a large number of pathways related to different mechanisms [[126]16,
[127]31]. In this experiment, we evaluated differential expression of
genes related to drought stress in the root tissue of wheat crop. The
responses of sensitive genotype (Marvdasht) were different from
tolerant landrace (L-82) in case of gene expression level and pathways
involved in response to drought. However, the number of DEGs was higher
in tolerant genotype but most of these genes (≈ 90%) have been
down-regulated ([128]Fig 2). On the other hand, the number of up and
down-regulated genes is almost equal in sensitive genotype. Based on
this finding, we hypothesize that the reaction of “L-82” happened at
the very initial stages of drought stress [[129]11]. However, the
up-regulation of genes related to “secondary metabolic process” (GO:
0019748) under well-watered condition in the “L-82” affirms its
intrinsic tolerance against stresses [[130]7, [131]32]. The only
enriched GO term in up-regulated DEGs of drought-tolerant genotype was
response to stress (GO:0006950), mostly oxidative-reduction to protect
the plants against oxidative stress [[132]10, [133]11, [134]33]. The
“first line of defense” against ROS and oxidative stress in plant cell
is hypothesized to be antioxidant enzymes [[135]32, [136]34].
Transcriptomic profiling of tolerant landrace revealed avoidance as the
dominant mechanisms in “L-82” against drought stress
Up-regulation of genes participating in phenylpropanoids
(TRIAE_CS42_7AL_TGACv1_556333_AA1761120) and flavonoids
(TRIAE_CS42_1AS_TGACv1_020323_AA0076650 and
TRIAE_CS42_1BS_TGACv1_049748_AA0160770) pathways was observed in L-82
landrace after drought treatment. Genes coding for these secondary
metabolites have the potential to protect plants against oxidative
stress [[137]34–[138]36]. The phenylpropanoid pathway involved in the
biosynthesis of different products such as lignin is critical in plant
adaptation to environmental adversities [[139]37] such as drought
stress [[140]38]. Reinforcement of cell walls for mechanical strength
and for water conductance will occur with increasing of lignin content
[[141]38, [142]39]. On the other hand, down-regulation of lignin
production as a product of phenylpropanoid compound
(TRIAE_CS42_2DL_TGACv1_158781_AA0526060 and
TRIAE_CS42_2AL_TGACv1_093330_AA0277780) may stimulate growing of roots
[[143]40]. It means that not only is the ability of plants to access
water from deep under the soil important [[144]19], but the roots also
need to have the adequate hydraulic characteristics to allow water
uptake conditions of scarcity [[145]41]. This could be due to a
differential cellular response of root tissues to conform the shape of
the root as suggested by previous studies [[146]42]. The importance of
scavenging of ROS generated during drought stress has been suggested as
the reason for up-regulation in flavonoids biosynthesis [[147]43,
[148]44] and these may play a role as antioxidants in “L-82”. Although
there are contrasting reports about up or down-regulation of flavonoids
in susceptible and tolerant genotypes [[149]32, [150]43], it has been
shown that the accumulation of these enzymes can be specific to
species, genotype, stress duration and severity [[151]32, [152]45].
Having considered the above results, we hypothesize that the drought
tolerance in “L-82” could be due to avoidance mechanisms like
optimization of water uptake [[153]46]. Morphological assessments of
this genotype under drought condition also revealed that this genotype
produced strong root system under drought condition ([154]Table 1).
However, some other aspects of drought-tolerant mechanisms like
oxidative stress defense (TRIAE_CS42_6BS_TGACv1_513206_AA1634480 and
TRIAE_CS42_6BS_TGACv1_513206_AA1634470) are obvious in “L-82” against
drought stress [[155]47]. Considering avoidance as the dominant
mechanisms in “L-82” against drought stress, we speculate that the high
number of down-regulated genes in this resistance line is related to
the management of energy resources driving growth and development in
response to environmental perturbations [[156]48].
Down-regulation of genes associated with indole-3-acetic acid (IAA)
(TRIAE_CS42_7DL_TGACv1_603709_AA1987720 and
TRIAE_CS42_1AL_TGACv1_001992_AA0037540) was also observed under drought
condition in “L-82”. It has generally been shown that auxin has a
negative effect on drought adaptation in plants [[157]19].
Down-regulation of IAA was shown to be associated with drought
tolerance in plants [[158]49]. A QTL study also showed that the
up-regulation of the DEEPER ROOTING gene in rice resulted in drought
avoidance and was negatively regulated by auxin [[159]50].
Pathway analysis also showed that thiamine metabolism was up-regulated
in “L-82” (TRIAE_CS42_4DS_TGACv1_361374_AA1166440) as a result of
drought treatment, which was validated by RT-qPCR (THIC). Accumulation
of thiamin in plants subjected to abiotic stress condition is reported
by several studies [[160]51, [161]52]. Interestingly, it is proved that
exogenous thiamin treatment in Arabidopsis improves root growth of
seedlings subjected to oxidative stress [[162]53]. It is obvious that
most of the mechanisms conferring drought resistance to “L-82” are
related to root growth enhancement in which more access to water is
possible through root penetration deep into the soil.
Drought stress causes endoplasmic reticulum (ER) protein folding
machinery, which reaches a limit and the demand for protein folding
exceeds its capacity. Increasing unfolded or miss-folded proteins in
the ER triggers an unfolded protein response. This results in
up-regulating the expression of genes encoding components of protein
folding machinery or the ER-associated degradation system [[163]54].
The same was detected in “L-82”
(TRIAE_CS42_4DS_TGACv1_362595_AA1180930) as the expression of genes
encoding heat shock protein 70 which is the most abundant chaperone
protein in the ER and oligosaccharyl transferase complex was
significantly up-regulated [[164]54].
Elevated root suberin (TRIAE_CS42_3AS_TGACv1_210771_AA0678600) is also
activated in both drought-treated genotypes. However, the fold change
in “L-82” was almost 1.5 times more than that in sensitive genotype. It
has been shown that root suberin in Arabidopsis is associated with “a
root-dependent increase in time to wilting during water stress”
[[165]55]. The researchers supposed that resistance to wilting could be
achieved through a reduction in water loss to the environment or via an
increase in the ability to take up water from the soil. They also
established that this delayed wilting was a root-dependent phenomenon
[[166]55].
Up-regulation of chitinase (TRIAE_CS42_3DL_TGACv1_250209_AA0864280) and
UDP-glucose 4-epimerase (UGE) (TRIAE_CS42_5DL_TGACv1_433530_AA1415640)
as two enzymes participating in “amino sugar and nucleotide sugar
metabolism” pathway happened in “L-82” under drought stress. There are
several observations suggesting that the UGE enzymes play a role in
plant root development and/or growth [[167]56, [168]57]. A possible
function of UGE in plant root growth after its increased expression
under drought stress is root thickness [[169]56]. In plants, it has
been observed that levels of chitinases are regulated by biotic and
abiotic stress such as pathogen infection, drought, cold, salt, heavy
metals and plant hormones [[170]58, [171]59]. The role of chitinases
has been also studied in grasses such as rye in response to cold and
drought stress [[172]60]. It has been shown that chitinases induced in
tomato plants make them tolerant to drought in comparison to the
susceptible genotype [[173]61].
B and D-genomes had a greater contribution in response to drought stress than
A-genome
In our study, most of the up-regulated genes in “L-82” belong to the B
genome ([174]Fig 5). Diploid genomes A and B (what is known belong to
Triticum urartu and Aegilops speltoides) formed the tetraploid wheat
Triticum turgidum through an alloploydization event during the
evolution processes. Following this step, with the contribution of the
D genome (Aegilops tauschii), hexaploid wheat (Triticum aestivum)
appeared [[175]62]. However, bread wheat has a hexaploid genome but it
tends towards functional diploidy [[176]63]. Moreover, it has been
shown that B-genome homoeoloci have more tendency for contribution in
gene expression than the other genome homoeoloci [[177]63]. In another
study, it has been revealed that chromosomes 3B, 5B and 2B have more
transcripts in both roots and leaves in response to drought stress
[[178]64]. However, most of the down-regulated genes in “L-82” belong
to D genome but there are some genes that up-regulated significantly in
this genome in response to drought stress, especially 3D, 4D and 5D
chromosome. It have been shown that D genome of Aegilops tauschii has
the potential of increasing drought adaptation in hexaploid wheat
[[179]65, [180]66]. Therefore, despite the occurrence of more
down-regulated genes of this genome in root tissue of drought-treated
landrace than in the other genomes, we observe that the genes that
contributed in drought tolerance/avoidance belong to this genome.
In summary, this study utilized RNA-Seq data in bread wheat genotypes
with contrasting drought tolerance appearance. In this study, an
Iranian bread wheat landrace was used which is highly adapted to the
adverse environmental condition. A review of the differential expressed
genes revealed that this landrace uses the avoidance mechanism for
resistance to drought. Considering this approach, the high number of
down-regulated genes in resistance line probably is to manage energy
resources for driving growth and development under drought stress
condition. In this study, the increased gene expression in roots after
drought stress was mostly associated with increasing root length in the
resistant genotype, presumably part of an adaptive response maintained
under intense environmental pressure. This study is a starting point to
explore the network of regulatory mechanisms in Iranian wheat landraces
under drought condition.
Conclusion
The response of landraces to environmental stresses has been always an
interesting topic for research given the high adaptability of these
genotypes to the inappropriate environmental condition. The present
study provided an opportunity to investigate the differential
expression of root transcripts of bread wheat landraces under
contrasting water regimes. RNA-Seq based transcriptomic analysis was
applied to find out drought-responsive genes in resistant (L-82) and
susceptible (Marvdasht) genotypes.
Morphological and physiological assessments point out that resistant
genotype had more adaptability to drought than susceptible genotype.
The first observations of transcriptome analysis showed that a large
number of DEGs were down-regulated in the resistant genotype.
Up-regulation of genes participating in phenylpropanoids and flavonoids
pathways as well as other genes associated with oxidative stress
defense was observed in resistant landrace after drought treatment.
Up-regulation of chitinase, UDP-glucose 4-epimerase, root suberin,
thiamine metabolism and Down-regulation of genes associated with
indole-3-acetic acid in drought-treated L-82, all are related to root
growth enhancement in which more access to water is possible through
deep root penetration into the soil. Increasing the expression of genes
associated with root growth, management of energy resources through
down-regulation of a large number of DEGs, and simultaneously
scavenging of ROS generated during drought stress, altogether affirmed
that the drought tolerance in “L-82” could be due to avoidance
mechanisms. The results of RNA-Seq analysis were confirmed by RT-qPCR
assay in both genotypes under contrasting water regimes. These results
indicated that this landrace could be used as a donor of appropriate
genes to improve the bread wheat varieties for drought resistance.
Supporting information
S1 Table. List of primers for real-time PCR and their related
sequences.
(XLSX)
[181]Click here for additional data file.^ (29.3KB, xlsx)
S1 Fig. Length distribution of transcripts from wheat libraries derived
from RNA-Seq data.
(TIF)
[182]Click here for additional data file.^ (4MB, tif)
S2 Fig. Percentage of expressed transcripts in each sample.
E and M represent “L-82” and “Marvdasht” genotype, respectively. D and
N means drought stress and normal condition.
(TIF)
[183]Click here for additional data file.^ (3.9MB, tif)
S3 Fig. Enriched gene ontology (GO) categories in DEGs of ED vs. EN
comparison (Tolerant genotype).
DEGs were analyzed using AgriGo and overrepresented terms in the three
main categories “Biological Process”, “Molecular Function”, and
“Cellular Component” were filtered using Fisher’s exact test and the
Benjamini-Hochberg multiple testing correction (Q-value < 0.05).
(TIF)
[184]Click here for additional data file.^ (9.2MB, tif)
S4 Fig. Enriched gene ontology (GO) categories in DEGs of MD vs. MN
comparison (Sensetive genotype).
DEGs were analyzed using AgriGo and overrepresented terms in the three
main categories “Biological Process”, “Molecular Function”, and
“Cellular Component” were filtered using Fisher’s exact test and the
Benjamini-Hochberg multiple testing correction (Q-value < 0.05).
(TIF)
[185]Click here for additional data file.^ (6MB, tif)
S5 Fig. Enriched gene ontology (GO) categories in down-regulated DEGs
of ED vs. EN comparison (drought tolerance genotype).
DEGs were analyzed using AgriGo and overrepresented terms in the three
main categories “Biological Process”, “Molecular Function”, and
“Cellular Component” were filtered using Fisher’s exact test and the
Benjamini-Hochberg multiple testing correction (Q-value < 0.05).
(TIF)
[186]Click here for additional data file.^ (5MB, tif)
S6 Fig. Enriched gene ontology (GO) categories in up-regulated DEGs of
MD vs. MN comparison (drought sensitive genotype).
DEGs were analyzed using AgriGo and overrepresented terms in the three
main categories “Biological Process”, “Molecular Function”, and
“Cellular Component” were filtered using Fisher’s exact test and the
Benjamini-Hochberg multiple testing correction (Q-value < 0.05).
(TIF)
[187]Click here for additional data file.^ (5.5MB, tif)
S7 Fig. Mapping all non-redundant DEGs to terms in the KEGG database
and found 47 and 46 pathways were enriched in ED vs. EN and MD vs. MN,
respectively.
These different pathways were assign to 16 clades under four major KEGG
categories namely, “metabolism”, “genetic information processing”,
“cellular processes” and “environmental information processing”.
(TIF)
[188]Click here for additional data file.^ (1.5MB, tif)
S8 Fig. Clustering all unique DEGs across different comparisons were
categorized them into six clades according to their expression
patterns.
The genes in cluster I (382) down-regulated in drought tolerance
genotype while had no differential expression in control plants.
Cluster II had genes (96) that down-regulated in control plants but
have not differential expression in drought tolerant genotypes. Cluster
III (35) and IV (6) comprised down and up-regulated genes with almost
the same expression profile across the comparisons. Cluster V (37)
exhibited DEGs which are up-regulated in drought tolerance genotype
while had no differential expression in control plants. Cluster VI
shows genes (99) that are up-regulated in control plants and have not
significant expression level in drought tolerance genotype.
(TIF)
[189]Click here for additional data file.^ (55MB, tif)
S9 Fig. The KEGG pathway analysis of each cluster was performed to
assign the related biological pathways for each group of DEGs.
Genes located in cluster I-VI were assigned to 31, 5, 2, 5, 16 and 37
KEGG pathways, respectively. Cluster I and II show enriched pathways
which are down-regulated in drought tolerance genotype and control
plants. Cluster III and IV comprised enriched pathways that related
genes down and up-regulated with almost the same expression level
across the comparisons. Cluster V and VI exhibited pathways which are
up-regulated in drought tolerance genotype and control plants,
respectively.
(TIF)
[190]Click here for additional data file.^ (4.7MB, tif)
Acknowledgments