Comparative Root Transcriptomics Provide Insights into Drought Adaptation Strategies in Chickpea (Cicer arietinum L.)

Drought adversely affects crop production across the globe. The root system immensely contributes to water management and the adaptability of plants to drought stress. In this study, drought-induced phenotypic and transcriptomic responses of two contrasting chickpea (Cicer arietinum L.) genotypes were compared at the vegetative, reproductive transition, and reproductive stages. At the vegetative stage, drought-tolerant genotype maintained higher root biomass, length, and surface area under drought stress as compared to sensitive genotype. However, at the reproductive stage, root length and surface area of tolerant genotype was lower but displayed higher root diameter than sensitive genotype. The shoot biomass of tolerant genotype was overall higher than the sensitive genotype under drought stress. RNA-seq analysis identified genotype- and developmental-stage specific differentially expressed genes (DEGs) in response to drought stress. At the vegetative stage, a total of 2161 and 1873 DEGs, and at reproductive stage 4109 and 3772 DEGs, were identified in the tolerant and sensitive genotypes, respectively. Gene ontology (GO) analysis revealed enrichment of biological categories related to cellular process, metabolic process, response to stimulus, response to abiotic stress, and response to hormones. Interestingly, the expression of stress-responsive transcription factors, kinases, ROS signaling and scavenging, transporters, root nodulation, and oxylipin biosynthesis genes were robustly upregulated in the tolerant genotype, possibly contributing to drought adaptation. Furthermore, activation/repression of hormone signaling and biosynthesis genes was observed. Overall, this study sheds new insights on drought tolerance mechanisms operating in roots with broader implications for chickpea improvement.


Introduction
Drought imposes a serious threat to crop growth and productivity [1]. Moreover, the intensity and frequency of drought periods are predicted to increase in the future climate [2]. Feeding the world's growing population under water limiting conditions will be a big challenge. To sustain agricultural patterns related to different functional categories were observed. We discuss the possible role of these genes in drought tolerance that can facilitate the development of drought-tolerant chickpea varieties in the future.

Phenotypic Responses
When exposed to drought stress, root-related traits such as dry weight, length, and surface area of the tolerant chickpea genotype (ICC8261) were significantly higher than the sensitive genotype (ICC283) at vegetative stage (VS) ( Figure 1A-C). In contrast, at the reproductive stage (RS), root length and surface area of the sensitive genotype were higher than the tolerant genotype ( Figure 1B,C). The average root diameter of the tolerant genotype was higher than the sensitive genotype at RS ( Figure 1D). Root volume and relative water content (RWC) were similar in both genotypes at all growth stages ( Figure 1E,F).
Int. J. Mol. Sci. 2020, 21, x FOR PEER REVIEW 3 of 20 genes in drought tolerance that can facilitate the development of drought-tolerant chickpea varieties in the future.

Phenotypic Responses
When exposed to drought stress, root-related traits such as dry weight, length, and surface area of the tolerant chickpea genotype (ICC8261) were significantly higher than the sensitive genotype (ICC283) at vegetative stage (VS) ( Figure 1A-C). In contrast, at the reproductive stage (RS), root length and surface area of the sensitive genotype were higher than the tolerant genotype ( Figure  1B,C). The average root diameter of the tolerant genotype was higher than the sensitive genotype at RS ( Figure 1D). Root volume and relative water content (RWC) were similar in both genotypes at all growth stages ( Figure 1E,F).  (SE). Physiological data were obtained from three independent biological replicates. Statistically significant differences between sensitive and tolerant genotypes obtained by one-way ANOVA at p < 0.05 are depicted by an asterisk (*). Furthermore, shoot dry biomass was significantly higher in the tolerant genotype under drought stress at all developmental stages ( Figure 1G). The chlorophyll content was significantly higher in the tolerant genotype during drought stress at RTS ( Figure 1H). The specific leaf area (SLA) was lower in the drought-tolerant genotype at RS ( Figure 1I). Overall, chickpea genotypes ICC8261 and ICC283 showed contrasting responses to drought stress, which is consistent with the previous findings suggesting that ICC8261 shows features of drought tolerance [14,16].

Transcriptional Responses and GO Analysis
To gain mechanistic insights into regulatory mechanisms underlying the drought stress response, we did RNA sequencing. The tolerant genotype showed a total of 2161 differentially expressed genes (DEGs) in response to drought stress at the vegetative stage, out of which 1214 DEGs Error bars represent standard errors (SE). Physiological data were obtained from three independent biological replicates. Statistically significant differences between sensitive and tolerant genotypes obtained by one-way ANOVA at p < 0.05 are depicted by an asterisk (*). Furthermore, shoot dry biomass was significantly higher in the tolerant genotype under drought stress at all developmental stages ( Figure 1G). The chlorophyll content was significantly higher in the tolerant genotype during drought stress at RTS ( Figure 1H). The specific leaf area (SLA) was lower in the drought-tolerant genotype at RS ( Figure 1I). Overall, chickpea genotypes ICC8261 and ICC283 showed contrasting responses to drought stress, which is consistent with the previous findings suggesting that ICC8261 shows features of drought tolerance [14,16].

Transcriptional Responses and GO Analysis
To gain mechanistic insights into regulatory mechanisms underlying the drought stress response, we did RNA sequencing. The tolerant genotype showed a total of 2161 differentially expressed genes (DEGs) in response to drought stress at the vegetative stage, out of which 1214 DEGs were upregulated  Figure 2A). Similarly, in the sensitive genotype, 1873 DEGs were obtained, out of which 821 DEGs were upregulated and 1052 DEGs were downregulated (Figure 2A). Venn diagram analysis showed that 261 DEGs were commonly upregulated and 262 DEGs were commonly downregulated in these genotypes ( Figure 2B) were upregulated and 947 DEGs were downregulated (Figure 2A). Similarly, in the sensitive genotype, 1873 DEGs were obtained, out of which 821 DEGs were upregulated and 1052 DEGs were downregulated ( Figure 2A). Venn diagram analysis showed that 261 DEGs were commonly upregulated and 262 DEGs were commonly downregulated in these genotypes ( Figure 2B). Interestingly, 61 DEGs showing upregulation in the tolerant genotype were downregulated in the sensitive genotype, and 20 DEGs showing downregulation in the tolerant genotype were upregulated in the sensitive genotype. At the reproductive stage, a total of 4109 DEGs were identified in the tolerant genotype in response to drought stress, out of which 2139 DEGs were upregulated and 1970 DEGs were downregulated ( Figure 2C). In the sensitive genotype, 3772 DEGs were identified wherein 2038 were up-and 1734 were downregulated. Venn diagram analysis showed that 1163 and 1059 DEGs were commonly up-and downregulated in both the genotypes, respectively ( Figure 2D). Out of these, 19 DEGs that were upregulated in the tolerant genotype were downregulated in the sensitive genotype, while 3 DEGs that showed downregulation in the tolerant genotype were upregulated in the sensitive genotype ( Figure 2D). The reliability of RNA-Seq gene expression was confirmed by qRT-PCR. The fold changes obtained from RNA-seq and qRT-PCR showed good correlation with each other, thus validating the RNA-seq data ( Figure S1).
Next, the functional over-representation of drought-responsive genes was performed by gene ontology (GO) analysis. GO terms classified genes into biological processes, molecular function, and cellular components (Figure 3). The most dominant terms recurring in the biological process category at vegetative and reproductive stages were 'cellular process', 'metabolic process', 'response to stimulus', 'response to abiotic stress', and 'response to hormones'. In the molecular function category, 'catalytic activity', 'binding', 'hydrolase', 'transporter activity', and 'antioxidant activity' terms were At the reproductive stage, a total of 4109 DEGs were identified in the tolerant genotype in response to drought stress, out of which 2139 DEGs were upregulated and 1970 DEGs were downregulated ( Figure 2C). In the sensitive genotype, 3772 DEGs were identified wherein 2038 were up-and 1734 were downregulated. Venn diagram analysis showed that 1163 and 1059 DEGs were commonly up-and downregulated in both the genotypes, respectively ( Figure 2D). Out of these, 19 DEGs that were upregulated in the tolerant genotype were downregulated in the sensitive genotype, while 3 DEGs that showed downregulation in the tolerant genotype were upregulated in the sensitive genotype ( Figure 2D). The reliability of RNA-Seq gene expression was confirmed by qRT-PCR. The fold changes obtained from RNA-seq and qRT-PCR showed good correlation with each other, thus validating the RNA-seq data ( Figure S1).
Next, the functional over-representation of drought-responsive genes was performed by gene ontology (GO) analysis. GO terms classified genes into biological processes, molecular function, and cellular components ( Figure 3). The most dominant terms recurring in the biological process category at vegetative and reproductive stages were 'cellular process', 'metabolic process', 'response to stimulus', 'response to abiotic stress', and 'response to hormones'. In the molecular function category, 'catalytic activity', 'binding', 'hydrolase', 'transporter activity', and 'antioxidant activity' terms were present. The cellular component category contained terms such as 'cell part', 'cytoplasm', 'intracellular', 'chloroplast', and 'cell wall'. Overall, these results suggest that roots underwent global transcriptional reprogramming during drought stress.

Root-Nodule Development Genes
Chickpea roots contain nodules that help in the establishment of a symbiotic association with nitrogen-fixing rhizobia. The expression of putative nodulation receptor kinase (LOC101507037) and calcium/calmodulin-dependent serine/threonine-protein kinase DMI3 (LOC101513751) were downregulated in the sensitive genotype both at VS and RS; however, the tolerant genotype showed marginal downregulation at RS only ( Figure 6, Table S3). In contrast, receptor-like kinases 3 (LOC101496137) and CLAVATA1 (LOC101488348) were upregulated at VS and RS in the tolerant genotype ( Figure 6).

Root-Nodule Development Genes
Chickpea roots contain nodules that help in the establishment of a symbiotic association with nitrogen-fixing rhizobia. The expression of putative nodulation receptor kinase (LOC101507037) and calcium/calmodulin-dependent serine/threonine-protein kinase DMI3 (LOC101513751) were downregulated in the sensitive genotype both at VS and RS; however, the tolerant genotype showed marginal downregulation at RS only ( Figure 6, Table S3). In contrast, receptor-like kinases 3 (LOC101496137) and CLAVATA1 (LOC101488348) were upregulated at VS and RS in the tolerant genotype ( Figure 6).

Discussion
Drought imposes a serious threat to agriculture. The root system plays a crucial role in the uptake of water and nutrients from the soil, thus helps in coping with drought stress. Root development and architecture are determined by intrinsic genetic properties and/or modulated by various environmental factors such as water availability [32,33]. So, we compared phenotypic and root-specific transcriptional responses of two contrasting chickpea genotypes at vegetative (VS), reproductive transition (RTS), and reproductive (RS) stages.
The tolerant genotype had higher root and shoot dry biomass under drought stress, possibly contributed by better water foraging capacity of roots and/or onset of resilience mechanisms. Indeed root length and surface area of tolerant genotype were higher at the vegetative stage; however, conservative growth was observed at the reproductive stage. In contrast, the root diameter showed the opposite trend. Longer roots help plants to forage water from deeper soils, and the larger surface area provides better opportunities for mycorrhizal colonization that facilitate nutrient acquisition [34]. The opposite response of root diameter could be explained by the fact that diameter controls the length and surface area for given biomass allocated to the root system, thus summarize the overall effect [35]. Chlorophyll content of the drought tolerant genotype was higher under drought stress at RTS, suggesting improved light-harvesting and energy production capacity to support growth and reproduction [36]. Also, lower specific leaf area in the tolerant genotype can improve water-use efficiency, as reported previously [37,38].
RNA sequencing allows genome-scale quantification of transcriptome changes underlying developmental transitions and stress responses [39]. In the present work, genes related to various biological functions were altered in response to drought stress. Transcription factor (TF) genes related to ABA-dependent (ABI5) and ABA-independent (DREB1A and DREB1C) pathways were upregulated in the tolerant genotype at RS, indicating robust induction of multiple regulatory networks to impart drought tolerance. Also, genotype-specific alterations in other TF genes such as AP2/ERF, bHLH, bZIP, DOF, and WRKY were observed. The AP2/ERF forms a large group of plant-specific TFs induced by multiple stresses and phytohormones [40,41]. ERFs were previously shown to be induced by wounding and biotic stresses [42]. Likewise, WRKY family TFs have been suggested to play a role in adaptation to biotic and abiotic stresses, and modulate stress-responsive signaling and ROS production [43]. WRKY33 is involved in the activation of peroxidases and glutathione S-transferases under drought stress, playing a role in scavenging ROS [44]. We found that bHLH family TF genes (VIP1, bHLH18-like, and bHLH93-like) were upregulated in the tolerant genotype. These TFs may contribute to improved defense; for instance, VIP1 activates stress-responsive genes via the mitogen-activated protein kinase pathway [45]. Also, other TF genes such as C2C2-Dof, MYB, BTB/POZ, and TAZ were upregulated in the tolerant genotype. Overall, upregulation of TFs in the tolerant chickpea genotype may activate stress-responsive genes conferring drought tolerance.
Protein kinases are integral components of signal transduction pathways playing an important role in the perception and activation of stress-responsive pathways [46]. The MKK1 gene was upregulated in tolerant genotype both at VS and RS upon drought exposure. It may contribute to the maintenance of cellular homeostasis by lowering ROS levels under stress [47,48]. CDPK family members (CDPK-SK5, CDPK3, and CDPK4) were upregulated in both the genotypes. CDPKs are activated by binding of intracellular Ca 2+ with a calmodulin-like domain that eventually regulates the downstream targets [49]. For instance, CDPK4 phosphorylates NADPH oxidase to regulate ROS production during stress [50][51][52]. We found upregulation of RLKs (IMK2, CRK25, and PERK4) in the tolerant genotype at VS, while PERK14 and HSL2 showed up-regulation in both genotypes at RS. RLKs get rapidly activated during the early stages of drought stress, as reported previously [53,54]. PERK4 modulates root tip growth via ABA signalling, and HSL2 regulates lateral root emergence by controlling cell separation [53,54]. The WAK-like 20 was upregulated in the tolerant genotype at VS, which possibly facilitate communication between the cell wall and cytoplasm, and regulate stress responses [55][56][57][58][59]. Therefore, their enhanced expression in the tolerant genotype may contribute to cell wall loosening that minimizes root growth inhibition under drought. Also, CIPK1-like was upregulated in the tolerant genotype at RS. The CIPK1 functions as a regulator of ABA-dependent and independent signaling pathways by interacting with calcium sensors [60]. Taken together, early and robust upregulation of protein kinases possibly contributes to drought tolerance by activation of downstream stress signaling pathways.
ROS act as both signaling and damaging molecules whose levels are determined by antioxidant molecules and enzymes [47]. The expression of many antioxidant enzymes, such as glutathione transferases (GST) and peroxidases (PER), was upregulated in the tolerant genotype. Simultaneous upregulation of GST and PER indicate that multiple ROS detoxification components efficiently scavenge ROS under drought; however, their downregulation in the sensitive genotype may lead to excessive ROS accumulation, causing root growth inhibition. The expression of RBOHH and RBOHD were downregulated in the tolerant genotype at VS, whereas RBOHE and RBOHA were upregulated at RS. Previously, downregulation of RBOH family genes has been reported under drought stress [61,62]. The down-regulation of RBOHD and RBOHH may cause inhibition of ABA-dependent ROS signaling, which allows marginal root growth inhibition in the tolerant genotype. Thus, the synergistic action of a battery of antioxidant enzymes confers stress tolerance in plants [63].
Stress-induced signaling molecules are transported via active or passive transporters. Active transport is mediated by ATP dependent transporters such as ATP-binding cassette (ABC), whereas passive transport occurs through ion-channels and carriers such as aquaporins [64][65][66][67].
ABC transporters were upregulated in the tolerant genotype at VS (ABCB15-like, ABCB19, and ABCC12-like) and RS (ABCB15-like, ABCB29, and ABCG22-like). Upregulation of ABC transporters can be linked to multiple functions, which include auxin transport, suberin formation, ABA signaling and transport, pollen development, leaf water retention, and stress tolerance [64][65][66][67]. Furthermore, sugar (N3, SWEET1-like, SWEET4) and nitrate (NRT1/PTR FAMILY) transporters were upregulated in both genotypes. The SWEET11 facilitates sugar supply to roots in a cell-type-specific manner and participates in the establishment of symbiotic association [68]. Furthermore, SWEET12 and SUC2 functions are associated with the reallocation of carbohydrates to roots for efficient root development under drought [69]. Also, NRT transporters are involved in drought stress responses by controlling ABA transport [70]. Aquaporins passively transport water molecules, and we found upregulation of PIP2-1-like and TIP2-2 in the tolerant genotype. Consistent with this, overexpression of OsPIP1-3 and SlTIP2-2 led to the maintenance of leaf water potential under drought stress in rice and tomato, respectively [71,72]. Also, ABA regulates aquaporin activity and hydraulic conductivity to maintain favorable plant water status [73]. Therefore, aquaporins could serve as ideal targets for engineering drought tolerance in chickpea.
Leguminous plants harbor nitrogen-fixing rhizobia in root-nodules, which facilitates the uptake of nutrients and water. It has been observed that Medicago with nodulated roots shows delayed leaf senescence as compared to non-nodulated plants under drought stress [74]. We also found downregulation of nodulation receptor kinase and calcium-dependent kinase DMI3 in the sensitive genotype. DMI3 acts upstream of nodulation-signaling pathway proteins (NSP1/NSP2) and is a regulator of Nod-factor genes [75,76]. Thus, their down-regulation in the sensitive genotype suggests a decrease in root nodulation, resulting in low nitrogen assimilation. In contrast, LYK3 and CLAVATA1 were upregulated in the tolerant genotype. LYK3 regulates rhizobial infection via Nod-factor genes [77], and CLAVATA1 regulates root length and nodulation [78]. Therefore, their upregulation may contribute to root nodulation and nitrogen assimilation under drought stress.
Cross-regulatory hormonal network modulates plant growth and stress response in plants [79]. Auxin regulates root and shoot growth under optimal and stressful conditions [80]. In this study, auxin biosynthesis gene (YUCCA2) was downregulated in the sensitive genotype at VS, indicating attenuation of auxin biosynthesis. Consistent with this inhibition of root biomass, length, and surface area were observed in the sensitive genotype. In contrast, PIN genes (PIN1b, PIN1c-like, and PIN4) were downregulated in the tolerant genotype at RS. PIN proteins act as auxin efflux carriers for basipetal auxin transport, which in turn controls lateral root formation [81]. Their downregulation may be associated with the reduction of root length and root surface area. Likewise, cytokinin controls lateral root, apical meristem, and vascular system development [82,83]. Cytokinin biosynthesis (IPT5) and signaling (AHK5) genes were downregulated in the tolerant genotype. It has been shown that AHK5 inhibits root elongation via ETR1-dependent ABA and ethylene signaling pathways [84]. Also, earlier transcriptome studies highlight the dual role of ethylene to regulate growth and defense responses during stress [85]. We found upregulation of ethylene-responsive transcription factors (ERF) such as ERF1-like, ERF113-like, ERF13-like, and ERF2-like in the tolerant genotype. ERFs integrate ethylene and jasmonic acid signaling pathways [86] and participate in stress responses [87,88]. Ethylene is known to control lateral root formation via modulating auxin biosynthesis and transport [89], and root and shoot growth via cross-talk with ABA which is a central mediator of drought stress-induced signaling [90]. The expression of ABA hydroxylases (CYP707A3-like and CYP707A1-like) was upregulated in the tolerant genotype, which could catalyze oxidative degradation of ABA, maintaining optimal ABA levels during drought [91]. The PYL genes were downregulated in tolerant genotype, while phosphatases (PP2C) were upregulated in the sensitive genotype. ABA binds with PYL/PYR/RCAR receptors to inhibit the activity of PP2Cs that negatively regulate ABA signaling through repression of SnRK2. PP2C inhibition leads to SnRK2 de-repression, eventually phosphorylating and activating downstream transcription factors [92]. The PYL4 downregulation and CYP707A3 activation in drought-tolerant chickpea genotypes suggest inhibition of ABA signaling at VS, whereas the up-regulation of NCED1 and PYL4 indicate activation of ABA signaling at RS. Also, ABA accumulation is known to inhibit root growth by promoting ethylene biosynthesis [93]. Gibberellic acid induces cell proliferation and elongation growth in plants [94,95]. DELLA proteins are negative regulators of GA signaling that act downstream of GA receptors. The binding of DELLA to GA receptor GID1 leads to their degradation and activation of GA function [96]. We observed downregulation GA receptor (GID1B-like) and DELLA genes in both the genotypes. This indicates that repression of GID1 and DELLA might result in reduced GA responses, thus attenuating root growth under drought stress. The role of lipid-derived phytohormone Jasmonate in abiotic stress tolerance is emerging [97]. We found that expression of JA biosynthesis genes such as LOX, LOX homologs, and TIFY/JAZ were significantly upregulated in the tolerant genotype, suggesting enhanced production of JA in this genotype. Previously, it was reported that activation of allene oxide synthase (AOS), an enzyme involved in JA biosynthesis, enhances drought tolerance in chickpea [97]. In addition, JA and ABA pathways are involved in the regulation of stomatal opening and closing to improve drought tolerance [98]. Overall, complex hormonal-crosstalk controls root growth under drought stress in chickpea.

Plant Material and Experimental Set-Up
The seeds of two chickpea (Cicer arietinum L.) cultivars viz. ICC283 (Desi type) and ICC8261 (Kabuli type), exhibiting contrasting responses to drought stress, were procured from ICRISAT (International Crops Research Institute for the Semi-Arid Tropics) Hyderabad, India. The genotype ICC283 is drought-sensitive, while ICC8261 is a drought-tolerant genotype. This phenotypic divergence in response to drought is attributed to their root system architecture, where the susceptible ICC283 genotype possesses shallow root systems, while the tolerant ICC8261 genotype possesses prolific root systems [14,16]. The seeds of these genotypes were sown in 27-cm polypropylene pots containing 9.5 kg of soil. Pots were placed in controlled growth conditions of the day (28 • C) and night (23 • C) temperatures, and randomly assigned as well-watered (WW) and drought-stressed (DS). The optimal water level of 95% ASWF (available soil water fraction) was maintained in WW pots by irrigating pots on alternate days. DS pots were maintained at 95% ASWF at sowing and were brought to 70% ASWF by dry-down fifteen days before the sampling. The pots were weighed every day and water was added to compensate the water lost through transpiration. The experiment was performed in a 3 × 2 × 2 (three time points, two genotypes, and two treatment conditions) completely randomized block design. Root tissue was harvested for physiological and RNA-seq analyses from WW and DS treatments at 30, 50, and 70 days after sowing (DAS). These time points reflect the vegetative stage (VS; 30 DAS), the reproductive transition stage (RTS; 50 DAS), and the reproductive stage (RS; 70 DAS). Root tissue harvested for RNA-seq was immediately put in liquid nitrogen. Three biological replicates were used to carry out physiological and transcriptional analyses.

Phenotypic Analysis
Phenotypic root and shoot traits were quantified in the sensitive and tolerant chickpea genotypes at three developmental stages viz. VS (30 DAS), RTS (50 DAS), and RS (70 DAS). Dry root and shoot biomass were determined after drying plant material at 65 • C for 3 days in the oven. To determine root-related traits such as length, surface area, diameter and volume harvested roots were washed with water and dried between the folds of filter paper. Then images were captured to compute these parameters. Relative water content (RWC) was calculated as RWC = [(fresh weight − dry weight)/(turgid weight − dry weight)] × 100. Leaf-related traits such as chlorophyll content (SPAD chlorophyll meter) and specific leaf area (SLA) was calculated as SLA = leaf area (cm 2 )/leaf dry weight (g). All the phenotypic analyses were performed in at least three replications, and data were reported as mean and standard error (SE). Means were statistically compared by one-way ANOVA at the probability level of 5% (p < 0.05) by using SPSS software (IBM, New York, NY, USA).

RNA Extraction and Library Construction
Total RNA was extracted with an RNeasy Mini Kit (Qiagen, Hilden, Germany). RNA quantification was done by using a Nanodrop Lite spectrophotometer (Thermo Fisher Scientific, Wilmington, MA, USA) and the integrity of RNA was checked with Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). For library preparation and sequencing, RNA samples with a 260/280 absorbance ratio between 2.1 to 2.2, 260/230 absorbance ratio between 2.0 to 2.5, and RIN (RNA integrity number) of more than 7.0 were used.
The enrichment of mRNA samples was done using a MicroPoly(A)Purist Kit (Thermo Fisher Scientific, Waltham, MA, USA) by following the guidelines of the manufacturer. The RNA-seq libraries were prepared using Ion Total RNA-Seq Kit v2 (Thermo Fisher Scientific, Waltham, MA, USA) and sequenced by using an Ion Proton System (Thermo Fisher Scientific, Waltham, MA, USA). The low-quality reads and primer/adapter sequences (QUADTrim; https://bitbucket.org/arobinson/ quadtrim) were removed to obtain single-end reads of length varying from 50-265 bp with an average of~80 bp. We sequenced thirty-three libraries in total, each genotype represented by two treatments, three time-points, and three biological replicates. The three samples of tolerant genotype (ICC8261) at the reproductive transition stage (RTS) were not sequenced due to RNA degradation.

Transcriptome Analysis
Sequenced reads were mapped to the Chickpea Reference Genome (Cicer arietinum v1.0) [28] by using Bowtie v2.2.3 [99] and Tophat v2.1.1 [100] with default settings. Then mapped reads were counted with HTSeq-count v0.61 [101], and differential expression analysis was performed by the EdgeR package [102]. This package runs in the R statistical environment, uses a count-based approach, and employs the over-dispersed Poisson model to account for both biological and technical variability. Finally, differentially expressed genes (DEGs) were selected by applying a cut-off of p-value ≤ 0.001, FDR ≤ 0.05, log2 fold change ≥ +1.0 and ≤ −1.0.