RNA-Seq analysis revealed genes associated with drought stress response in kabuli chickpea (Cicer arietinum L.).

Drought is the most important constraint that effects chickpea production globally. RNA-Seq has great potential to dissect the molecular mechanisms of tolerance to environmental stresses. Transcriptome profiles in roots and shoots of two contrasting Iranian kabuli chickpea genotypes (Bivanij and Hashem) were investigated under water-limited conditions at early flowering stage using RNA-Seq approach. A total of 4,572 differentially expressed genes (DEGs) were identified. Of these, 261 and 169 drought stress responsive genes were identified in the shoots and the roots, respectively, and 17 genes were common in the shoots and the roots. Gene Ontology (GO) analysis revealed several sub-categories related to the stress, including response to stress, defense response and response to stimulus in the tolerant genotype Bivanij as compared to the sensitive genotype Hashem under drought stress. In addition, several Transcription factors (TFs) were identified in major metabolic pathways such as, ABA, proline and flavonoid biosynthesis. Furthermore, a number of the DEGs were observed in "QTL-hotspot" regions which were reported earlier in chickpea. Drought tolerance dissection in the genotypes revealed that the genes and the pathways involved in shoots of Bivanij were the most important factor to make a difference between the genotypes for drought tolerance. The identified TFs in the experiment, particularly those which were up-regulated in shoots of Bivanij during drought stress, were potential candidates for enhancing tolerance to drought.

Introduction

Chickpea (Cicer arietinum L.) is the second most important cool season grain legume crop cultivated on 13.98 M ha globally in 59 countries with a total production of 13.73 M tons [1]. It is a self-pollinated crop with a basic chromosome number eight (2n = 2x = 16) and genome size of ~738 Mb [2]. It is cultivated on residual soil moisture by resource poor farmers especially in the arid and semi-arid regions of South Asia and Sub-Saharan Africa. Major chickpea producing countries include India, Australia, Turkey, Myanmar, Pakistan, Ethiopia, Iran, etc. In terms of chickpea production Iran ranks seventh, while productivity (500 kg/ha) is much less than the world’s average productivity. Desi (light to dark brown in color) and kabuli (white or beige colored seed) are two types of chickpea cultivated globally. Like other legumes, chickpea enhances the soil fertility by fixing atmospheric nitrogen and also serves as a valuable source of proteins, vitamins and essential amino acids [3]. The productivity of chickpea has been low due to several abiotic and biotic stresses.

Among abiotic stresses, drought alone leads to >50% annual yield losses in chickpea. End season drought or terminal drought delays flowering and reduces growth duration particularly at reproductive stage and finally leads to yield losses [4-5]. Considering drastic effects of drought at reproductive stages, identification of drought tolerance genes during flowering is essential for understanding molecular mechanisms of tolerance and plant adaptation. Advances in genomics research during last decade facilitated development of several molecular markers [6], genetic maps [7] which increased our understanding of the genetics of complex traits like drought. Several approaches like linkage mapping [8] and genome-wide association mapping [9] were employed in past to dissect drought tolerance mechanism in chickpea. Further, efforts were also made to fine map this genomic region to few kilobases [10-Sci Rep. 2015 ">11]. Unprecedented developments in next generation sequencing (NGS) technologies not only decoded the genome architecture of several important crop plants [12-13] but also provided insights into the proteome and metabolome of model plant species [14]. In the case of chickpea using NGS technologies, draft genomes of both desi and kabuli chickpea have been deciphered [2,15]. In addition, whole genome resequencing (WGRS) of parental lines of mapping populations provided large scale single nucleotide polymorphisms (SNPs) data for trait improvement [16]. Insights into the temporal and spatial trends in diversity of released varieties have been reported using WGRS approach [17]. Although large amount of transcriptome data on desi genotypes were generated in recent years [18-22], specific role of stress responsive genes and their interactions are still unknown in kabuli chickpea genotypes.

RNA-Seq approach, that uses deep sequencing technologies, is a cost effective way for transcriptome studies [23]. Furthermore, RNA-Seq has been the technology of choice for identification of novel genes and isoforms, detection of variants including expressed SNPs, InDels, SSRs and gene fusion events [24]. RNA-Seq approach enabled a clear understanding of differentially expressed genes and physiological responses under drought stress in sorghum [25]. In case of barley, genes responsible for reproductive success under drought stress were reported [26]. High temperature and drought stress-related genes were reported in case of Brassica juncea [27].

In the present study, RNA-Seq was employed to investigate transcriptome profiles in drought-responsive contrasting genotypes of Iranian kabuli chickpea under drought stress in root and shoot tissues at early flowering stage. RNA-Seq on Illumina platform provided a thorough scenario on the whole chickpea transcriptome in response to drought stress. Several categories of key genes involved in drought response have been identified.

Materials and methods

Plant Material, growth conditions and drought stress treatment

Two contrasting Iranian chickpea kabuli genotypes Bivanij (drought tolerant) and Hashem (drought sensitive) were used in this study. The seeds were obtained from Sararood Rainfed Agriculture Research Station, Kermanshah, Iran. Two seeds were sown in pots, (0.21 m deep with 0.25 m diameter) containing 9 kg of vertisol soil and fertilized with sterilized farm yard manure as one part for 20 parts soil (v/v), in a greenhouse (ICRISAT, Hyderabad, India). The plants were thinned to one healthy plant per pot after 8 days of germination. Five biological replicates of each genotype and treatment were maintained under glasshouse conditions. The temperature inside the greenhouse were at a maximum of 24–28°C and a minimum of 12–22°C, and a day-time relative humidity of 30–70% with a maximum solar radiation incidence of >1200 μE m−2 s−1. The plants were grown under well-watered conditions for 4 weeks before initiating a dry-down experiment. One day before imposing stress, all the pots were saturated with water and allowed to drain excess water in 24 hours to maintain the field capacity so that the soil moisture amount in each pot were uniform. Then each pot was weighed to know the amount of water at field capacity. All pots were weighted every day on a sensitive balance with a resolution of plus or minus 0.1 g. The control plants were then watered to restore the soil to about 80% of field capacity by adding a pre-calculated amount of water while the stress plants were remained stressed for 15 days until the stress level reached to 20% of their field capacity. All the root and the shoot tissue samples were harvested in at least three biological replicates. The tissues were washed thoroughly with 0.1% DEPC water, frozen in liquid nitrogen and stored at -80°C until RNA extraction.

RNA extraction, Illumina sequencing and data quality control

Total RNA was extracted from 24 samples (two genotypes × two treatments × two tissues × three biological replicates) using NucleoSpin RNA Plant kit (MACHEREY-NAGEL GmbH & Co. KG, Germany), according to manufacturer’s instructions. Quantity and quality of the RNA was determined using Nanodrop 8000 Spectrophotometer (ThemoScientific, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, USA). The RNA (with RIN ≥8) was pooled from the three biological replicates of each sample and then a total of eight samples (namely Bivanij root drought control (BRDC), Bivanij shoot drought control (BSDC), Bivanij root drought stress (BRDS), Bivanij shoot drought stress (BSDS), Hashem root drought control (HRDC), Hashem shoot drought control (HSDC), Hashem root drought stress (HRDS) and Hashem shoot drought stress (HSDS) were processed further using Illumina TrueSeq RNA Sample Prep kit for cDNA library preparation. In brief, it includes purification of the polyadenylated RNA, the RNA fragmentation, cDNA synthesis and adaptor and barcode ligation. After cluster generation using Illumina cBot, the libraries were sequenced on Illumina HiSeq 2500 platform to generate >30 million 125 bp length paired-end (PE) reads for each sample. The high-quality reads were obtained after several steps of quality checks which included trimming, removal of adaptor/primer and low quality reads using Trimmomatic v 0.35 [28] and NGS-QCbox [29].

Reference-based mapping and assembly

The transcriptome data was analyzed using Tuxedo pipeline [30]. The filtered high-quality reads of all samples were mapped on kabuli chickpea reference genome (v1.0) using TopHat v2.1.0 [31] with default parameters. The mapped reads from each sample along with the genome GFF were used to perform reference annotation based transcript (RABT) assembly using Cufflinks v2.0.2 [32,33]. These Cufflink assemblies were then compared and merged using Cuffmerge to remove transfags and generate a consensus assembly for downstream analysis [32].

Identification of differentially expressed genes (DEGs), GO and pathway analysis

The normalized expression of all the samples was estimated in FPKM (fragments per kilobase of exon per million fragments mapped). The changes in the relative abundance of the genes between the different genotypes/treatments/tissues were estimated usingCuffdiff [34]. The genes exhibiting significant differences (at least two-fold change and a FPKM of > = 3 for either of the sample) were considered to be differentially expressed genes. Further, CummeRbund was used to plot and visualize abundance and differential expression results from Cuffdiff. The heatmaps showing the expression profiles were generated using MultiExperiment Viewer v4.9 [35]. To understand the functional classification of DEGs, the DEGS were subjected to BLASTX similarity search against NCBI non-redundant protein database with an E-value cut off of ≤10−5 followed by their annotation using Blast2go [36]. Pathway analysis of the DEGs was carried out using KEGG database. For identification of transcription factor encoding genes, the DEGs were searched against Plant transcription factor database (PlantTFDB 4.0) with an E-value cut off of ≤10−5.

Quantitative real-time PCR (qRT-PCR) analysis

To validate RNA-Seq results, 12 genes (S1 Table) selected from the DEGs analysis were subjected to quantitative real-time PCR (qRT-PCR). The primer pairs were designed using Primer Express (v3.0) software (Applied Biosystems, Foster City, CA). The RNA samples used for sequencing (mentioned above) from three biological replicates were used for the validation. The total RNA was reverse-transcribed using SuperScript III First-Strand synthesis kit (Invitrogen, life technologies, USA). After quality control by Nanodrop spectrophotometer, the cDNA samples were diluted and normalized with the house keeping gene Glyceraldehyde 3-phosphate dehydrogenase (GAPDH). The qRT-PCR analysis was performed using Applied Biosystems 7500 Real-Time PCR System with SYBR green chemistry (Applied Biosystems, USA) in two technical replicates. In brief, PCR reactions were carried out in 10 μl reaction volume containing 30 ng of first strand cDNA, 1X PCR buffer, 125 mM dNTPs, 1.5 mM MgCl2, 0.2 mM primers and 1U Taq polymerase and programed as follows: 50°C for 2 min and denaturation at 95°C for 10 min followed by 40 cycles of denaturation at 95°C for 15 sec and annealing and extension at 60°C for 1 min. The data obtained from different PCR runs or cDNA samples was analyzed using the mean of the CT values of the three biological replicates that were normalized to the mean CT values of the house keeping gene GAPDH. The expression values were calculated using the 2_ΔΔCt method, student's t-test was used to calculate significance and relative transcription levels are presented graphically [37].

Results

High throughput sequencing and assembly

A total of 642.93 million reads were generated from the transcriptome sequencing of eight samples with about 60.8 to 94.8 million reads per sample. The reads with adapter contamination and low base quality (≤ Q20) were removed using Trimmomatic and NGS-QCbox [29]. As a result, in total 593.94 million (92.38%) high quality reads were obtained. A total of 568.8 million reads (95.76%) were mapped on to chickpea reference genome (v1.0) [2] using TopHat2 software. On average, 71.1 million reads (95.7%) were mapped to the reference genome of chickpea for each sample (Table 1).

Table 1

Summary of Illumina sequencing data and mapped reads for the samples.

Genotype	Tissue	Treatment	Sample ID	Total reads	Mapped reads	Mapping rate (%)	Unique match	Multi-position match
Bivanij	Root	Control	BRDC	91,293,742	80,263,790	95.4	77,898,377	23,654,13
		Drought stress	BRDS	71,897,392	63,330,795	95.3	61,694,984	16,358,11
	Shoot	Control	BSDC	84,473,708	75,724,423	96.0	73,832,159	18,922,64
		Drought stress	BSDS	94,898,074	83,051,032	95.1	80,903,347	21,476,85
Hashem	Root	Control	HRDC	79,714,416	70,462,846	95.8	68,448,980	20,138,66
		Drought stress	HRDS	60,810,664	53,408,459	95.8	52,031,036	13,774,23
	Shoot	Control	HSDC	84,992,750	76,347,429	96.2	74,519,352	18,280,77
		Drought stress	HSDS	74,858,446	66,211,881	96.7	32,458,029	33,753,852
Total				642,939,192	568,800,655
Average				80,367,399	71,100,081.80	95.7

Global gene expression analysis

A gene was considered to be expressed in a sample if its FPKM ≥1. A total of 22,363 genes were expressed in at least one of the eight samples analyzed. Among these, 7,589 genes were found to be constitutively expressed in all samples with coefficient of variation less than 10%. The global expression profiles of all samples grouped root and shoot samples separately in accordance to their phenotypic drought tolerance/sensitive expression. In roots, two control and two drought stress samples were grouped together, whereas the correlation between the drought stress samples was less than the control samples in the shoots (Fig 1A). The eight samples were categorized based on highly expressed genes. The expression values of 80 highly expressed genes (with FPKM value >400) among the samples were used to develop the heatmap (Fig 1B). For better visualization of the heatmap, log10 transformed FPKM for the genes in all the samples were calculated. The 80 highly expressed genes were placed in two large clusters containing 25 and 55 members for individual cluster. The expression profiles were different between the clusters. Many sub-clusters were located per cluster based on the difference among the genotypes, the tissues and the treatments.

Fig 1

Correlation of the expression values of the samples.

(A) Dendogram depicting correlation among different samples based on global expression profiles. (B) Heatmap grouped the eight samples based on log10 transformed FPKM (>400) of the highly expressed genes. Color scale shows high and low expressions as red and green, respectively.

Differentially expressed genes (DEGs)

A total of 4,572 differentially expressed genes (DEGs) were identified by comparing the eight samples using Cuffdiff, of which 428 (9.37%) were novel. Of the 4,572, 94.1% (4,302) DEGs mapped uniquely and 5.9% (270) DEGs mapped to multiple locations in the genome. The highest number of DEGs (1,468) were identified between the roots and the shoots of Hashem under the control condition (HRDC vs HSDC) and the lowest number of DEGs (59) were in comparison of the roots of Bivanij and Hashem under the control condition (BRDC vs HRDC) (S1 Fig). For better understanding, differential gene expression analysis was performed based on the tissues, the treatments and the genotypes. Under the drought stress, DEGs in the roots of Bivanij were higher as compared to Hashem (891 vs 781) and were similar with DEGs in the shoots (507 vs 332). Overall, the number of down-regulated genes under drought stress conditions were higher as compared to the control conditions (Fig 2A). In addition, higher number of genes were differentially expressed between the tolerant and the sensitive genotypes in the drought stress as compared to the control condition. In both the control and the drought stress conditions, DEGs were more in the shoots than the roots. Under the drought stress, 127 and 52 genes were up-regulated in roots of Hashem and Bivanij, respectively. While, 155 and 121 genes were up-regulated in the shoots of Hashem and Bivanij, respectively (Fig 2B). Finally, DEGs between the tissues in two genotypes under both the control and the drought stress conditions were compared. As shown, the DEGs between the roots and the shoots in the sensitive genotype were more than those in tolerant genotype, and in both the genotypes, the DEGs were reduced during the drought stress as compared to the control condition. Overall, in all comparisons, number of the up-regulated genes were more than down-regulated genes (Fig 2C). In addition, overlapped differentially expressed genes between Bivanij and Hashem under the control and the drought stress conditions were analyzed (Fig 2D). In comparison the blue and the orange ovals, drought-responsive genes were identified, which were specifically expressed in the roots, including 171 (121 up-regulated in Hashem and 50 up-regulated in Bivanij), 51 (28 up-regulated in Hashem and 23 up-regulated in Bivanij) and eight drought-responsive, control and common genes, respectively. The other comparison was performed between green and red ovals in order to identify drought responsive genes, which were specifically expressed in the shoots. The results exhibited that 265 genes (153 up-regulated in Hashem and 112 up-regulated in Bivanij) were specifically expressed under the drought stress, 70 genes (28 up-regulated in Hashem and 42 up-regulated in Bivanij) were expressed under the control condition and 11 genes were common. Finally, 261, 169 and 17 drought responsive genes were expressed in the shoots, the roots and common genes between the tissues, respectively (S2 Table and S3 Table). Interestingly, only one gene (Ca_02523), an integral component of membrane and transporter activity was showing an overlap in all combinations.

Fig 2

Number of differentially expressed genes.

Comparison between: (A) The treatments (B) The genotypes and (C) The tissues. (D) Venn diagram shows the overlapped genes by comparing the genotypes under the control and the drought treatments in the roots and the shoots.

Functional classification of DEGs

Gene ontology (GO) analysis was performed to functionally categorize the 4,572 differentially expressed genes into three principal categories: cellular component, molecular function, and biological process (S2 Fig). In cellular component category, cell and cell part sub-categories were found to be mostly enriched and in molecular function category, binding, catalytic activity followed by transporter and transcription regulator were the most enriched categories. A large number of GO sub-categories, such as metabolic and cellular process, biological regulation, pigmentation, localization and response to stimulus were found to be significantly enriched in biological process category. GO analysis was also performed specifically for the drought responsive genes in the category of biological process (Fig 3). The sub-category, metabolic process with highest number of genes as well as oxidation-reduction were commonly found in all the conditions, while other biological processes were differently enriched. In the roots, several GO sub-categories associated with RNA transcription and cell wall were involved in Hashem, whereas localization, transport and carbohydrate metabolic processes were found in Bivanij. In addition to metabolic process and oxidation-reduction, primary metabolic process was highly enriched in shoots of Hashem, followed by catabolic, cellular and carbohydrate metabolic processes. Although a higher number of enriched sub-categories were in Hashem in both roots and shoots, however a number of major sub-categories, including response to stress, defense response and response to biotic stimulus were significantly and specifically enriched in the shoots of the tolerant genotype, indicating the vital role of Bivanij for plant drought tolerance.

Fig 3

GO enrichment analysis of DEGs identified between the tolerant and the sensitive genotypes under the drought stress in the roots and the shoots.

Differential expression of genes associated with KEGG metabolic pathways

KEGG pathway analysis was performed in order to compare the metabolic pathways between Bivanij and Hashem genotypes under the drought stress in the roots and the shoots (Fig 4). The enzymes with most frequency were monooxygenase (13.3%), aldolase (8%) and lactase (6%) in the shoots and adenylyl transferase (10%) and RNA polymerase (10%) in the roots. Totally, the results of KEGG pathway enrichment were in agreement with GO analysis. Various pathways of carbohydrate metabolism, biosynthesis of other secondary metabolites, energy metabolism, xenobiotics biodegradation and metabolism and nucleotide metabolism were enriched in the shoots with the highest number of the genes, whereas the pathways associated with carbohydrate metabolism were frequent in the roots.

Fig 4

KEGG significantly metabolic pathways in the roots and the shoots obtained from comparing the genotypes under the drought stress.

Differentially expressed transcription factors (TFs)

A total of 1,806 transcription factor related genes were identified from all the conditions and the top 20 TF families are shown in S3 Fig. The largest members of the TFs belonged to bHLH family (101), followed by ERF (87), kinase superfamily (76), NAC (74), MYB (72), WRKY (72), etc. The distribution of the differentially expressed TFs between Bivanij and Hashem under the drought stress was separately investigated in the roots and the shoots (Fig 5A–5E, S4 Table and S5 Table). Overall, 65 DEGs encoding TFs were found in the roots, including 42 (17 TF families) up-regulated in Hashem and 23 (16 TF families) up-regulated in Bivanij. In the shoots, 93 differentially expressed TFs were found between the genotypes during the stress condition, including 48 (24 TF families) up-regulated in Hashem and 45 (21 TF families) up-regulated in Bivanij. Six common differentially expressed TFs were also identified in the roots and the shoots. Interestingly, bHLH (Ca_09798) and MYB-related (Ca_19895) were induced in the roots of Hashem and repressed in the shoots. The distribution of the members of the differentially expressed TF families between the genotypes was different based on the tissue types and the genotypes. In the roots, the TF families, ERF (14%), MYB (12%), MYB-related (12%) and B3 (12%) were up-regulated in Hashem and HSF (17%), bHLH (13%), C3H (9%) and NF-YA (9%) were up-regulated in Bivanij with the highest members. However, in the shoots, the families MYB-related (11%), bHLH (9%), ERF (9%) and C2H2 (9%) with the high members were up-regulated in Hashem, whereas the members B3 (13%), NAC (11%) and MYB (9%) were up-regulated in Bivanij. Expression patterns of NAC transcription factor across different samples through heatmap which clearly indicated different expression level of NAC TF family in different tissues, treatments and genotypes.

Fig 5

Distribution of TF families in comparison between the genotypes under the drought stress.

(A) Up-regulated in the roots of Hashem (B) Up-regulated in the shoots of Hashem (C) Up-regulated in the roots of Bivanij and (D) Up-regulated in the shoots of Bivanij. (E) Heatmap showing expression profile of NAC transcription factor family across different samples.

DEGs coincident with drought tolerance-related QTLs

In the present study, we found two (Ca_04551 and Ca_04552) and three (Ca_04561, Ca_04564 and Ca_04569) DEGs located in “QTL-hotspot_a” (of 15 genes) and “QTL-hotspot_b” (of 11 genes) regions, respectively, reported earlier by Kale et al. [11] (S4 Fig, S6 Table). No significant difference observed in the expression patterns of these five DEGs between the tolerant and the sensitive genotypes in the roots, however Ca_04551 and Ca_04564 were up-regulated in the shoots of Hashem. Ca_04569 was up-regulated only in the roots of Bivanij under drought stress and was annotated as ‘inositol polyphosphate-related phosphatase’. Ca_04551 was annotated as, ‘aldo/keto reductase’ and ‘potassium channel, voltage-dependent, beta subunit, KCNAB-related’. Ca_04552 was involved in ‘chlorophyll A-B binding protein’. The functions of Ca_04561 were associated with ‘E3 ubiquitin-protein ligase Msl2, zinc RING finger’. Finally, Ca_04564 was found as ‘leucine-rich repeat’. Among them, Ca_04551, Ca_04564 and Ca_04569 were TFs related to bHLH, S1Fa and M-type families, respectively.

Validation of differentially expressed genes

In order to validate the sequencing results, quantitative real time PCR (qRT-PCR) analysis was performed using 12 drought stress responsive genes selected from the in silico analysis of the DEGs.

Comparison of the expression values of the RNA-Seq and the qRT-PCR results under the control and the drought stress conditions indicated a strong and consistent correlation (r2 = 0.80) between the RNA-Seq and the qRT-PCR methods (Fig 6A). For instance, Ca_04358 and Ca_04561 showed up-regulation under drought stress in both the genotypes studied. On the other hand, two genes Ca_04355 and Ca_04561 showed down-regulation under drought stress conditions. These results are support in silico data. Further, the expression values of the 12 genes for RNA-Seq and qRT-PCR are shown as graphical representation (Fig 6B). Most of them revealed similar expression patterns. Interestingly, the expression of Ca_15236 as an important gene involved in the drought stress was almost same (-2.99 fold change in RNA-Seq and -2.97 fold change in qRT-PCR) between the two approaches.

Fig 6

Validation of the RNA-Seq using the qRT-PCR.

(A) Close correlation (r2 = 0.80) of 12 selected genes between two approaches. (B) Comparison of the RNA-Seq and the qRT-PCR considering the expression of the 12 selected genes in some combinations for instance.

Discussion

Occurrence of drought during flowering time (early drought) and late at pod filling stage (terminal or end season drought) were reported to have significant impacts on chickpea yields. Drought being a complex trait, therefore understanding genetics of trait and genes associated with the drought tolerance has been the major area of focus in case of chickpea. The knowledge about these target drought responsive genes would enable development of drought tolerant cultivars through molecular breeding or transgenic approaches. Genomic regions/QTLs for drought tolerance have been reported using bi-parental mapping populations as well as diverse germplasm lines [38,8]. Attempts were also made to understand drought stress transcriptome [39,18,21,22], however these studies were on desi genotypes. In present study, drought tolerance was investigated in two Iranian kabuli chickpea genotypes Bivanij (tolerant) and Hashem (susceptible) at early flowering stage under drought stress conditions using RNA-Seq approach to understand the global transcriptome.

The GO enrichment analysis could functionally classify the DEGs and about 10% of the whole DEGs were assigned in response to stimulus. The most important stress-related GO terms such as response to stress, defense response and response to biotic stimulus were found to be specifically involved for up-regulated genes in the shoots of the tolerant genotype as compared to the sensitive genotype. Generally, the GO analysis suggested that the DEGs in the shoots were responsible for creating difference in drought tolerance between the genotypes. Transcriptome studies conducted in chickpea desi type have also noted the same GO terms in response to abiotic stresses [18,21,22].

Plant hormones regulate growth and development and are involved in various environmental stresses through signal transduction. Abscisic acid (ABA) is a key hormone in plants and its roles, such as seed dormancy and germination, modulation of roots architecture, senescence, stomata regulation and abiotic stress responses have been reviewed by Sah et al. [40]. ABA can also induce transcription factors and subsequent stress responsive genes, which confer tolerance to abiotic stresses, such as drought [41,42]. In this study, the gene (Ca_00449) involved in carotenoid biosynthesis pathway in shoots, which produces ABA was identified. This gene annotated as TF (FAR1 family) was up-regulated in Bivanij. The enzyme zeta-carotene desaturase (EC:1.3.5.6) was controlled by this gene and catalyzes the conversion of neurosporene to lycopene. Several ABA biosynthesis-deficient mutants have been shown to be sensitive to environmental stresses [43]. Overall, manipulation of ABA biosynthesis genes are important tools to enhance the tolerance to stresses. In addition, a number of ABA and ethylene responsive genes were differentially expressed between the two genotypes. Auxin and cytokinin hormone biosynthesis pathways were also identified to show differential regulation between the genotypes under the drought stress in the shoots and the roots, respectively. Both of them were up-regulated in Bivanij, indicating the successes of this genotype for the growth and the development during the drought stress.

Proline is a primary metabolite and accumulates as an osmotic adjustment during dehydration [44]. Further, proline protects macromolecules and enzymes under stress conditions and provide energy for growth and survival [45]. Arginine and proline metabolism is one of the proline biosynthesis pathways, which was identified in the present study. Remarkably, two enzymes were identified in proline biosynthesis, including glutamate 5-kinase (EC:2.7.2.11) and glutamate-5-semialdehyde dehydrogenase (EC:1.2.1.41) induced by a B3 transcription factor. This gene (Ca_24241) was up-regulated in the shoots of Bivanij as compared to Hashem under the drought stress treatment.

The genes of three pathways in energy metabolism, including photosynthesis, oxidative phosphorylation and methane metabolism were up-regulated in the shoots of Hashem. Considerably, this can be an adaptive response to reduce energy consumption in Bivanij under the drought condition. Increased levels of carbon fixation and photosynthesis due to open stomata led to reduction of osmotic pressure and thereby increasing the rate of transpiration in Hashem (as shown in Sade et al. [46]). Reduction of oxidative phosphorylation and photosynthesis have been reported in Bombax ceiba due to increased ABA, stomatal closure and increased level of proline under drought stress [47]. Soluble sugars, such as starch and sucrose are known to regulate osmotic pressure of plants [48,49]. In this study, the tolerant and the sensitive genotypes exhibited different responses in the shoots and the roots during starch and sucrose metabolism. In the shoots, two genes were up-regulated in Hashem and one gene was up-regulated in Bivanij, whereas in the roots, three genes were induced in Bivanij and only one gene was induced in Hashem. Increased amount of soluble sugars can be involved in growth and development or osmotic adjustment and enhancing stress tolerance [50,51].

Drought stress damages the photosynthetic pigments and results in reduction of chlorophyll content [52]. The KEGG pathway of porphyrin and chlorophyll metabolism conducts the biosynthesis of chlorophyll a and b. In the present study, there was no significant difference between two genotypes for the genes in that pathway, while the comparison of the tissues under drought stress identified several genes up-regulated in the shoots of both genotypes. It seems that increased degradation of chlorophyll content in the shoots as compared to the roots led to enhanced expression of the genes encoding the enzymes of chlorophyll biosynthesis pathway under drought stress.

The genes encoding biosynthesis pathways of secondary metabolites were mostly induced in Bivanij except the phenylpropanoid biosynthesis, which was also induced in Hashem. Flavonoids are a large group of secondary metabolites, which act as antioxidant and protect plants against various biotic and abiotic stresses [53]. In this study, shikimate O-hydroxy cinnamoyl transferase enzyme (EC:2.3.1.133) involved in flavonoid biosynthesis was induced by ARF transcription factor (auxin response factor). This gene (Ca_05702) was up-regulated in shoots of Bivanij as compared to Hashem under the drought stress. Overexpression of gene AtMYB12 involved in flavonoid biosynthesis enhanced drought and salt tolerance through accumulation of flavonoids in transgenic Arabidopsis plants [54].

Reactive oxygen species (ROS) are generated in cells during stress conditions, and enzymes like catalase, superoxide dismutase and peroxidase play an important role in detoxifying these compounds [55]. Four genes involved in hydrogen peroxide catabolic process were identified in present study, three of them were up-regulated in Hashem and only one gene (Ca_04125) was up-regulated in Bivanij. These results indicated the complexity of tolerance to abiotic stresses, such as drought.

Plants respond to diverse stresses through perceiving the stresses by receptors or sensors at membrane level. Then the secondary messages, such as calcium ion, inositol phosphates, ROS, cyclic nucleotides and ABA are involved to transduce the signals into nucleus. Finally, the signals induce or repress various transcription factors (TFs) and consequently, many stress-related downstream genes are regulated through specific binding the TFs with their corresponding cis-elements in promoters [56,57]. In the present study, various TF families were identified, and their distribution and frequency were found to be greatly different, depending on the tissues, the treatments and the genotypes. 21 genes containing stress-related annotations were identified by comparing Hashem and Bivanij under the drought stress. Four and 17 differentially expressed genes were found in the roots and the shoots, respectively, suggesting the higher roles of the shoots to make different responses to drought stress between Hashem and Bivanij. Remarkably, 14 genes were up-regulated in Bivanij and only seven genes were up-regulated in Hashem, indicating the potential of Bivanij to be considered as the tolerant genotype (S7 Table). Moreover, the stress-related genes were involved in various stresses, suggesting a crosstalk among different stresses. Recently Nakashima et al. [58] confirmed presence of interactions between TFs and regulatory networks involved in diverse stress responses, including drought, cold and heat. Out of 21 genes with stress-related annotations, nine genes (three in the roots and six in the shoots) were considered as TFs. The TFs, C3H (2), NAC (2) and MYB-related (1) were up-regulated in Bivanij under the drought stress, whereas M-type (2), MYB-related (1) and G2-like (1) were up-regulated in Hashem.

The gene families such as AP2/ERF and heat shock protein (HSP) are crucial protein factors involved in many developmental processes and stress responses. HSP90 as a highly conserved chaperone protein assists to refold, stabilize and maintain other proteins, which were damaged under stress conditions [59,60]. Agarwal et al. [61] dissected AP2/ERF and HSP90 genes families in chickpea, pigeonpea and common bean under various biotic and abiotic stresses. Interestingly, four genes from AP2/ERF (Ca_08430, Ca_09578, Ca_23799 and Ca_15031) and two genes from HSP90 (Ca_25811and Ca_25602) were found commonly with our differentially expressed genes; however, none of them were differentially expressed between Hashem and Bivanij.

Considering the vital roles of TFs in environmental stresses, several key TFs, such as AP2/ERF [59], AREB/ABF [62], MYB [63], WRKY [64], NAC [65] and bZIP [66] have been reported to improve drought tolerance in different plants through genetic engineering approach. The TF, CIPK25 isolated from chickpea enhanced the root growth and improved drought and salinity stresses in tobacco [67]. Several target stress responsive genes, such as RD29A, RD10, COR15A, COR47, KIN1 and DREB2A have been reported to be induced in Arabidopsis due to the transgene CarNAC4 isolated from chickpea [68]. The expression levels of several CaNAC genes have been assessed between two chickpea contrasting genotypes under dehydration stress and the TF CaNAC16 was introduced as potential candidate gene, contributing to the better drought tolerance of Hashem and ILC482 cultivars [69]. Noteworthy, we found the gene CaNAC16 (Ca_18090) up-regulated in Bivanijas compared to Hashem under the drought stress, indicating the vital role of this TF during plant drought tolerance.

Remarkably, five DEGs identified in this study were found in two drought tolerance-related “QTL-hotspot” regions in chickpea. The Ca_04561 gene encoding E3 ubiquitin-protein ligase, has been reported to be involved in the increasing of proline contents and enhancing of drought tolerance in Arabidopsis [70]. Finally, expression profiling through qRT-PCR of 12 selected genes supported the results obtained through the RNA-Seq analysis.

Conclusion

In summary, the transcriptome analysis revealed complex genes and networks involved in chickpea drought tolerance. The gene expression values were different across the tissues, followed by the treatments and the genotypes. Although, a high number of genes were differentially expressed in the roots of both genotypes under the stress condition, however the comparison between the genotypes exhibited the important role of the shoots for conferring higher tolerance in Bivanij. In other words, the genotypes employed same genes and mechanisms in the roots as compared to the shoots. GO and KEGG pathway analysis identified several molecular mechanisms involved in the stress response and their corresponding drought-related pathways. The role of key transcription factors in the drought tolerance mechanism has been deciphered and these genes will serve as useful targets for both future research and breeding for drought tolerance/abiotic stresses in chickpea.

List of primers used in the experiment.

(DOC)

Click here for additional data file.

Drought responsive genes in the root tissues.

(XLS)

Click here for additional data file.

Drought responsive genes in the shoot tissues.

(XLS)

Click here for additional data file.

Drought responsive genes encoding TFs in the root tissues.

(XLS)

Click here for additional data file.

Drought responsive genes encoding TFs in the shoot tissues.

(XLS)

Click here for additional data file.

DEGs in “QTL-hotspot” reported in Kale et al. 2015.

(DOC)

Click here for additional data file.

Stress-related annotations for differentially expressed genes between the genotypes under drought stress.

The genes were also annotated for some important characteristics involved in various stresses.

(DOC)

Click here for additional data file.

Number of up and down-regulated genes for 12 comparative combinations of the samples.

(DOC)

Click here for additional data file.

GO analysis of 4572 differentially expressed genes in the experiment.

(DOC)

Click here for additional data file.

Distribution of 20 top TF families identified in 4572 DEGs.

(DOC)

Click here for additional data file.

Differentially expressed genes identified in the “QTL-hotspot_a” and “QTL-hotspot_b” in the present study.

The genes corresponding to the present DEGs are indicated with the red arrows.

(DOC)

Click here for additional data file.