Characterization of Dynamic Regulatory Gene and Protein Networks in Wheat Roots Upon Perceiving Water Deficit Through Comparative Transcriptomics Survey

A well-developed root system benefits host plants by optimizing water absorption and nutrient uptake and thereby increases plant productivity. In this study we have characterized the root transcriptome using RNA-seq and subsequential functional analysis in a set of drought tolerant and susceptible genotypes. The goal of the study was to elucidate and characterize water deficit-responsive genes in wheat landraces that had been through long-term field and biochemical screening for drought tolerance. The results confirm genotype differences in water-deficit tolerance in line with earlier results from field trials. The transcriptomics survey highlighted a total of 14,187 differentially expressed genes (DEGs) that responded to water deficit. The characterization of these genes shows that all chromosomes contribute to water-deficit tolerance, but to different degrees, and the B genome showed higher involvement than the A and D genomes. The DEGs were mainly mapped to flavonoid, phenylpropanoid, and diterpenoid biosynthesis pathways, as well as glutathione metabolism and hormone signaling. Furthermore, extracellular region, apoplast, cell periphery, and external encapsulating structure were the main water deficit-responsive cellular components in roots. A total of 1,377 DEGs were also predicted to function as transcription factors (TFs) from different families regulating downstream cascades. TFs from the AP2/ERF-ERF, MYB-related, B3, WRKY, Tify, and NAC families were the main genotype-specific regulatory factors. To further characterize the dynamic biosynthetic pathways, protein-protein interaction (PPI) networks were constructed using significant KEGG proteins and putative TFs. In PPIs, enzymes from the CYP450, TaABA8OH2, PAL, and GST families play important roles in water-deficit tolerance in connection with MYB13-1, MADS-box, and NAC transcription factors.

Introduction

A well-developed root system is an essential component of plant productivity, as roots contribute to a host of plant growth functions (Khan et al., 2016). Roots function by anchoring plants to the soil but they also provide a platform for water and nutrient uptake and benefit plants by facilitating interactions with different symbiosis-forming micro-organisms in the rhizosphere (Comas et al., 2013; Fang et al., 2017). The spatial distribution of roots is referred to as the root system architecture (RSA) and this is known to dynamically change when encountering external elements (Lynch, 1995). An optimized RSA, both in terms of anatomy and growth attributes, can increase water absorption from different layers of soil and subsequently increase plant productivity (Wasson et al., 2012).

Despite difficulties in obtaining reliable root measurements, there have been attempts to study RSA which has resulted in a better understanding of differences in root system morphology between different species (Hochholdinger and Tuberosa, 2009; Bishopp and Lynch, 2015; Maccaferri et al., 2016). These studies suggest that plant roots generally comprise both embryo-formed roots, which include primary and seminal roots, and post embryonically-formed roots, including adventitious roots (Lynch and Brown, 2012; Khan et al., 2016). The growth features of these roots, such as elongation, longevity, and lateral dispersion are controlled by many genes and their interactions with the environment, as well as hormone homeostasis, play an important role in root development (Lynch and Brown, 2012; Smith and De Smet, 2012; Paez-Garcia et al., 2015). It has been suggested that auxin produced in both root and shoot affects the expression level of genes that are involved in cell division (Overvoorde et al., 2010). Cell-cycle genes, including Cyclin-dependent-kinases (CDKs) and D-type cyclin (CYCD), are activated by auxin and their suppression leads to reduced root branching (Nieuwland et al., 2009; Khan et al., 2016).

However, there is still a clear lack of genomic knowledge on roots compared to aboveground organs, although access to high throughput sequencing technologies, such as RNA-sequencing and previously microarrays has facilitated genomics-based comparative studies especially in model species and crop plants (Wang et al., 2009; Ma et al., 2017; Xiong et al., 2017). Wheat is one of these important crop species that has been feeding a considerable portion of humans for over 5,000 years (Peng et al., 2011). Meanwhile, around 50% of wheat arable lands are currently affected by drought and future climate change may increase it to more than 60% (Fahad et al., 2017; Hernandez-Ochoa et al., 2018; Pessarakli, 2019; Saddique et al., 2020).

Roots are integral to the detection networks for drought stress and this induce changes in hormone signaling, carbon allocation, protein, and carbohydrate metabolisms which, in turn, affects plant biomass production (Giehl et al., 2014; Paez-Garcia et al., 2015). A better understanding of the genomic regions controlling root development would therefore assist in identifying key biosynthetic pathways involved in root growth. In this paper we present the results of (I) biochemical and root morphology studies in wheat genotypes with contrasting drought response phenotypes, (II) perform comparative transcriptomics of root cells, (III) assess functional enrichment of water deficit-responsive biosynthetic pathways and regulatory elements, (IV) and create a protein-protein interaction network to elucidate underlying genomic loci and dynamic elements in root development under water deficit conditions.

Materials and Methods

Plant Materials and Water Treatments

A collection of 320 Iranian bread wheat accessions was evaluated during the 2016–17 and 2017–18 cropping seasons for a number of agronomic traits at the Agricultural Research Lands of the Department of Agronomy and Plant Breeding, University of Tehran (Rahimi et al., 2019). From this germplasm, 100 landraces were selected based on drought tolerance indices and their performance under drought conditions, ranging from susceptible to tolerant (Supplementary Table 1). By applying polyethylene glycol 6000 (PEG 6000) on seeds, we simulated water deficit conditions at the seedling growth stage. The growth parameters of seedlings assessed included fresh and dry matter, length of roots and shoot, percentage of germination and germination speed. The genotypes were also grown in the greenhouse and exposed to stress conditions [50% field capacity (FC)] 4 weeks after planting. Leaves samples were collected 1 week after exposure to stress conditions to study antioxidant capacity.

Assessment of Hydrogen Peroxide (H₂O₂) and Malondialdehyde (MDA) Content

To measure H₂O₂ content, 0.40 g of fresh leaves were homogenized using 6 ml of 0.1% (w/v) trichloroacetic acid (TCA). Samples were centrifuged at 12,000 × g and 0.6 ml of the resulting supernatant was transferred to a 2 ml tube mixed with 0.5 ml of 10 mM phosphate buffer and 0.9 ml 1 M KI. The solutions were kept under dark condition for 1 h after which the absorbance of the solution at 390 nm was scored (Loreto and Velikova, 2001). To measure MDA, 0.2 g of fresh leaves were ground by adding 5 ml 0.1% TCA. After centrifuging the homogenate at 8,000 × g for 5 min, 1 ml of the resulting supernatant was pipetted into a new tube and mixed with 4.5 ml of 20% TCA. The solution was centrifuged at 4,000 × g for 10 min and the absorbance of the solution was read at 532 nm (Heath and Packer, 1968).

Measurement of Proline Content and Superoxide Dismutase (SOD) Activity

To measure proline content we used a ninhydrin reaction buffer based on the method of Bates et al. (1973). The proteins of leaves were extracted using the Bradford (1976) method with a phosphate buffer (pH 7.0). The resulting supernatant was further used to determine soluble proteins and antioxidant enzyme activity. The activity of SOD was measured using the procedure from Dhindsa et al. (1982), where the enzyme activity is measured by its ability to inhibit photochemical reduction. An amount of SOD which inhibits the photochemical reduction of nitroblue tetrazolium by 50% is referred to one unit of the enzyme.

Preparing RNA Samples

After field trials and in-situ studies, 19 landraces were selected, ranging from extremely susceptible to very drought tolerant. These landraces were then used for further investigations at the Swedish University of Agricultural Sciences, Uppsala, Sweden. Seeds of these genotypes were kindly provided by the International Maize and Wheat Improvement Center (CIMMYT). Isolated landraces were grown in the greenhouse by applying 16/8 h light and dark and were then exposed to either normal irrigation or water deficit conditions 4 weeks after planting. Roots of drought stressed-plants and their corresponding controls from the normal irrigation scheme were sampled at three time-points; 9, 18, and 25 days after applying stress (DA). Water stress treatment involved supplying no water to pots in the stress treatment while control pots were irrigated regularly every 4 days until sampling. At each time point, three biological replicates were collected for each treatment, resulting in a total of 342 samples that were evaluated for root attributes. Root length, root and shoot biomass and root/shoot ratio of all samples were recorded, and all samples were subsequently frozen in liquid nitrogen and stored at −70°C. After assessing grain yield under drought conditions, seedling traits, antioxidant activity, osmoprotectant capacity, and root features, three landraces, including two tolerant and one susceptible, were selected for RNA extraction. We tried to select genetically differentiated tolerant landraces to get insights into possible differences in biosynthesis pathways used under severe conditions. Overall, we used three landraces, PI627299 (tolerant), PI627038 (tolerant), and PI624837 (susceptible), two water treatments (normal irrigation and water deficit stress) and two time points (nine and 18 DA) with three biological replicates for a total of 36 unique sample combinations that were used for RNA extraction.

RNA Isolation

Total RNA was extracted using a QIAGEN RNeasy^®Plant Mini Kit. Briefly, 100 mg of plant material was disrupted in liquid nitrogen and transferred to an RNase-free 2 ml microcentrifuge tube. 450 μl of RLT Buffer was added and the sample was vortexed vigorously. The lysate was transferred to a QIAshredder spin column and centrifuged for 2 min at 13,000 × g. The supernatant was transferred to new microcentrifuge tube and 0.5 volume of ethanol (96%) was added. 650 μl of the resulting sample was transferred to an RNeasy Mini spin column placed in a 2 ml collection tube and centrifuged at 8,000 × g for 15 s. The flow-through was discarded and 700 μl Buffer RW1 was added to the RNeasy Mini spin and centrifuged at 8,000 × g for 15 s, then 500 μl Buffer RPE was added and centrifuged at 8,000 × g for 15. This step repeated two times. Finally, total RNA was collected in 1.5 ml collection tube using 40 μl of RNase-free water and stored at −70°C before sequencing. The quantity and quality of all samples were measured on an Agilent 2100 Bioanalyzer.

Library Construction, Transcriptome Sequencing, and Reads Alignment

The sequencing libraries were prepared at the SNP&SEQ Technology Platform at Uppsala University, Sweden. Libraries were prepared from 1 μg total RNA using the TruSeq stranded mRNA library preparation kit (Cat# 20020594/5, Illumina Inc.) including polyA selection. The library preparation was performed according to the manufacturers' protocol (#1000000040498). The libraries were sequenced on a NovaSeq SP flowcell using paired end reads and a read length of 50 bp. Raw reads were filtered for adapter sequences, low-quality sequences and ambiguous nucleotides. All downstream analyses were performed using clean data. Before aligning reads, the wheat reference genome, released by the International Wheat Genome Sequencing Consortium (IWGSC, available at (ftp://ftp.ensemblgenomes.org/pub/plants/release-43/fasta/triticum_aestivum/dna/), was downloaded and indexed using HISAT2 (https://ccb.jhu.edu/software/hisat2/). All reads then were aligned using HISAT2 and the featurCounts software (http://bioinf.wehi.edu.au/featureCounts/) was used to count the number of reads which mapped to each gene using the aligned reads and an annotation file for the wheat genome (https://plants.ensembl.org/Triticum_aestivum/Info/Index). Primary controls were performed on raw counts by filtering genes not detected across all samples. After filtering, different normalization approaches were assessed and sample qualities were evaluated by plotting the distribution of counts based on their mean variances. Correlations between samples were visualized using a variance-stabilizing-transformation (VST) to quantify the correlation between biological replicates. Moreover, the variability among samples was assessed by performing a principal component analysis (PCA) on normalized read counts.

Differential Expression Analysis and Visualization

To assess differential expression, we used the DESeq2 package (Love et al., 2014) implemented in Rstudio (Team, 2020). For this purpose, filtered raw counts were used in a GLM model where genotype, condition and time-point were considered as covariates in the model. DESeq2 internally corrects counts for sequencing depth and RNA compositional biases using the Median of ratios method (Team, 2020). After estimating the size factor, gene-wide dispersion was performed to measure variation in the dataset. A Negative Binomial model was fitted for each gene and tests for differential expression were performed based on the values of the log-fold change using a Wald test under the null hypothesis that groups are not different and logFC = 0. Four main comparisons were made between genotypes, conditions, and time-points to identify differentially expressed genes (DEGs) at p-value < 0.05. To display mean expression vs log fold change for all genes in each comparison we visualized MA plots through the DESeq2 R-package. Furthermore, the EnhancedVolcano R-package (Blighe et al., 2019) was used to plot log fold change vs adjusted p-values. A heatmap was used to visualize clusters of DEGs based on their relative expression levels. The location of highly significant genes were plotted at their relevant chromosomes using the shinyCircos online tools (http://shinycircos.ncpgr.cn/).

GO Enrichment, Functional Analysis, and Predication of Regulatory Elements

Gene ontology enrichment analyses were first carried out using the g-profiler database (Raudvere et al., 2019) using GO terms from wheat to determine the molecular function (MF), biological process (BP), cellular components (CC), and potential pathways. All tests were performed at a Bonferroni corrected p-value of 0.05. Based on the list of DEGs from each comparison, sequences of homologs rice genes were obtained using BioMart (http://plants.ensembl.org/biomart/martview/8a426f1eb50080f5f0e2770c9d719a44). These sequences were also used for GO enrichment analyses using the KOBAS software (Xie et al., 2011) to test for statistically enriched pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/) database. A list of wheat gene sequences for DEGs was prepared for each comparison using the BioMart database and used as input on iTAK (Zheng et al., 2016) and PlantTFcat (Dai et al., 2013), two tools used to predict transcription factors (TFs), transcriptional regulators (TRs), and protein kinases (PKs).

Protein-Protein Interaction (PPI) Network

To construct the PPI network we used the Search Tool for the Retrieval of Interacting Genes/Proteins database or STRING v11 (Szklarczyk et al., 2021). Before that, we extracted peptide sequences of DEGs enriched in significant biosynthetic pathways and peptide sequences of candidate transcription factors from BioMart and used a list of these sequences as the input for STRING to find putative interactions. These networks include known interactions (from curated databases and experimentally determined), predicted interactions (gene neighborhood, gene fusions, and gene co-occurrence), and other interactions (textmining, co-expression and protein homology) retrieved at high confidence interaction score of 0.7.

Modeling 3D Protein Structure

After detecting key proteins involved in PPI networks, we used SWISS-MODEL homology modeling to predict 3D structures and quality assessments (Waterhouse et al., 2018). The SWISS-MODEL pipeline relies on ProMod3 (Studer et al., 2021) and comprises the following steps: identification of structural template(s), alignment of target sequence and template structure(s), model-building, and model quality evaluation. Through an automatic pipeline, suitable templates identifies based on BLAST (Camacho et al., 2009), and HHblits (Steinegger et al., 2019), so that the template's quality has been predicted from features of the target-template alignment. A set of 50 top-ranked templates is then selected from multiple templates according to a score obtained by combining sequence coverage and sequence similarity. Candidate models are then built based on the target-template alignment using ProMod3. To select the final models, we considered GMQE (Global Model Quality Estimation) and QMEAN (Qualitative Model Energy Analysis) scores. GMQE is a quality estimation that combines properties from the target–template alignment and the template structure. The GMQE score is expressed as a number between 0 and 1, reflecting the expected accuracy of a model built with that alignment and template, normalized by the coverage of the target sequence. Higher numbers indicate higher model reliability. On the other hand, the global and per-residue model quality were assessed through QMEAN scoring function (Studer et al., 2020).

Validation of RNA-Seq Analysis Using RT-qPCR

To validate the result from the RNA-seq analyses, a set of nine genes were randomly selected among DEGs (Supplementary Table 2) and tested using RT-qPCR. In brief, 1.5 μg of total RNA was treated with dsDNase provided with the Thermo Scientific™ Maxima™ First Strand cDNA Synthesis Kit for RT-qPCR and used for cDNA synthesis using the following thermocycler program: incubation for 10 min at 25°C followed by 15 min at 50°C. The reaction was terminated by heating at 85°C for 5 min. The RT-qPCR reactions were performed using Thermo Scientific Luminaris HiGreen High ROX qPCR Master Mix. The reactions solution was prepared in a volume of 10 μl by using 5 μl Master Mix (2X), 0.3 μl forward and reverse primers, 1 μg cDNA and nuclease-free water up to 10 μl. A three-step cycling protocol was used with the following thermocycler program: UDG pre-treatment at 50°C for 2 min, initial denaturation at 95°C for 10 min, 40 cycles of denaturation at 92°C for 15 s, annealing at 60°C for 30 s and extension at 72°C for 30 s. To control the specificity of all products a melting curve analysis was used by increasing the temperature from 60 to 95°C and read every 0.5°C. To normalize Ct values of each reaction, an internal reference gene (actin) was used in the relative expression equation.

Result

Early Growth Under Simulated Water Deficit and Root Attributes at Tillering

Water deficit reduced all germination-related traits at−10 bar, although this decrease was not significant in tolerant genotypes (Table 1). Root length under stress conditions for the two tolerant genotypes, PI627038 and PI627299, decreased ~15% relative to control conditions whereas the susceptible genotype, PI624837, suffered a reduction of ~60%. The total fresh and dry weight of seedlings was significantly reduced in all genotypes under stress conditions compared to the control. The highest total germination percentage (TGP) and normal germination percentage (NGP) was observed for PI627299 under control condition (100%). However, under simulated water deficit through PEG application, these values dropped by 90%. The susceptible genotype, PI624837, had the lowest germination percentage, both for total and for normal germination. To test the germination ability of the tolerant genotypes, a severe water stress of−12 bar was also used and the results are shown in Figure 1A. Under these conditions, the susceptible genotype PI624837, showed no germination at all. The speed of germination was calculated based on the coefficient of velocity of germination (CVG) and tolerant genotypes showed relatively high values even under a water stress of−10 bar (Figure 1B).

TABLE 1

Table 1. Assessment of tolerant and susceptible genotypes according to germination parameters under simulated water deficit using PEG 6000.

FIGURE 1

Figure 1. Germination study of selected wheat genotypes. (A) Response to simulated water deficit at−12 bar, the tolerant genotype PI627038 is showing root development including radicle, seminal roots; and shoot development including coleoptile emergence, while no germination observed in susceptible genotype PI624837. (B) Coefficient velocity of germination (CVG) at−10 bar. (C) Root length. (D) Root weight. (E) Shoot weight. (F) Shoot/root ratio. T, drought-tolerant; S, drought-susceptible; DA, days after applying water-deficit stress. Treatments are grouped based on Duncan multiple range test.

On the other hand, after a long-term application of water deficit at tillering stage, we measured root features at three time points. All four studied traits were significantly influenced by both stress treatment and sampling time. However, the studied genotypes showed different responses to control and water deficit conditions. The two tolerant genotypes had greater root length (Figure 1C), root biomass (Figure 1D), shoot weight (Figure 1E), and root/shoot ratio (Figure 1F) compared to susceptible genotype, especially under stress treatments. PI627038 (T) had the greatest root length, root and shoot biomass under normal and stress conditions at the 25 DA time point. Water deficit induced root growth and increased the root/shoot ratio for all plants although the increase for the two tolerant genotypes were greater than for the susceptible genotype.

Biochemical Assessment of Selected Landraces

We assessed the antioxidant capacity of the three selected genotypes after applying water deficit (Supplementary Figure 1). We observed significant differences in proline, MDA, and H₂O₂ content and in SOD activity between the tolerant and susceptible genotypes. PI627038 and PI627299, showed the highest proline content and SOD activity under water deficit treatments (18.75 mg g⁻¹ FW and 0.0683 μmol of guaiacol oxidized min⁻¹ mg⁻¹ protein, respectively). Water deficit increased H₂O₂ in all genotypes and the level of cell membrane damage increased significantly. However, damage to phospholipid membranes was greater in PI624837 (S) (46.35 μmol g⁻¹ FW, MDA) compared to the tolerant genotypes.

RNA-Seq Experiment, Gene Dispersion, and DEGs

A total of 1.06 billion clean sequencing reads were generated from the 36 cDNA libraries. The mapping rate of reads ranged from 78.1 to 84.5% across libraries without considering multi-mapped reads (average = 81.7%). The filtered and normalized data from VST were used for data exploration. The three biological replicates of each sample were highly correlated (r > 0.95) a consistently clustered into the same groups (Supplementary Figure 2). A principal component analysis was used to determine how variation in gene expression resulted from variation due to genotypes, conditions, and sampling time-points (Figure 2A). The PCA plots indicates that the three selected genotypes have highly differentiated expression profiles, confirming the results from field trials and other biochemical evaluations.

FIGURE 2

Figure 2. The exploited result of VST normalization method and Venn diagrams. (A) PCA plot of samples based on three covariates including genotypes, environments and time-points that showing differentiation between tolerant (green and blue spots) and susceptible genotypes (red spots). Venn diagram of up (B) and down-regulated (C) genes in four main comparisons. G1, PI627299 (T); G2, PI627038 (T); G3, PI624837 (S); Normal, well-watering (control); Stress, water-deficit stress; 9DA, 9 days after applying water-deficit stress; 18DA, 18 days after applying water-deficit stress; T, drought-tolerant; S, drought-susceptible.

The gene dispersion procedure implemented in the DESeq2 package was used to estimate the variation in the dataset. This method is more reliable than using the variance and/or standard deviation of gene expression, as these measures are dependent on the mean expression level (Supplementary Figure 3). A plot of the relationship between mean counts and variance suggest that the variables are positively correlated, meaning that highly expressed genes also have higher expression variance. One way to control for this mean-variance association is to calculate the coefficient of variation which standardizes the variation by mean expression level. However, the assumption in DESeq2 is that genes with similar expression levels have similar dispersion values. Dispersion estimates are computed gene for gene using a maximum likelihood method (the internal version of dispersion in DESeq2) and gene-wise estimates are then shrunk to the fitted curve of the previously estimated dispersion. This method helps in reducing false positives in the DGEs analyses.

Putative DEGs were identified by considering a false discovery rate (FDR) of 0.05 and were further filtered to identify genes showing a FDR < 0.01 and having an absolute fold change > 2 (Figures 2B,C). In total, 14,187 unique genes were differentially expressed across the four main comparisons, with 4,738 DEGs related to the contrast between control and stress condition (1,815 up-regulated and 2,923 down-regulated), 4,008 DEGs were found when contrasting individuals PI627299 vs. PI624837 (1,999 up-regulated and 2,009 down-regulated), 3,613 DEGs for contrasting PI627038 vs. PI624837 (1,833 up-regulated and 1,780 down-regulated), and 1,828 genes when contrasting time-points 9 DA vs. 18 DA (1,210 up-regulated and 618 down-regulated). Overall, different comparisons showed some degree of overlap in both up and down-regulated genes, with the highest number of up and down-regulated overlapping genes observed between tolerant genotypes (1,053 and 584, respectively). A total of 23 up-regulated and 87 down-regulated genes were consistently expressed in all comparisons.

The results of the DEG analyses were further visualized using volcano plots, in which genes are sorted by log fold changes and their adjusted p-values (Figure 3). In all comparisons we observed plenty of highly significant genes. For example, in the comparison between tolerant and susceptible genotypes, TraesCS1B02G459100, TraesCS4D02G360300, TraesCS4D02G357100, and TraesCS1A02G335200 were highly up-regulated while TraesCS2A02G216200, TraesCS2B02G63800, TraesCS5A02G482400, TraesCS2B02G484200, and TraesCS2B02G614100 were down-regulated (Figures 3A,B). When comparing the stress versus control condition and time-points the upregulated genes were TraesCS6B02G383500, TraesCSU02G049400, TraesCS4B02G312200, and TraesCS3B02G256400 while TraesCS1B02G290000, TraesCS7D02G201400, TraesCS5D02G046500, TraesCS1D02G280200, TraesCS1B02G290000, and TraesCS1A02G281000 were among the strongly down-regulated genes (Figures 3C,D). Visualizing the genomic location of the these genes confirmed that all chromosomes harbor genes contributing to the water deficit response (Figure 3E). A total of 1,604 genes were located on chromosomes from the B genome and 932 genes were located on chromosomes from the D genome, representing the largest and smallest number of genes when comparing genotypes PI627299 (T) and PI624837 (S). Similarly, when comparing genotypes PI627038 (T) and PI624837 (S), the B genome had the greatest proportion of genes responding to water deficit. Moreover, in the comparison between water deficit and control conditions, the D genome contained 1,611 significant genes, while the number of genes in A and B genomes were roughly equal (1,512 and 1,525 for A and B genomes, respectively). In the comparison between time-points, 639 genes responding to water deficit were located on the D genome. No drought-responsive genes located in the centromeric regions of the different chromosomes were observed. Important up and down regulated genes are distinguished in blue and red colors, respectively (Figure 3F).

FIGURE 3

Figure 3. Volcano plot and chromosome wise distribution of highly significant up and down-regulated genes. (A) PI627299 (T) vs. PI624837 (S). (B) PI627038 (T) vs. PI624837 (S). (C) Water deficit vs. control. (D) Time-point 18 DA vs. 9 DA. (E) Circos plot of all significant genes, in comparison 1/ PI627299 (T) vs. PI624837 (S) (purple spots), in comparison 2/ PI627038 (T) vs. PI624837 (S) (light blue spots), in comparison 3/ condition water deficit vs. control (green spots), and in comparison 4/ time-point 18 DA vs. 9 DA (red spots). (F) Chromosomal region of important key genes and their encoding proteins are shown inside the plot.

A heat map of top 40 most significant genes in each comparison were plotted across all samples, showing both up and down-regulated genes (Supplementary Figure 4). Interestingly, TraesCS1A02G335200, TraesCS2B02G027300, TraesCS3B02G037400, TraesCS5B02G302500, and TraesCSB02G450800 showed different pattern of expression between the two tolerant genotypes for comparisons 1 and 2. A large fraction of the 40 selected genes showed differential expression patterns between tolerant and susceptible genotypes. With an increasing number of days after applying water stress, most genes tended to have reduced transcript levels in the different samples. The copy number of key genes are also provided in Supplementary Tables 3, 4.

Transcription Factors Involved in Water-Deficit Tolerance

A total of 1,377 DGEs were predicted to belong to different families of transcription factors (TFs, Table 2). A set of 443 genes encoding TFs, that showed different expression patterns between normal and stress treatments, were identified. The number of TFs identified in the comparisons between tolerant and susceptible landraces were approximately the same (328 and 356 TFs, respectively). Moreover, 250 TFs were identified when the period of water deficit was increased from 9 to 18 days. Most of these TFs belong to the AP2/ERF-ERF, MYB-related, bHLH, B3, and MADS-M-type classes as genotype specific TFs. However, MYB, NAC, WRKY, C2H2, and bZIP were also major classes of TFs specific to short/long-term water deficit. Moreover, transcriptional regulators were estimated, and these largely consisted of mTERF, SNF2, TRAF, and TAZ (Supplementary Table 5).

TABLE 2

Table 2. A summary of identified TFs among the different sets of DEGs in wheat.

Functional Analysis Using g-Profiler

A functional analysis of DEGs with g-Profiler showed that a total of 45, 32, and 8, GO terms were involved in molecular function (MF), biological process (BP), cellular components (CC), respectively and 3 KEGG pathways were significantly enriched in the comparison between PI627299 (T) vs. PI624837 (S) (Supplementary Figure 5A). Comparing PI627038 (T) vs. PI624837 (S) revealed that 36, 27, and 6 GO terms were associated with MF, BP, and CC, respectively (Supplementary Figure 5B). In comparison between water deficit and normal conditions 92, 85, and 18 GO terms were involved in MF, BP, CC, and 5 KEGG pathways were significantly enriched (Supplementary Figure 6A). Finally, comparing time points 18 vs. 9 DA showed that 38, 32, and 9 GO terms were involved in MF, BP, CC, and 7 KEGG pathways were enriched (Supplementary Figure 6B). These findings show that extracellular region, apoplast, cell periphery, external encapsulating structure, cell wall, plasma membrane were the main cellular components affected by water deficit conditions. Furthermore, genes in box C/D snoRNP complex, ASTRA complex, intrinsic component of membrane, integral component of membrane, anchored component of membrane, extracellular space, photosystem showed significant changes in the previously mentioned cellular components.

Functional Annotation of DEGs and KEGG Pathways on KOBAS

The KEGG pathways and GO terms related to DEGs were obtained from the gene enrichment analyses using KOBAS and are listed in Supplementary Tables 6–8. Five pathways, glutathione metabolism, flavonoid biosynthesis, phenylpropanoid biosynthesis, diterpenoid biosynthesis, valine, leucine, and isoleucine degradation were significantly enriched when comparing tolerant and susceptible genotypes (Supplementary Table 6). Six pathways, phenylpropanoid biosynthesis, glutathione metabolism, flavonoid biosynthesis, plant hormone signal transduction, starch and sucrose metabolism and carotenoid biosynthesis pathways were significantly enriched when comparing water deficit with normal conditions (Supplementary Table 7). Finally, five pathways, glutathione metabolism, flavonoid biosynthesis, phenylpropanoid biosynthesis, zeatin biosynthesis, and photosynthesis process, were significantly enriched between the two time points (Supplementary Table 8). These results suggest that phenylpropanoid biosynthesis, glutathione metabolism, flavonoid biosynthesis, photosynthesis process and valine, leucine and isoleucine degradation were strongly affected pathways and mechanisms in this study. In the flavonoid biosynthesis pathway, trans-cinnamate 4-monooxygenase, flavonoid 3'-monooxygenase, anthocyanidin synthase, anthocyanidin reductase were the main enzymes showing different activities between one of the tolerant landraces and the susceptible genotype (Figure 4). In the valine, leucine and isoleucine degradation process, a branched-chain amino acid aminotransferase, 2-oxoisovalerate dehydrogenase E1 component beta subunit, isovaleryl-CoA dehydrogenase, 3-hydroxyisobutyryl-CoA hydrolase, acetyl-CoA C-acetyltransferase and alanine-glyoxylate transaminase/(R)-3-amino-2-methylpropionate-pyruvate transaminase were the enzymes that showed the largest difference between the tolerant and susceptible genotype (Figure 5). In the phenylpropanoid biosynthesis pathway, the activity of phenylalanine ammonia-lyase, trans-cinnamate 4-monooxygenase, 4-coumarate–CoA ligase, ferulate-5-hydroxylase, cinnamyl-alcohol dehydrogenase, beta-glucosidase and peroxidase were considerably altered when plants were exposed to water deficit (Figure 6). Interestingly, some enriched genes identified through the functional analysis were involved in photosynthesis (Figure 7). Photosystem II 10kDa protein, photosystem II oxygen-evolving enhancer protein 1, photosystem II oxygen-evolving enhancer protein 2, photosystem II oxygen-evolving enhancer protein 3, photosystem II PsbY protein, photosystem II Psb27 protein, photosystem I subunit X and photosystem I subunit PsaO were affected after long-term water deficit. Although there are many other genes and enzymes identified for the different pathways or mechanisms, not all of them have so far been characterized, neither in rice nor in wheat. The rest of significant pathways in the different comparisons are presented in Supplementary Figures 7–12.

FIGURE 4

Figure 4. Flavonoid biosynthesis pathway. (1.14.14.91) Cytochrome P450 CYP73A100. (1.14.14.82) Flavonoid 3'-monooxygenase. (1.14.20.4) Leucoanthocyanidin dioxygenase. (1.3.1.77) Anthocyanidin reductase. Red color represents significant and characterized enzymes based on Oryza sativa annotation. Green color represents significant and non-characterized enzymes. White color represents characterized enzymes based on Arabidopsis thaliana annotation.

FIGURE 5

Figure 5. Valine, leucine and isoleucine degradation pathway. (2.6.1.42) Branched-chain-amino-acid aminotransferase 5, chloroplastic isoform X1. (1.2.4.4) 2-oxoisovalerate dehydrogenase subunit beta 1, mitochondrial. (1.3.8.4) Isovaleryl-CoA dehydrogenase mitochondrial. (2.3.1.9) Acetyl-CoA acetyltransferase, cytosolic 1. (2.6.1.40) Alanine–glyoxylate aminotransferase 2 homolog 1, mitochondrial. Red color represents significant and characterized enzymes based on Oryza sativa annotation. Green color represents significant and non-characterized enzymes. White color represents characterized enzymes based on Arabidopsis thaliana annotation.

FIGURE 6

Figure 6. Phenylpropanoid biosynthesis pathway. (4.3.1.24) Phenylalanine ammonia-lyase. (1.14.14.91) Cytochrome P450 CYP73A100. (F5H) Cytochrome P450 84A1. (3.2.1.21) Beta-glucosidase 1. (1.1.1.195) Cinnamyl alcohol dehydrogenase 2. (1.11.1.7) Peroxidase 70-like. Red color represents significant and characterized enzymes based on Oryza sativa annotation. Green color represents significant and non-characterized enzymes. White color represents characterized enzymes based on Arabidopsis thaliana annotation.

FIGURE 7

Figure 7. Water deficit-affected subunits of photosystems in photosynthesis chain. (PsbO) Oxygen-evolving enhancer protein 1. (PsbP) Oxygen-evolving enhancer protein 2. (PsbQ) Photosynthetic NDH subunit of lumenal location 2. (PsbR) Photosystem II 10 kDa polypeptide. (PsbY) Photosystem II core complex proteins psbY. (Psb27) Photosystem II repair protein PSB27-H1. (PsaK) Photosystem I reaction center subunit psaK. (PsaO) Photosystem I subunit O. Red color represents significant and characterized enzymes based on Oryza sativa annotation. Green color represents significant and non-characterized enzymes. White color represents characterized enzymes based on Arabidopsis thaliana annotation.

PPI Network and 3D Protein Structure

To study the interaction between proteins from transcription factors families and proteins enriched in highly significant KEGG pathways, three PPIs were constructed using the STRING database. We used peptide sequences of TFs and KEGG genes as input for STRING retrieval. Genotype specific and drought-responsive interaction networks predicted 329, 212, and 259 node proteins in an enrichment p-value of <1.0e−16 with a high confidence interaction score of 0.7 (Table 3). This shows all three interaction networks are significantly enriched and that the proteins are expected to be biologically connected, confirming the result of the functional analysis. Networks were clustered using the specified Markov Cluster Algorithm (MCL), which is an approach for both weighted and unweighted networks. Detailed information about these three PPIs is presented in Supplementary Tables 9–11.

TABLE 3

Table 3. PPI network statistics specified to tolerant genotypes and water deficit conditions.

Genotype-Specific Protein Networks

Nineteen connected and disconnected clusters were obtained for PI627299 (T) vs PI624837 (S), although five clusters showed stronger interactions (Figure 8). The first cluster is showing enzymes responsible for the phenolic component's pathways (salmon-colored), and the key elements are two uncharacterized proteins belongs to the cytochrome P450 (CYP450) family located on chromosomes 3A and 5A (TraesCS3A02G136100 and TraesCS5A02G534000, respectively). These key proteins are interacting with phenylalanine ammonia-lyase (PAL, TraesCS2D02G377500) and 23 other characterized or uncharacterized proteins. The second important interaction is colored in red, and describes interactions involving glutathione S-transferase 1 (GST1, TraesCS4B02G059300, blue colored), and indicates important enzymes in glutathione metabolism. Furthermore, interactions of important transcription factors, such as MADS-box, BZIP, NAC and HsfC1e (Heat shock factor C1e,) with other proteins were also identified in several clusters.

FIGURE 8

Figure 8. Interaction networks of related DEGs in significant KEGG pathways and putative TFs identified by comparing PI627299 (T) vs. PI624837 (S). e.g., TraesCS2D02G377500 (Phenylalanine ammonia-lyase), TraesCS5B02G11510 (MADS-box transcription factor TaAGL31), TraesCS5D02G401500 (VRN-D1), TraesCS3B02G255600 (Heat shock factor C1e), TraesCS4B02G059300 (Glutathione S-transferase 1).

The second enriched PPI network (PI624837 (S) vs. PI627038 (T)) showed similar interactions, while highlighting the interaction of two additional important TFs, MYB13-1 and GATA, with other proteins (Supplementary Figure 13). The key protein in this cluster was an uncharacterized protein (TraesCSU02G031300) containing a plant homeodomain (PHD). PHD is a C4HC3 zinc-finger-like motif that has been found in some nuclear proteins that are potentially involved in epigenetic and chromatin-mediated transcriptional regulation. This key protein has also a strong interaction with the MYB13-1 transcription factor, which is located at the primary shell of the network. Another important cluster colored in brown, and three key proteins including TaANS-A2 (TraesCS6A02G041800), TaANS-D1 (TraesCS6D02G004300), and an uncharacterized protein (TraesCS2B02G613000) belonging to the iron/ascorbate-dependent oxidoreductase family interacted with other 19 proteins. Both TaANS-A2 and TaANS-D1 are anthocyanidin synthase from the iron/ascorbate-dependent oxidoreductase family.

General Protein Networks Involved in Water-Deficit Tolerance

In the last PPI, the connection between encoded proteins from significant KEGG pathways and putative TFs obtained from water deficit vs control comparison was generated at a high confidence score and grouped into 11 important clusters (Figure 9). The first cluster contained 32 biologically connected proteins mostly with phenylalanine ammonia-lyase activity (red-colored) and cytochrome P450 family (e.g., TraesCS3A02G083600, TraesCS1D02G039300, TraesCS5B02G468300, TraesCS1D02G019400, etc.). Two uncharacterized proteins (TraesCS2B02G291100 and TraesCS6A02G266700) with AMP-binding domain were key elements and interacting with other proteins in this cluster. There was a potential connection between the first cluster and the second cluster through TaANS-A2 (TraesCS6A02G041800), TaANS-D1 (TraesCS6A02G001500), and another uncharacterized protein (TraesCS2B02G613000). In the third important cluster, the interaction between enzymes with sucrose synthase, alpha-amylase and beta-amylase activities were generated, where sucrose synthase or sucrose-cleaving enzyme provides UDP-glucose and fructose for various metabolic pathways was the key factor. However, there were also important regulatory elements in each cluster such as BZIP transcription factor (cluster 4, green-yellow), putative serine/threonine protein kinase (cluster 5, lime green), abscisic acid-inducible protein kinase (cluster 5, lime green), SNF1-type serine/threonine protein kinase (cluster 5, lime green) and ABA 8'-hydroxylase (cluster 11, pink).

FIGURE 9

Figure 9. Interaction networks of related DEGs in significant KEGG pathways and putative TFs identified by comparing water deficit vs. control. e.g., TraesCS2D02G377500 (Phenylalanine ammonia-lyase), TraesCS6A02G041800 (Anthocyanidin synthase/TaANS-A2), TraesCS6A02G001500 (Anthocyanidin synthase/TaANS-D1), TraesCS7A02G539000 (Putative serine/threonine protein kinase), TraesCS2B02G174400 (1-Cys peroxiredoxin PER1), TraesCS5D02G244900 (TaABA8OH2 or ABA 8'-hydroxylase).

Prediction of the 3D Structure of Key Proteins

The SWISS-MODEL template library was searched with BLAST and HHblits for evolutionary related structures matching the target sequences of key proteins involved in PPIs. For genotype-specific networks, the 3D structures of available core proteins (according to UniProt ID) were predicted as follows (Figure 10): W5D4C8/monooxygenase (TraesCS3B02G154000), W5AA79/R1 protein (TraesCS2A02G312700), Q0WYI7/anthocyanidin synthase (TraesCS6A02G041800), and S4V9N0/MYB13 (TraesCS3A02G535100). The SWISS-MODEL template library was also searched for related structures matching the target sequences of key proteins in general drought-responsive networks (Figure 11) as follows: Q43210/phenylalanine ammonia-lyase (TraesCS2A02G381100), Q6W8Q2/1-Cys peroxiredoxin PER1 (TraesCS2B02G174400), Q02066/abscisic acid-inducible protein kinase (TraesCS2A02G493800), and T2HNE3/ABA 8'-hydroxylase or (+)-abscisic acid 8'-hydroxylase (TraesCS5B02G236500). Generally, higher GMQE scores observed for the predicted models in this network compared to genotype-specific networks. Templates with a similarity of more than 60% were selected to build the models. According to GMQE and QMEAN scores, predicted models were evaluated and the most reliable 3D structures with higher GMQE and QMEAN Z score close to zero were extracted for further characterization. The “Local Quality” plot shows, for each residue of the model (x-axis), the expected similarity to the native structure (y-axis). In the “Comparison” plot, model quality scores of individual models are related to scores obtained for experimental structures of similar size. The x-axis shows protein length (number of residues). The y-axis is the normalized QMEAN score. Every dot represents one experimental protein structure. Black dots are experimental structures with a normalized QMEAN score within 1 standard deviation of the mean (|Z-score| between 0 and 1), experimental structures with a |Z-score| between 1 and 2 are gray. The experimental structure that is even further from the mean is light gray. The actual model is represented as a red star.

FIGURE 10

Figure 10. Swiss model 3D structure prediction and model accuracy based on QMEAN for following proteins (according to UniProt ID): (A) W5D4C8 (monooxygenase), (B) W5AA79 (R1 protein), (C) Q0WYI7 (anthocyanidin synthase), and (D) S4V9N0 (MYB13). The “Local Quality” plot shows, for each residue of the model (reported on the x-axis), the expected similarity to the native structure (y-axis). In the “Comparison” plot, model quality scores of individual models are related to scores obtained for experimental structures of similar size.

FIGURE 11

Figure 11. Swiss model 3D structure prediction and model accuracy based on QMEAN for following proteins (according to UniProt ID): (A) Q43210 (Phenylalanine ammonia-lyase), (B) Q6W8Q2 (1-Cys peroxiredoxin PER1), (C) Q02066 (Abscisic acid-inducible protein kinase), and (D) T2HNE3 (ABA 8'-hydroxylase). The “Local Quality” plot shows, for each residue of the model (reported on the x-axis), the expected similarity to the native structure (y-axis). In the “Comparison” plot, model quality scores of individual models are related to scores obtained for experimental structures of similar size.

Experimental Validation of DEGs via qRT-PCR

The results of the RNA-seq analyses were validated by qRT-PCR using a total of nine randomly selected genes (five up-regulated and four down-regulated) that were normalized to the expression of the housekeeping gene actin (Supplementary Figure 14). The up-regulated genes were involved in rRNA N-glycosylase activity, catalytic activity, peroxidase activity, transmembrane transporter activity, ADP binding, and the down-regulated genes were involved in carbohydrate binding, nucleic acid binding, inorganic phosphate transmembrane transporter activity, methyltransferase activity, respectively. The expression pattern (up and down-regulations) observed in the qRT-PCR experiment was highly similar to the observed RNA-seq log2 fold change for all selected genes (RSquare and adjusted RSquare were 0.90 and 0.89, respectively). This shows that the results from the transcriptomics survey are highly reliable.

Discussion

Plant productivity is mainly affected by water deficit during different developmental stages, especially in arid and semi-arid areas. In the current study, we have compared drought tolerant and susceptible wheat genotypes at germination (Zadoks 05–09) and tillering (Zadoks 20–25) stages in terms of important root attributes and antioxidant capacity. We have further evaluated whole root transcriptomes to elucidate specific and general genes and protein pathways that contribute to water-deficit tolerance in plants, and which are interesting targets for molecular breeding programs or for manipulation with genetic engineering approaches.

Root Development and Radical Scavenging Activity of Contrasting Genotypes

Studying wheat landraces at an early stage using PEG 6000 at−10 bar showed a reduction in germination percentage, root length, and root biomass in both tolerant and susceptible genotypes, although tolerant genotypes showed smaller reductions in all traits. Tolerant genotypes were able to easily germinate even under severe water stress conditions (−10 and−12 bar) and showed higher germination speeds. Our study at tillering stage, showed that water deficit is followed by the accumulation of reactive oxygen species (ROS), such as H₂O₂, in leaves. H₂O₂ then further induce lipid peroxidation in phospholipid membranes yielding MDA (Czarnocka and Karpiński, 2018). MDA is the result of peroxidation of poly-unsaturated fatty acids which seems to be a signaling mechanism in response to stresses (De Dios Alché, 2019). In this study, the concentration of both H₂O₂ and MDA in PI627038 (T) and PI627299 (T) was steadily lower than PI624837 (S). After the production of ROS, plants deploy radical scavenging systems, either through enzymatic or non-enzymatic mechanisms. It has been shown, antioxidant enzymes, such as SOD, could function in radical scavenging efficiently enough to reduce oxidative damage (Ashraf, 2009; Caverzan et al., 2016; Kapoor et al., 2019). In our study, SOD showed higher activity after water deficit, although greater activities were observed in tolerant genotypes.

Another important plant response to stress is the accumulation of a range of metabolites e.g., amino acids. Proline is an amino acid that plays a crucial role in plant metabolism, especially under stressful conditions, as it contributes to antioxidant defense, signaling, metal chelating and most importantly acts as an osmolyte (Hayat et al., 2012; Kahraman et al., 2019). Our findings suggest that also proline content increase under stress treatment, compared to control conditions. We observed significant differences in the accumulation of proline between tolerant and susceptible genotypes under stress conditions. The overproduction of proline under stressful conditions would further convey stress tolerance through maintaining osmotic and ROS balance in cells, and ultimately prevents electrolyte leakage and oxidative damage in plants (Hayat et al., 2012; Liang et al., 2013; Forlani et al., 2019).

Previous studies suggest that plants usually extend lateral and embryo-formed roots when encountering drought stress conditions (Kashiwagi et al., 2005; Moumeni et al., 2011; Ahmadi et al., 2018; Bristiel et al., 2019). In the present study, tolerant genotypes showed deeper roots and more root biomass under both control and water deficit conditions and especially after 18 and 25 days of stress treatment. The total fresh biomass of roots and shoots generally declined during the water stress experiment, indicating the negative effect of water deficiency on the development of roots and shoots in wheat. The root/shoot ratio tends to increase under stress conditions, suggesting that plants invest more assimilates into developing the length and density of roots by extending lateral and embryo roots. However, the pilot investigation showed that more root biomass in tolerant genotypes was due to an increase in the number of lateral and adventitious root, rather than a production of deeper roots. In addition to environmental factors, root development mainly relies on hormone homeostasis, including auxin and abscisic acid (ABA) content. Auxin contributes to cell elongation by modulating the H⁺ secretion through adjusting PM H⁺-ATPase enzyme activity (Rober-Kleber et al., 2003; Fleming, 2006; Staal et al., 2011; Yuan et al., 2017; Ogura et al., 2019). Furthermore, ABA will accumulates in root tips after sensing drought and is crucial for maintaining root elongation (Sengupta et al., 2011; Xu et al., 2013). ABA functions in cell wall extensibility and regulate indole-3-acetic acid (IAA) transport, which further enhance root elongation and the generation of hairy roots (Yamaguchi and Sharp, 2010; Xu et al., 2013; Ma et al., 2019; Rosales et al., 2019).

Differential Gene Expression and Pathway Enrichment

In addition to understanding of hormone homeostasis, dissecting the expression pattern of other biological pathways involved in root development and in plant response to water deficit will further provide fundamental information for manipulating and exploiting these biological processes in order to improve the adaptability of plants. Our initial hypothesis was that differences in root morphology and antioxidant capacity between drought tolerant and susceptible genotypes is correlated with the differential expression of key genes, including transcription factors, and stress response pathways. We therefore used RNA-seq to enable us to detect DEGs and then used functional enrichment derived from these genes to build and characterize putative gene and protein interaction networks that contribute to root development. DEGs were extracted by pairwise comparisons between wheat genotypes, different drought conditions and time-points during the drought period. These analyses illuminate genetic differences between genotypes as well as responsive genes to drought. Finally, this should enable us to identify genes that are consistently drought-responsive under long-term water shortage in wheat.

We observed distinct differences in the transcriptome comparisons. A total of 14,187 genes were identified as being differentially expressed across comparisons, highlighting genes that potentially contribute to drought tolerance through up or down regulations. Some genes overlapped in both comparisons while others were unique to different tolerant genotypes. The number of down-regulated genes were generally higher than the number of up-regulated genes. Interestingly, we observed a significantly higher number of genotype-specific DEGs in the susceptible genotype under stress and non-stress conditions compared to the tolerant genotypes, in line with earlier observations by Fracasso et al. (2016) and Mia et al. (2020) who also reported similar expression pattern changes in susceptible materials. The DEGs that were down-regulated showed GO enrichment associated to catalytic activities, oxidative reactions, ion binding, protein kinase activity and antioxidant activity under both short-term and long-term water deficit. In contrast, up-regulated genes in the susceptible genotype under short and long-term water deficit were associated with transmembrane transporter activity, transporter activity and response to water processes. These findings suggest a lower focus on drought tolerance mechanisms in susceptible plants and that they invest more assimilate into the biosynthesis of secondary metabolites and metabolic pathways, resulting in higher energy consumption and consequently lower adaptation to drought (Mia et al., 2020).

Down-regulated genes under short-term water deficit in tolerant plants were significantly associated with nicotianamine synthase activity, tricarboxylic acid biosynthetic process and cellular biogenic amine biosynthetic process while up-regulated genes were mainly linked to peroxidase activity, oxidoreductase activity, antioxidant activity, heme binding, hydrogen peroxide catabolic process, response to oxidative stress, response to water and reactive oxygen species metabolic process. Interestingly, by increasing the time of water deficit from 9 to 18 days, we observed a downregulation in catalytic activity, oxidoreductase activity, heme binding, peroxidase activity, and antioxidant activity also in tolerant genotypes. At the same ime, up-regulated DEGs in tolerant genotypes were associated with phosphatase activity, protein serine/threonine phosphatase activity, phosphoric ester hydrolase activity, chitin-binding, DNA-binding transcription factor activity and sequence-specific DNA binding. These results suggest the existence of a water-deficit response threshold in these plants. In other words, tolerant genotypes would up-regulate genes and activate more pathways connected to water-deficit tolerance upon required and then they will active other regulatory gene networks involved in processes rather than tolerance. A further analysis of up-regulated genes showed that the most common responsive pathways to water deficit were the flavonoid pathway (pathway ID ko 00941), glutathione metabolism (ko 00480), phenylpropanoid biosynthesis (ko 00940), and diterpenoid biosynthesis (ko 00904). Accumulation of secondary metabolites such as flavonoids and phenylpropanoids occur when plants are exposed to drought and other stresses (Sharma et al., 2019). These products are part of the antioxidant defense system and scavenge reactive oxygen species (ROS) that are induced upon oxidative stress (Ma et al., 2014; Tattini et al., 2015; Niinemets, 2016). Phenylpropanoids are among the largest group of plant secondary metabolites and are formed from aromatic amino acids like phenylalanine and tyrosine. The final products in this pathway are flavonoids, phenolic acids, monolignols, coumarins, and stilbenes. Crucial regulatory enzymes, such as phenylalanine ammonia lyase (PAL), chalcone synthase (CHS), and transcription factors belonging to the MYB family, WD40-type or basic helix-loop-helix (bHLH) classes are well-known members of the phenylpropanoids pathway and could therefore be keye targets in genetic engineering experiments aimed at increasing wheat adaptation to drought (Falcone Ferreyra et al., 2012; Sharma et al., 2020; Verma et al., 2020). In our study, 33 genes encoding enzymes such as phenylalanine ammonia-lyase (K10775), cinnamoyl-CoA reductase 1 (K09753), peroxidase 70-like (K00430), and beta-glucosidase 1 (K05350) were up-regulated in the phenylpropanoid pathway. Seven genes encoding enzymes, such as ent-copalyl diphosphate synthase 1 (K04120), ent-kaurene oxidase 2 isoform X3 (K04122), ent-kaurene oxidase-like 5 (K21719), and 9-beta-pimara-7,15-diene oxidase (K16085) were upregulated in the diterpenoid pathway. In a previous study on roots of two pairs of drought-tolerant and susceptible rice near-isogenic lines, Moumeni et al. (2011) showed that a set of highly up regulated genes were mostly involved in secondary metabolism, response to stimulus, amino acid metabolism and signal transduction, while down-regulated genes were mainly responsible in photosynthesis and cell wall growth.

Predication of Regulatory Elements

Further analyses of regulatory factors highlighted the role of, for example, AP2/ERF-ERF, MYB-related, bHLH, B3, MADS-M-type, MYB, NAC, and WRKY as specific TFs affecting downstream genes and biosynthesis pathways. However, AP2/ERF-ERF, MYB-related, bHLH, and B3 transcription factor families were genotype-specific regulatory elements under water-deficit treatments in the present study. Transcription factors or trans-acting elements are the core of regulatory networks of downstream genes through binding to cis-acting elements in the promoters of target genes and thereby play a crucial role in modulating gene expression (Yamaguchi-Shinozaki and Shinozaki, 2006; Todaka et al., 2012; Nakashima et al., 2014; Xie et al., 2019). Among the TFs, the APETALA 2/ethylene-responsive element binding factor (AP2/ERF) which contain a conserved DNA-binding domain is well-known to harbor a large group of TFs in plants (Sakuma et al., 2002; Xie et al., 2019). This family includes four main subfamilies: the AP2, ERF, RAV, and DREB (dehydration-responsive element-binding protein), though the DREB subfamily is more well-characterized and family members are known to be involved in response to abiotic stress responses and especially drought stress (Yamaguchi-Shinozaki and Shinozaki, 2006; Mizoi et al., 2012; Xie et al., 2019). However, this response occurs via binding of these TFs to cis-acting dehydration-responsive element/C-repeat (DRE/CRT) domain in the promoters of target genes (Sakuma et al., 2002; Mizoi et al., 2012). Our results show that 19 and 16 TFs belonging to the AP2/ERF family members were involved in response to water deficit in tolerant genotypes.

The second important TF group we identified was the MYB-related subfamily with a single SANT domain (17 and 11 members in tolerant genotypes, respectively). Generally, MYB family genes can be classified into four subfamilies according to the number of repeats in the sequence (ranging from 1 to 4) including one repeat/MYB-related (1R-MYB), two repeats (R2R3-MYB), three repeats (3R-MYB), and four repeats or 4R-MYB (Ambawat et al., 2013; Xiong et al., 2014; Baldoni et al., 2015). Although it has been suggested R2R3-MYB type are primarily involved in response to various abiotic stresses (Xiong et al., 2014), our findings suggest that 1R-MYB could also contribute to drought-responsive regulatory networks. In line with other TFs, MYB proteins have two unique regions, a conserved N-terminal (DNA-binding domain) and a diverse C-terminal that determines the regulatory potential of the protein (Ambawat et al., 2013). Previous research have shown these TFs are mostly involved in the ABA-signaling pathway, phenylpropanoid metabolism, auxin signaling, secondary metabolite pathways and interestingly the activation of lateral root meristem under water deficit conditions (Seo and Park, 2009; Ambawat et al., 2013; Baldoni et al., 2015; Zhang et al., 2020).

Another important genotype-specific TF family we detected in the present study was bHLH or basic helix-loop-helix protein, which are known to contribute to a wide range of process in plants, such as plant metabolism, light signal transduction, photomorphogenesis and, more importantly, response to various stresses (Mao et al., 2017; Sun et al., 2018). This TF family contains a highly conserved basic/helix-loop-helix domain, which is formed by two different sections located at the N-terminal and C-terminal of the bHLH (Sun et al., 2018). The first part is a basic amino acid region (at N-terminus) and the second part is the helix-loop-helix region/HLH (at C-terminus). These TFs identify and bind to the target sequences using the basic amino acid region, while the HLH domain provides domain dimerization and further triggers protein-protein interactions (Niu et al., 2017; Tanabe et al., 2019). This in turn lead to the establishment of homodimeric or heterodimeric complexes which eventually regulates gene transcription (Sun et al., 2018; Tanabe et al., 2019). Previous studies suggest that members of this transcription factor family improve drought tolerance via inducing genes involved in ABA signaling, osmotic capacity, ROS scavenging, flavonoids pathway genes and launching other signaling cascades (Dong et al., 2014; Ji et al., 2016; Wang et al., 2016). Our results identified 12 and 8 TFs from the bHLH family that are involved in the drought response cascade in tolerant genotypes. Janiak et al. (2019) showed 13 transcription factors from DREB, WRKYs, AP2, NACs, MYB or MYB-related, bHLH, and bZIP TFs were the main TFs families affected under mild drought stress in barely as well.

Protein-Protein Interactions

The prediction of probable interactions between DEG-encoded proteins and transcription factors through co-expression, co-occurrence, fusions networks or other analyses provides a functionally characterized framework for further studies. Our results suggest that PPI networks in tolerant wheat landraces were relatively similar, although some unique differences were also observed. In both genotypes, the interactions between enzymes involved in phenolic component pathways, antioxidant defense system and glutathione metabolism were the main source of similarity. Two important proteins belonging to the CYP450 family were connected to PAL isozymes in the first cluster of PPI1 (Figure 8). Cytochromes P450 (CYPs) members form one of the largest and most important enzymatic protein families in plants (around 1% of total gene annotations in plant species) and they are involved in various processes, such as developmental, physiological, signaling, and defense systems in plants (Mizutani and Ohta, 2010; Wei and Chen, 2018; Magwanga et al., 2019). The crucial role of CYPs in primary and secondary metabolism occurs by catalyzing a broad range of monooxygenation/hydroxylation reactions (Pandian et al., 2020). The catalytic function of CYP73, CYP98, CYP84, CYP96B4, CYP709B1, CYP709B2, and CYP709B3 in the phenylpropanoid pathway, that provide structural phenolic compounds e.g., lignin and suberin, antioxidants components like polyphenols, and other protectant compounds such as flavonoids, has previously been reported for several plant species (Mizutani and Ohta, 2010; Mao et al., 2013; Wei and Chen, 2018). Duan et al. (2017) showed overexpression of SoCYP85A1 (Spinach Cytochrome p450) gene in transgenic tobacco could increase primary root length and lateral roots that resulted in enhanced drought tolerance of transgenic plants with respect to wild type. More recently, overexpression of CsCYT75B1 (Citrus Cytochrome P450) gene in transgenic Arabidopsis significantly increased the total contents of flavonoid and ultimate antioxidant capacity (Rao et al., 2020). In addition to mentioned CYP450 proteins in PPI1, another three CYP450 showed interactions with proteins belong to the iron/ascorbate-dependent oxidoreductase family in PPI2. In PPI3, we observed other isozymes of CYP450 interacting with beta-carotene hydroxylase A1, beta-carotene hydroxylase D1, and proteins in the carotenoid pathway. TaABA8OH2 or ABA 8'-hydroxylase belongs to the CYP707 family and catalyzes the conversion of ABA to 8'-hydroxy ABA and then phaseic acid (PA) through isomerization under drought conditions (Pandian et al., 2020). A dynamic balance between the biosynthesis and catabolism of ABA is interceded via ABA 8'-hydroxylase and determines changes in downstream processes controlled by ABA (Zheng et al., 2012; Pandian et al., 2020).

In the present study, GST enzymes were also significantly enriched for interacted with other proteins, in which GSTA1 showed a concrete relationship with a 219 amino acid (aa) protein as a mediator of other GSTs. These enzymes contribute to many functions related to stress tolerance through the detoxification of endobiotic and xenobiotic compounds. This occurs via conjugating reduced glutathione (tripeptide y-Glu-Cys-Gly) to hydrophobic electrophiles (Ji et al., 2010; Rezaei et al., 2013; Islam et al., 2019). On the other hand, drought-responsive transcription factors and protein kinases were significantly enriched in all networks. For TFs, MYB13-1, MADS box, GATA, bZIP, GRAS family, NAC2B, and the GATA transcription factors all demonstrated crucial interactions. The putative serine/threonine protein kinase (PK4), abscisic acid-inducible protein kinase, and SNF1-type serine-threonine protein kinase also showed important connections in several clusters.

Conclusion

We report putative mechanisms and biosynthetic pathways that are responsive in wheat roots in response to water-deficit stress. Investigation of contrasting drought tolerant and susceptible wheat genotypes allowed us to elucidate transcriptome changes, dynamic regulatory networks, and protein-protein interactions that are important under adverse conditions. In total, more than 14,000 genes were differentially expressed (DEG) and most of these genes mapped to the B genome of wheat. A significant proportion of these DEGs were enriched in flavonoid biosynthesis, phenylpropanoid biosynthesis, glutathione metabolism, diterpenoid biosynthesis, and hormone signaling pathways. Further, inspection of the DEGs identified 1,377 DEGs that potentially function as transcription factors and moderate downstream cascades in response to water deficit. These TFs were dominantly from AP2/ERF-ERF, MYB-related, bHLH, WRKY, and NAC families. We also modeled three protein-protein interaction networks (PPIs), utilizing significant KEGG proteins and putative TFs to predict and characterize interaction networks. It appears that several enzymes from the CYP450, TaABA8OH2, PAL, and GST families play conclusive roles in water-deficit tolerance in connection with MYB13-1, MADS-box, and NAC transcription factors and the protein kinase superfamily.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: ENA under the accession number PRJEB44859 (https://www.ebi.ac.uk/ena/browser/view/PRJEB44859).

Author Contributions

YR: designed and performed the experiments, analyzed the datasets, and wrote the manuscript. PI: assisted in designing the RNA-seq experiment and root studies panel, transcriptome analysis, and writing the manuscript. MB: assisted in biochemical and germination studies and edited the manuscript. HA: assisted in designing the RNA-seq experiment and edited the manuscript. AT: edited the manuscript. SK: assisted in biochemical and germination studies. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Swedish University of Agricultural Sciences, Uppsala, Sweden through a grant from the Swedish Research Council to PI. The funding body was not involved in the material preparation, designing the study, data analysis, or the writing the manuscript. The SNP&SEQ Platform was also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. All computer-based analyses were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Uppmax, partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank the University of Tehran and CIMMYT for providing plant materials in this study. Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala. The facility is part of the National Genomics Infrastructure (NGI) Sweden and Science for Life Laboratory.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.710867/full#supplementary-material

References

Ahmadi, J., Pour-Aboughadareh, A., Fabriki-Ourang, S., Mehrabi, A.-A., and Siddique, K. H. (2018). Screening wheat germplasm for seedling root architectural traits under contrasting water regimes: potential sources of variability for drought adaptation. Archiv. Agron. Soil Sci. 64, 1351–1365. doi: 10.1080/03650340.2018.1432855

CrossRef Full Text | Google Scholar

Ambawat, S., Sharma, P., Yadav, N. R., and Yadav, R. C. (2013). MYB transcription factor genes as regulators for plant responses: an overview. Physiol. Mol. Biol. Plants 19, 307–321. doi: 10.1007/s12298-013-0179-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashraf, M. (2009). Biotechnological approach of improving plant salt tolerance using antioxidants as markers. Biotechnol. Adv. 27, 84–93. doi: 10.1016/j.biotechadv.2008.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Baldoni, E., Genga, A., and Cominelli, E. (2015). Plant MYB transcription factors: their role in drought response mechanisms. Int. J. Mol. Sci. 16, 15811–15851. doi: 10.3390/ijms160715811

PubMed Abstract | CrossRef Full Text | Google Scholar

Bates, L. S., Waldren, R. P., and Teare, I. (1973). Rapid determination of free proline for water-stress studies. Plant Soil 39, 205–207. doi: 10.1007/BF00018060

CrossRef Full Text | Google Scholar

Bishopp, A., and Lynch, J. P. (2015). The hidden half of crop yields. Nat. Plants 1, 1–2. doi: 10.1038/nplants.2015.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Blighe, K., Rana, S., and Lewis, M. (2019). EnhancedVolcano: Publication-Ready Volcano Plots With Enhanced Colouring and Labeling. Available online at: https://bioconductor.org/packages/release/bioc/vignettes/EnhancedVolcano/inst/doc/EnhancedVolcano.html

Bradford, M. M. (1976). A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal. Biochem. 72, 248–254. doi: 10.1016/0003-2697(76)90527-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Bristiel, P., Roumet, C., Violle, C., and Volaire, F. (2019). Coping with drought: root trait variability within the perennial grass Dactylis glomerata captures a trade-off between dehydration avoidance and dehydration tolerance. Plant Soil 434, 327–342. doi: 10.1007/s11104-018-3854-8

CrossRef Full Text | Google Scholar

Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. (2009). BLAST+: architecture and applications. BMC Bioinform. 10, 1–9. doi: 10.1186/1471-2105-10-421

PubMed Abstract | CrossRef Full Text | Google Scholar

Caverzan, A., Casassola, A., and Brammer, S. P. (2016). Reactive oxygen species and antioxidant enzymes involved in plant tolerance to stress. Shanker AK & Shanker C. Abiotic and biotic stress in plants-Recent advances and future perspectives. Publisher InTech 2016, 463–480. doi: 10.5772/61368

CrossRef Full Text | Google Scholar

Comas, L., Becker, S., Cruz, V. M. V., Byrne, P. F., and Dierig, D. A. (2013). Root traits contributing to plant productivity under drought. Front. Plant Sci. 4:442. doi: 10.3389/fpls.2013.00442

PubMed Abstract | CrossRef Full Text | Google Scholar

Czarnocka, W., and Karpiński, S. (2018). Friend or foe? Reactive oxygen species production, scavenging and signaling in plant response to environmental stresses. Free Rad. Biol. Med. 122, 4–20. doi: 10.1016/j.freeradbiomed.2018.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, X., Sinharoy, S., Udvardi, M., and Zhao, P. X. (2013). PlantTFcat: an online plant transcription factor and transcriptional regulator categorization and analysis tool. BMC Bioinform. 14, 1–6. doi: 10.1186/1471-2105-14-321

PubMed Abstract | CrossRef Full Text | Google Scholar

De Dios Alché, J. (2019). A concise appraisal of lipid oxidation and lipoxidation in higher plants. Redox Biol. 23:101136. doi: 10.1016/j.redox.2019.101136

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhindsa, R. S., Plumb-Dhindsa, P. L., and Reid, D. M. (1982). Leaf senescence and lipid peroxidation: effects of some phytohormones, and scavengers of free radicals and singlet oxygen. Physiol. Plant. 56, 453–457. doi: 10.1111/j.1399-3054.1982.tb04539.x

CrossRef Full Text | Google Scholar

Dong, Y., Wang, C., Han, X., Tang, S., Liu, S., Xia, X., et al. (2014). A novel bHLH transcription factor PebHLH35 from Populus euphratica confers drought tolerance through regulating stomatal development, photosynthesis and growth in Arabidopsis. Biochem. Biophys. Res. Commun. 450, 453–458. doi: 10.1016/j.bbrc.2014.05.139

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, F., Ding, J., Lee, D., Lu, X., Feng, Y., and Song, W. (2017). Overexpression of SoCYP85A1, a spinach cytochrome p450 gene in transgenic tobacco enhances root development and drought stress tolerance. Front. Plant Sci. 8:1909. doi: 10.3389/fpls.2017.01909

PubMed Abstract | CrossRef Full Text | Google Scholar

Fahad, S., Bajwa, A. A., Nazir, U., Anjum, S. A., Farooq, A., Zohaib, A., et al. (2017). Crop production under drought and heat stress: plant responses and management options. Front. Plant Sci. 8:1147. doi: 10.3389/fpls.2017.01147

PubMed Abstract | CrossRef Full Text | Google Scholar

Falcone Ferreyra, M. L., Rius, S., and Casati, P. (2012). Flavonoids: biosynthesis, biological functions, and biotechnological applications. Front. Plant Sci. 3:222. doi: 10.3389/fpls.2012.00222

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, Y., Du, Y., Wang, J., Wu, A., Qiao, S., Xu, B., et al. (2017). Moderate drought stress affected root growth and grain yield in old, modern and newly released cultivars of winter wheat. Front. Plant Sci. 8:672. doi: 10.3389/fpls.2017.00672

PubMed Abstract | CrossRef Full Text | Google Scholar

Fleming, A. J. (2006). Plant signalling: the inexorable rise of auxin. Trends Cell Biol. 16, 397–402. doi: 10.1016/j.tcb.2006.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Forlani, G., Trovato, M., Funck, D., and Signorelli, S. (2019). “Regulation of proline accumulation and its molecular and physiological functions in stress defence,” in Osmoprotectant-Mediated Abiotic Stress Tolerance in Plants, (eds) M. A. Hossain, V. Kumar, D. J. Burritt, M. Fujita, and P. S. A. Mäkelä (Springer), 73–97. doi: 10.1007/978-3-030-27423-8_3

CrossRef Full Text | Google Scholar

Fracasso, A., Trindade, L. M., and Amaducci, S. (2016). Drought stress tolerance strategies revealed by RNA-Seq in two sorghum genotypes with contrasting WUE. BMC Plant Biol. 16, 1–18. doi: 10.1186/s12870-016-0800-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Giehl, R. F., Gruber, B. D., and Von Wirén, N. (2014). It's time to make changes: modulation of root system architecture by nutrient signals. J. Exp. Bot. 65, 769–778. doi: 10.1093/jxb/ert421

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayat, S., Hayat, Q., Alyemeni, M. N., Wani, A. S., Pichtel, J., and Ahmad, A. (2012). Role of proline under changing environments: a review. Plant Signal. Behav. 7, 1456–1466. doi: 10.4161/psb.21949

PubMed Abstract | CrossRef Full Text | Google Scholar

Heath, R. L., and Packer, L. (1968). Photoperoxidation in isolated chloroplasts: I. Kinetics and stoichiometry of fatty acid peroxidation. Archiv. Biochem. Biophys. 125, 189–198. doi: 10.1016/0003-9861(68)90654-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernandez-Ochoa, I. M., Asseng, S., Kassie, B. T., Xiong, W., Robertson, R., Pequeno, D. N. L., et al. (2018). Climate change impact on Mexico wheat production. Agri. Forest Meteorol. 263, 373–387. doi: 10.1016/j.agrformet.2018.09.008

CrossRef Full Text | Google Scholar

Hochholdinger, F., and Tuberosa, R. (2009). Genetic and genomic dissection of maize root development and architecture. Curr. Opin. Plant Biol. 12, 172–177. doi: 10.1016/j.pbi.2008.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Islam, S., Sajib, S. D., Jui, Z. S., Arabia, S., Islam, T., and Ghosh, A. (2019). Genome-wide identification of glutathione S-transferase gene family in pepper, its classification, and expression profiling under different anatomical and environmental conditions. Sci. Rep. 9, 1–15. doi: 10.1038/s41598-019-45320-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Janiak, A., Kwasniewski, M., Sowa, M., Kuczyńska, A., Mikołajczak, K., Ogrodowicz, P., et al. (2019). Insights into barley root transcriptome under mild drought stress with an emphasis on gene expression regulatory mechanisms. Int. J. Mol. Sci. 20:6139. doi: 10.3390/ijms20246139

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, W., Zhu, Y., Li, Y., Yang, L., Zhao, X., Cai, H., et al. (2010). Over-expression of a glutathione S-transferase gene, GsGST, from wild soybean (Glycine soja) enhances drought and salt tolerance in transgenic tobacco. Biotechnol. Lett. 32, 1173–1179. doi: 10.1007/s10529-010-0269-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, X., Nie, X., Liu, Y., Zheng, L., Zhao, H., Zhang, B., et al. (2016). A bHLH gene from Tamarix hispida improves abiotic stress tolerance by enhancing osmotic potential and decreasing reactive oxygen species accumulation. Tree Physiol. 36, 193–207. doi: 10.1093/treephys/tpv139

PubMed Abstract | CrossRef Full Text | Google Scholar

Kahraman, M., Sevim, G., and Bor, M. (2019). “The role of proline, glycinebetaine, and trehalose in stress-responsive gene expression,” in Osmoprotectant-Mediated Abiotic Stress Tolerance in Plants, (eds) M. Hossain, V. Kumar, D. Burritt, M. Fujita, and P. Mäkelä (Cham: Springer), 241–256. doi: 10.1007/978-3-030-27423-8_11

CrossRef Full Text | Google Scholar

Kapoor, D., Singh, S., Kumar, V., Romero, R., Prasad, R., and Singh, J. (2019). Antioxidant enzymes regulation in plants in reference to reactive oxygen species (ROS) and reactive nitrogen species (RNS). Plant Gene 19:100182. doi: 10.1016/j.plgene.2019.100182

CrossRef Full Text | Google Scholar

Kashiwagi, J., Krishnamurthy, L., Upadhyaya, H. D., Krishna, H., Chandra, S., Vadez, V., et al. (2005). Genetic variability of drought-avoidance root traits in the mini-core germplasm collection of chickpea (Cicer arietinum L.). Euphytica 146, 213–222. doi: 10.1007/s10681-005-9007-1

CrossRef Full Text | Google Scholar

Khan, M., Gemenet, D. C., and Villordon, A. (2016). Root system architecture and abiotic stress tolerance: current knowledge in root and tuber crops. Front. Plant Sci. 7:1584. doi: 10.3389/fpls.2016.01584

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, X., Zhang, L., Natarajan, S. K., and Becker, D. F. (2013). Proline mechanisms of stress survival. Antioxid. Redox Signal. 19, 998–1011. doi: 10.1089/ars.2012.5074

PubMed Abstract | CrossRef Full Text | Google Scholar

Loreto, F., and Velikova, V. (2001). Isoprene produced by leaves protects the photosynthetic apparatus against ozone damage, quenches ozone products, and reduces lipid peroxidation of cellular membranes. Plant Physiol. 127, 1781–1787. doi: 10.1104/pp.010497

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M., Anders, S., and Huber, W. (2014). Differential analysis of count data–the DESeq2 package. Genome Biol. 15:10.1186. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, J. (1995). Root architecture and plant productivity. Plant Physiol. 109:7. doi: 10.1104/pp.109.1.7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, J. P., and Brown, K. M. (2012). New roots for agriculture: exploiting the root phenome. Philos. Trans. Royal Soc. B Biol. Sci. 367, 1598–1604. doi: 10.1098/rstb.2011.0243

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, D., Sun, D., Wang, C., Li, Y., and Guo, T. (2014). Expression of flavonoid biosynthesis genes and accumulation of flavonoid in wheat leaves in response to drought stress. Plant Physiol. Biochem. 80, 60–66. doi: 10.1016/j.plaphy.2014.03.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J., Li, R., Wang, H., Li, D., Wang, X., Zhang, Y., et al. (2017). Transcriptomics analyses reveal wheat responses to drought stress during reproductive stages under field conditions. Front. Plant Sci. 8:592. doi: 10.3389/fpls.2017.00592

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, X., Zhang, X., Yang, L., Tang, M., Wang, K., Wang, L., et al. (2019). Hydrogen peroxide plays an important role in PERK4-mediated abscisic acid-regulated root growth in Arabidopsis. Funct. Plant Biol. 46, 165–174. doi: 10.1071/FP18219

PubMed Abstract | CrossRef Full Text | Google Scholar

Maccaferri, M., El-Feki, W., Nazemi, G., Salvi, S., Canè, M. A., Colalongo, M. C., et al. (2016). Prioritizing quantitative trait loci for root system architecture in tetraploid wheat. J. Exp. Bot. 67, 1161–1178. doi: 10.1093/jxb/erw039

PubMed Abstract | CrossRef Full Text | Google Scholar

Magwanga, R. O., Lu, P., Kirungu, J. N., Dong, Q., Cai, X., Zhou, Z., et al. (2019). Knockdown of cytochrome P450 genes Gh_D07G1197 and Gh_A13G2057 on chromosomes D07 and A13 reveals their putative role in enhancing drought and salt stress tolerance in Gossypium hirsutum. Genes 10:226. doi: 10.3390/genes10030226

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, G., Seebeck, T., Schrenker, D., and Yu, O. (2013). CYP709B3, a cytochrome P450 monooxygenase gene involved in salt tolerance in Arabidopsis thaliana. BMC Plant Biol. 13, 1–13. doi: 10.1186/1471-2229-13-169

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, K., Dong, Q., Li, C., Liu, C., and Ma, F. (2017). Genome wide identification and characterization of apple bHLH transcription factors and expression analysis in response to drought and salt stress. Front. Plant Sci. 8:480. doi: 10.3389/fpls.2017.00480

PubMed Abstract | CrossRef Full Text | Google Scholar

Mia, M. S., Liu, H., Wang, X., Zhang, C., and Yan, G. (2020). Root transcriptome profiling of contrasting wheat genotypes provides an insight to their adaptive strategies to water deficit. Sci. Rep. 10, 1–11. doi: 10.1038/s41598-020-61680-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizoi, J., Shinozaki, K., and Yamaguchi-Shinozaki, K. (2012). AP2/ERF family transcription factors in plant abiotic stress responses. Biochim. Biophys. Acta 1819, 86–96. doi: 10.1016/j.bbagrm.2011.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizutani, M., and Ohta, D. (2010). Diversification of P450 genes during land plant evolution. Annu. Rev. Plant Biol. 61, 291–315. doi: 10.1146/annurev-arplant-042809-112305

PubMed Abstract | CrossRef Full Text | Google Scholar

Moumeni, A., Satoh, K., Kondoh, H., Asano, T., Hosaka, A., Venuprasad, R., et al. (2011). Comparative analysis of root transcriptome profiles of two pairs of drought-tolerant and susceptible rice near-isogenic lines under different drought stress. BMC Plant Biol. 11, 1–17. doi: 10.1186/1471-2229-11-174

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakashima, K., Yamaguchi-Shinozaki, K., and Shinozaki, K. (2014). The transcriptional regulatory network in the drought response and its crosstalk in abiotic stress responses including drought, cold, and heat. Front. Plant Sci. 5:170. doi: 10.3389/fpls.2014.00170

PubMed Abstract | CrossRef Full Text | Google Scholar

Nieuwland, J., Maughan, S., Dewitte, W., Scofield, S., Sanz, L., and Murray, J. A. (2009). The D-type cyclin CYCD4; 1 modulates lateral root density in Arabidopsis by affecting the basal meristem region. Proc. Nat. Acad. Sci. U. S. A. 106, 22528–22533. doi: 10.1073/pnas.0906354106

PubMed Abstract | CrossRef Full Text | Google Scholar

Niinemets, Ü. (2016). Uncovering the hidden facets of drought stress: secondary metabolites make the difference. Tree Physiol. 36, 129–132. doi: 10.1093/treephys/tpv128

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, X., Guan, Y., Chen, S., and Li, H. (2017). Genome-wide analysis of basic helix-loop-helix (bHLH) transcription factors in Brachypodium distachyon. BMC Genomics 18, 1–20. doi: 10.1186/s12864-017-4044-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogura, T., Goeschl, C., Filiault, D., Mirea, M., Slovak, R., Wolhrab, B., et al. (2019). Root system depth in Arabidopsis is shaped by EXOCYST70A3 via the dynamic modulation of auxin transport. Cell 178, 400–412.e416. doi: 10.1016/j.cell.2019.06.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Overvoorde, P., Fukaki, H., and Beeckman, T. (2010). Auxin control of root development. Cold Spring Harb. Perspect. Biol. 2:a001537. doi: 10.1101/cshperspect.a001537

PubMed Abstract | CrossRef Full Text | Google Scholar

Paez-Garcia, A., Motes, C. M., Scheible, W.-R., Chen, R., Blancaflor, E. B., and Monteros, M. J. (2015). Root traits and phenotyping strategies for plant improvement. Plants 4, 334–355. doi: 10.3390/plants4020334

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandian, B. A., Sathishraj, R., Djanaguiraman, M., Prasad, P., and Jugulam, M. (2020). Role of cytochrome P450 enzymes in plant stress response. Antioxidants 9:454. doi: 10.3390/antiox9050454

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, J. H., Sun, D., and Nevo, E. (2011). Domestication evolution, genetics and genomics in wheat. Mol. Breeding 28, 281–301. doi: 10.1007/s11032-011-9608-4

CrossRef Full Text | Google Scholar

Pessarakli, M. (2019). Handbook of Plant and Crop Stress. Boca Raton, FL: CRC Press. doi: 10.1201/9781351104609

CrossRef Full Text | Google Scholar

Rahimi, Y., Bihamta, M. R., Taleei, A., Alipour, H., and Ingvarsson, P. K. (2019). Genome-wide association study of agronomic traits in bread wheat reveals novel putative alleles for future breeding programs. BMC Plant Biol. 19, 1–19. doi: 10.1186/s12870-019-2165-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Rao, M. J., Xu, Y., Tang, X., Huang, Y., Liu, J., Deng, X., et al. (2020). CsCYT75B1, a citrus cytochrome P450 gene, is involved in accumulation of antioxidant flavonoids and induces drought tolerance in transgenic Arabidopsis. Antioxidants 9:161. doi: 10.3390/antiox9020161

PubMed Abstract | CrossRef Full Text | Google Scholar

Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson, H., et al. (2019). g: Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucl. Acids Res. 47, W191–W198. doi: 10.1093/nar/gkz369

PubMed Abstract | CrossRef Full Text | Google Scholar

Rezaei, M. K., Shobbar, Z.-S., Shahbazi, M., Abedini, R., and Zare, S. (2013). Glutathione S-transferase (GST) family in barley: identification of members, enzyme activity, and gene expression pattern. J. Plant Physiol. 170, 1277–1284. doi: 10.1016/j.jplph.2013.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Rober-Kleber, N., Albrechtová, J. T., Fleig, S., Huck, N., Michalke, W., Wagner, E., et al. (2003). Plasma membrane H+-ATPase is involved in auxin-mediated cell elongation during wheat embryo development. Plant Physiol. 131, 1302–1312. doi: 10.1104/pp.013466

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosales, M. A., Maurel, C., and Nacry, P. (2019). Abscisic acid coordinates dose-dependent developmental and hydraulic responses of roots to water deficit. Plant Physiol. 180, 2198–2211. doi: 10.1104/pp.18.01546

PubMed Abstract | CrossRef Full Text | Google Scholar

Saddique, Q., Li Liu, D., Wang, B., Feng, P., He, J., Ajaz, A., et al. (2020). Modelling future climate change impacts on winter wheat yield and water use: a case study in Guanzhong Plain, northwestern China. Eur. J. Agron. 119:126113. doi: 10.1016/j.eja.2020.126113

CrossRef Full Text | Google Scholar

Sakuma, Y., Liu, Q., Dubouzet, J. G., Abe, H., Shinozaki, K., and Yamaguchi-Shinozaki, K. (2002). DNA-binding specificity of the ERF/AP2 domain of Arabidopsis DREBs, transcription factors involved in dehydration-and cold-inducible gene expression. Biochem. Biophys. Res. Commun. 290, 998–1009. doi: 10.1006/bbrc.2001.6299

PubMed Abstract | CrossRef Full Text | Google Scholar

Sengupta, D., Kannan, M., and Reddy, A. R. (2011). A root proteomics-based insight reveals dynamic regulation of root proteins under progressive drought stress and recovery in Vigna radiata (L.) Wilczek. Planta 233, 1111–1127. doi: 10.1007/s00425-011-1365-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Seo, P. J., and Park, C.-M. (2009). Auxin homeostasis during lateral root development under drought condition. Plant Signal. Behav. 4, 1002–1004. doi: 10.4161/psb.4.10.9716

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, A., Badola, P. K., Bhatia, C., Sharma, D., and Trivedi, P. K. (2020). Primary transcript of miR858 encodes regulatory peptide and controls flavonoid biosynthesis and development in Arabidopsis. Nat. Plants 6, 1262–1274. doi: 10.1038/s41477-020-00769-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, A., Shahzad, B., Rehman, A., Bhardwaj, R., Landi, M., and Zheng, B. (2019). Response of phenylpropanoid pathway and the role of polyphenols in plants under abiotic stress. Molecules 24:2452. doi: 10.3390/molecules24132452

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S., and De Smet, I. (2012). Root system architecture: insights from Arabidopsis and cereal crops. Philos. Trans. R Soc. Lond. B Biol. Sci. 367, 1441–1452. doi: 10.1098/rstb.2011.0234

PubMed Abstract | CrossRef Full Text | Google Scholar

Staal, M., De Cnodder, T., Simon, D., Vandenbussche, F., Van Der Straeten, D., Verbelen, J.-P., et al. (2011). Apoplastic alkalinization is instrumental for the inhibition of cell elongation in the Arabidopsis root by the ethylene precursor 1-aminocyclopropane-1-carboxylic acid. Plant Physiol. 155, 2049–2055. doi: 10.1104/pp.110.168476

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinegger, M., Meier, M., Mirdita, M., Vöhringer, H., Haunsberger, S. J., and Söding, J. (2019). HH-suite3 for fast remote homology detection and deep protein annotation. BMC Bioinform. 20, 1–15. doi: 10.1186/s12859-019-3019-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Studer, G., Rempfer, C., Waterhouse, A. M., Gumienny, R., Haas, J., and Schwede, T. (2020). QMEANDisCo—distance constraints applied on model quality estimation. Bioinformatics 36, 1765–1771. doi: 10.1093/bioinformatics/btz828

PubMed Abstract | CrossRef Full Text | Google Scholar

Studer, G., Tauriello, G., Bienert, S., Biasini, M., Johner, N., and Schwede, T. (2021). ProMod3—A versatile homology modelling toolbox. PLoS Comput. Biol. 17:e1008667. doi: 10.1371/journal.pcbi.1008667

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Wang, Y., and Sui, N. (2018). Transcriptional regulation of bHLH during plant response to stress. Biochem. Biophys. Res. Commun. 503, 397–401. doi: 10.1016/j.bbrc.2018.07.123

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., et al. (2021). The STRING database in 2021: customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 49, D605–D612. doi: 10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanabe, N., Noshi, M., Mori, D., Nozawa, K., Tamoi, M., and Shigeoka, S. (2019). The basic helix-loop-helix transcription factor, bHLH11 functions in the iron-uptake system in Arabidopsis thaliana. J. Plant Res. 132, 93–105. doi: 10.1007/s10265-018-1068-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Tattini, M., Loreto, F., Fini, A., Guidi, L., Brunetti, C., Velikova, V., et al. (2015). Isoprenoids and phenylpropanoids are part of the antioxidant defense orchestrated daily by drought-stressed P latanus× acerifolia plants during Mediterranean summers. New Phytol. 207, 613–626. doi: 10.1111/nph.13380

PubMed Abstract | CrossRef Full Text | Google Scholar

Team, R. (2020). RStudio: Integrated Development for R. Boston, MA: RStudio, PBC. Available online at: http://www.rstudio.com

Todaka, D., Nakashima, K., Shinozaki, K., and Yamaguchi-Shinozaki, K. (2012). Toward understanding transcriptional regulatory networks in abiotic stress responses and tolerance in rice. Rice 5, 1–9. doi: 10.1186/1939-8433-5-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Verma, M. L., Sharma, S., Saini, R., Rani, V., and Kushwaha, R. (2020). “Bioflavonoids: synthesis, functions and biotechnological applications,” in Biotechnological Production of Bioactive Compounds, (eds) M. L. Verma, and A. K. Chandel (Elsevier), 69–105. doi: 10.1016/B978-0-444-64323-0.00003-5

CrossRef Full Text | Google Scholar

Wang, F., Zhu, H., Chen, D., Li, Z., Peng, R., and Yao, Q. (2016). A grape bHLH transcription factor gene, VvbHLH1, increases the accumulation of flavonoids and enhances salt and drought tolerance in transgenic Arabidopsis thaliana. Plant Cell Tissue Org. Cult. 125, 387–398. doi: 10.1007/s11240-016-0953-1

CrossRef Full Text | Google Scholar

Wang, Z., Gerstein, M., and Snyder, M. (2009). RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57–63. doi: 10.1038/nrg2484

PubMed Abstract | CrossRef Full Text | Google Scholar

Wasson, A. P., Richards, R., Chatrath, R., Misra, S., Prasad, S. S., Rebetzke, G., et al. (2012). Traits and selection strategies to improve root systems and water uptake in water-limited wheat crops. J. Exp. Bot. 63, 3485–3498. doi: 10.1093/jxb/ers111

PubMed Abstract | CrossRef Full Text | Google Scholar

Waterhouse, A., Bertoni, M., Bienert, S., Studer, G., Tauriello, G., Gumienny, R., et al. (2018). SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 46, W296–W303. doi: 10.1093/nar/gky427

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, K., and Chen, H. (2018). Global identification, structural analysis and expression characterization of cytochrome P450 monooxygenase superfamily in rice. BMC Genom. 19, 1–18. doi: 10.1186/s12864-017-4425-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., et al. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. doi: 10.1093/nar/gkr483

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Z., Nolan, T. M., Jiang, H., and Yin, Y. (2019). AP2/ERF transcription factor regulatory networks in hormone and abiotic stress responses in Arabidopsis. Front. Plant Sci. 10:228. doi: 10.3389/fpls.2019.00228

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, H., Guo, H., Xie, Y., Zhao, L., Gu, J., Zhao, S., et al. (2017). RNAseq analysis reveals pathways and candidate genes associated with salinity tolerance in a spaceflight-induced wheat mutant. Sci. Rep. 7, 1–13. doi: 10.1038/s41598-017-03024-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, H., Li, J., Liu, P., Duan, J., Zhao, Y., Guo, X., et al. (2014). Overexpression of OsMYB48-1, a novel MYB-related transcription factor, enhances drought and salinity tolerance in rice. PLoS ONE 9:e92913. doi: 10.1371/journal.pone.0092913

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, W., Jia, L., Shi, W., Liang, J., Zhou, F., Li, Q., et al. (2013). Abscisic acid accumulation modulates auxin transport in the root tip to enhance proton secretion for maintaining root growth under moderate water stress. New Phytol. 197, 139–150. doi: 10.1111/nph.12004

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaguchi, M., and Sharp, R. E. (2010). Complexity and coordination of root growth at low water potentials: recent advances from transcriptomic and proteomic analyses. Plant Cell Environ. 33, 590–603. doi: 10.1111/j.1365-3040.2009.02064.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaguchi-Shinozaki, K., and Shinozaki, K. (2006). Transcriptional regulatory networks in cellular responses and tolerance to dehydration and cold stresses. Annu. Rev. Plant Biol. 57, 781–803. doi: 10.1146/annurev.arplant.57.032905.105444

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, W., Zhang, D., Song, T., Xu, F., Lin, S., Xu, W., et al. (2017). Arabidopsis plasma membrane H+-ATPase genes AHA2 and AHA7 have distinct and overlapping roles in the modulation of root tip H+ efflux in response to low-phosphorus stress. J. Exp. Bot. 68, 1731–1741. doi: 10.1093/jxb/erx040

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P., Wang, R., Yang, X., Ju, Q., Li, W., Lü, S., et al. (2020). The R2R3-MYB transcription factor AtMYB49 modulates salt tolerance in Arabidopsis by modulating the cuticle formation and antioxidant defence. Plant Cell Environ. 43, 1925–1943. doi: 10.1111/pce.13784

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Y., Huang, Y., Xian, W., Wang, J., and Liao, H. (2012). Identification and expression analysis of the Glycine max CYP707A gene family in response to drought and salt stresses. Ann. Bot. 110, 743–756. doi: 10.1093/aob/mcs133

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, Y., Jiao, C., Sun, H., Rosli, H. G., Pombo, M. A., Zhang, P., et al. (2016). iTAK: a program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases. Mol. Plant 9, 1667–1670. doi: 10.1016/j.molp.2016.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar