Transcriptome analysis of lentil (Lens culinaris Medikus) in response to seedling drought stress

Background Drought stress is one of the most harmful abiotic stresses in crop plants. As a moderately drought tolerant crop, lentil is a major crop in rainfed areas and a suitable candidate for drought stress tolerance research work. Screening for drought tolerance stress under hydroponic conditions at seedling stage with air exposure is an efficient technique to select genotypes with contrasting traits. Transcriptome analysis provides valuable resources, especially for lentil, as here the information on complete genome sequence is not available. Hence, the present studies were carried out. Results This study was undertaken to understand the biochemical mechanisms and transcriptome changes involved in imparting adaptation to drought stress at seedling stage in drought-tolerant (PDL-2) and drought-sensitive (JL-3) cultivars. Among different physiological and biochemical parameters, a significant increase was recorded in proline, glycine betaine contents and activities of SOD, APX and GPX in PDL-2 compared to JL-3while chlorophyll, RWC and catalase activity decreased significantly in JL-3. Transcriptome changes between the PDL-2 and JL-3 under drought stress were evaluated using Illumina HiSeq 2500 platform. Total number of bases ranged from 5.1 to 6.7 Gb. Sequence analysis of control and drought treated cDNA libraries of PDL-2 and JL-3 produced 74032, 75500, 78328 and 81523 contigs, respectively with respective N50 value of 2011, 2008, 2000 and 1991. Differential gene expression of drought treated genotypes along with their controls revealed a total of 11,435 upregulated and 6,934 downregulated transcripts. For functional classification of DEGs, KEGG pathway annotation analysis extracted a total of 413 GO annotation terms where 176 were within molecular process, 128 in cellular and 109 in biological process groups. Conclusion The transcriptional profiles provide a foundation for deciphering the underlying mechanism for drought tolerance in lentil. Transcriptional regulation, signal transduction and secondary metabolism in two genotypes revealed significant differences at seedling stage under severe drought. Our finding suggests role of candidate genes for improving drought tolerance in lentil. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3596-7) contains supplementary material, which is available to authorized users.


Background
Lentil (Lens culinaris Medikus), a self pollinating crop with an approximate genome size of 4 Gbp [1] is an important legume which provides quality protein, carbohydrates, fibre and minerals for the humans and fodder for livestock. It is a moderately drought tolerant crop [2], but the yield is drastically reduced with increased drought stress. As water availability is important for crop growth and productivity, drought stress at critical stage with high severity can impose a threat to world food security. It ranks as the single most common cause of severe food shortages mostly in the developing parts of the world and represents far-reaching natural trigger of malnutrition and famine [3]. The mechanisms for drought tolerance in plants are very complex and highly variable [4]. Although, the tolerance to drought stress in lentil varies considerably among genotypes yet genomic information pertaining to drought stress is limited in this crop. For germplasm enhancement and to develop hardy lentil plants with tolerance to drought stress, understanding of molecular mechanisms governing response towards drought stress is necessary.
Drought stress induces a number of biochemical and physiological responses which are controlled by a number of genes at molecular, cellular and whole plant level that helps in maintaining water and ionichomeostasis and protect the plant from wilting and destined death. This can be achieved by maintaining osmotic compatibility within the cell, reconstruction of primary and secondary metabolism and by restoration of proteins in their native folded tertiary structure. Most well reported mechanisms of drought stress tolerance has been related to accumulation of n metabolites like proline, glycine betaine, soluble carbohydrates antioxidants etc. which help in maintaining the vital properties of the cell [5][6][7][8]. Several major classes of genes have been documented whose expressions are altered under drought stress. Most prominent among them are those involved in cellular metabolism, including cellular detoxification e.g., aldehyde dehydrogenase family genes; genes involved in cellular transport and signal transduction e.g. for ABA responses; genes encoding transcription factors which are involved in transcriptional regulation and genes for hydrophilic and heat-soluble proteins e.g., late embryogenesis abundant (lea) genes, etc. [9]. Functional genes comprised of heat shock proteins, facilitating protein refolding and stabilizing polypeptides and membranes under drought stress. Endogenous abscisic acid (ABA), content increases under drought stress which protects the plant from immediate desiccation by stomatal closure [10,11]. ABA has been shown to regulate expression of few genes under drought stress [12].
Recent developments of Next Generation Sequencing (NGS) technologies have enabled mass sequencing of genomes and transcriptomes, which produce a vast array of genomic information [13]. Genome wide expression studies provide to breeders a framework of dataset to understand the molecular basis of complex traits. Using NGS technology, Bett et al. have used lentil cultivar, CDC Redberry to develop an initial draft of 23x coverage which covered over half the lentil genome (2.7 Gb of the expected 4.3 Gb) [14]. Single nucleotide polymorphisms (SNPs) derived from this draft are available for use in molecular breeding of lentil. Kaur et al. have also performed sequence analysis in lentil using second generation sequencing technology and have developed a collection of expressed sequence tags (ESTs) [15]. In other legumes like chickpea, transcriptome analysis under drought stress has already been undertaken. Hiremath et al. used Roche/454 and Illumina Solexa to identify drought responsive genes and gene based molecular markers including simple sequence repeats (SSRs), SNPs and conserved ortholog set (COS) in chickpea [16]. A number of transcription factor families and defence related genes were identified in peanut under drought stress, using transcriptome analysis [17]. Wu et al. have found differentially expressed genes (DEGs) between terminal drought and optimal irrigation treatments in two different genotypes of common bean i.e. Long 22-0579 and Naihua which were functionally associated with drought stress [18]. However, information regarding transcriptomic changes under drought stress in lentil is very limited; therefore, to deduce pathways involved in drought stress response, expression study in contrasting genotypes of lentil is essential. Significant changes in gene expression are difficult to appraise without comparison. Therefore, PDL-2 and JL-3 which are two contrasting genotypes for drought tolerance, PDL-2 being drought tolerant and JL-3, drought sensitive were used in this study [19]. This study was undertaken to understand the biochemical mechanisms associated with adaptation to drought stress at seedling stage and to identify differentially expressed genes in contrasting lentil genotypes under drought stress using IlluminaHi-Seq2500 platform.

Plant material, cultivation and drought stress treatment
Two lentil genotypes: drought tolerant (PDL-2) and drought sensitive (JL-3) were included in this study. These two genotypes were selected on the basis of previous studies, reporting PDL-2 as a drought tolerant breeding line and JL-3 was drought sensitive on the basis of seedling survivability [19]. PDL-2, a breeding line derived from a cross between ILL-590 and ILL-7663, was obtained from International Center for Agricultural Research in the Dry Areas (ICARDA), Syria. JL-3 is a released variety of Central Zone of India and it is selected from landrace of Sagar district of Madhya Pradesh, India.

Hydroponic experiment
The hydroponic experiment was conducted at National Phytotron Facility, Indian Agricultural Research Institute (IARI), New Delhi, India in a completely randomized block design with three replications. Air temperature in the controlled environment was 22/18°C (2°C) day/ night; photoperiod was 10/14 h light/dark; and the relative humidity was approximately 45%.
Drought tolerance was evaluated by the protocol of Singh et al. in a nutrient solution culture [20]. Seeds were disinfected with 1% sodium hypochlorite for 2-3 min and rinsed thoroughly with distilled water and then germinated on filter paper. One week old seedlings were transferred to hydroponic medium (KNO 3 [21]. Drought stress was given one week after transplantation of seedlings. Two environments were created: Drought stress-the seedlings were exposed to air for 4 h for a period of 3 d. Control-the plants were kept in the nutrient solution for the entire period (3 d) of development without interruption. The pH of nutrient solution was adjusted at 6.5 with 1 M HCl or 1 M KOH. The solution was regularly aerated by an aquarium air pump and was replaced on alternate days. After 3 d of treatment, both the genotypes were rated for drought tolerance based on scale suggested by Singh et al. [20]. Data on relative water content (RWC), membrane stability index (MSI), photosynthetic pigments (total chlorophyll), proline, glycine betaine (GB), lipid peroxidation and antioxidant activities were taken after 3 d of drought stress treatment as described previously by Singh et al. [19,22,23]. The experiments were conducted in a completely randomized block design with three replications comprising twelve seedlings per replication.

Total ribonucleic acid (RNA) extraction and purification
Leaf samples were taken from the 12 seedlings for each genotype under control and drought stress treatments. Total RNA was extracted from leaves of drought tolerant line, PDL-2 and drought sensitive variety, JL-3 from both control and treated samples using QIAGEN RNeasy Plant Mini Kit. Ribosomal RNA was removed from total RNA by using Epicentre RiboZerorRNA removal Kit and Agencourt RNA clean XP Kit. Equal amounts of total RNA extracted from each seedling of each genotype were pooled together.
Complementary DNA (cDNA) library construction and sequencing Strand specific cDNA library was constructed using Tru-Seq RNA library preparation kit following the manufacturer's instructions. The quality of cDNA libraries was tested using Agilent 2100 bioanalyzer and quantified cDNA was subjected to purification using AMPure XP beads. Purified cDNA was then end repaired, using 3′ to 5′ exonuclease activity of End Repair Mix which removed 3′ overhangs and filled 5′ overhangs through its polymerase activity. cDNA was then polyadenylated and multiple indexing adapters were ligated to its end. For enrichment of cDNA in the library, polymerase chain reaction (PCR) was performed which selectively amplified those fragments that have adapter molecules on both the ends. The established cDNA libraries were sequenced using Illumina HiSeq2500 platform (Illumina Inc., San Diego. CA, USA) to generate 2 × 100 base paired end-reads.

Sequence analysis
Raw Fastq files obtained from the sequencer were checked for quality parameters of the sequences viz. base quality score distribution, average base content per read and GC distribution in the reads. The base quality score distribution and average base content per read were inferred using Phred quality score. The raw Fastq files were trimmed before performing de novo transcriptome assembly. First two bases and last ten bases were removed from all the reads. Fastq-mcf tool was used for removing adapter sequences. Reads of average quality score of less than 20 were also filtered out. Trinity with default options were used to assemble trimmed reads. Transcripts > = 200 bp were focussed for transcription expression estimation and downstream annotation. The trimmed reads were aligned to the assembled transcriptome (length > =200 bp) using Bowtie programme.

Screening and annotation of DEGs
Differential gene expression studies were performed using DESeq program. Transcripts having read counts > =1 and adjusted p value < = 0.05 were chosen for differential gene expression analysis. The read counts, transcript expression in Fragments Per Kilobase of transcript per Million mapped reads (FPKM) were evaluated in each stage. Edward plots were plotted to elucidate comparison of contigs in all the samples, both for upregulated and downregulated contigs. DEGs between treatments were also identified based on the FPKM value using with the edgeR program. We combined the statistical test with the multiple-hypothesis-testing correction method Benjamini and Hochberg [24], which calculates the False Discovery Rate (FDR), to qualify statistically significant, differentially expressed genes by avoiding inflation of type-1 errors.
The assembled transcript was annotated using CANoPI (Contig Annotator Pipeline). Assembled transcripts were compared with National Center for Biotechnology Information (NCBI) non-redundant protein database using BLASTx (Basic Local Alignment Search Toolx) programme. Matches with E-value < = 10 -5 and similarity score > = 40% were used for further annotation. For organism annotation, the top BLASTX hit of each transcript was studied and organism name was extracted. Gene and protein annotation was done as perthe matched transcripts. Among the total significant BLASTX hit transcripts 94,694 were annotated using UniProt (Universal Protein Resource) database and for the remaining ones, NCBI predicted protein annotation was done. Transcripts with proper gene name from UniProt and differentially expressed transcripts were shown with heat map. Gene Ontology (GO) terms for transcripts were extracted and were categorized into molecular function, biological process and cellular component categories. Enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways was performed using KOBAS (KEGG Orthology Based Annotation System). Reactome analysis was performed to predict the pathway information for DEGs.

Differential expression verification by Quantitative realtime PCR (qRT-PCR)
Accuracy of transcriptome sequencing data was validated by qRT-PCR. RNA was extracted from leaves of four biological replicates using QiagenRNeasy Plant Mini Kit. Quantification of RNA was done using Nano drop Spectrophotometer. Total RNA sample was reverse transcribed using Biorad cDNA synthesis Kit. Variable amount of total RNA for different samples was taken for final reaction volume of 20 μl. Normalization of all the cDNA samples was done so as to equalize concentration of 90 ng/μl. The primer sequences were designed using Primer3Plus and are listed in Table 1. The β-tubulin gene was used as a reference gene to normalize all the data. The 25 μl PCR reaction mixture was comprised of 4 μl diluted cDNA, 4 μl each of forward and reverse primers, 12 μl Evagreen dye. PCR amplification was carried out at 50 o C for 2 min, 95 o C for 10 min and 40 cycles consisting of 95 o C for 15 s and 60 o C for 1 min as run parameters. The relative quantification method 2 -(ΔΔCT) was used to calculate and calibrate the expression level of target genes in different treatments.

Filtering and alignment of SSRsand SNPs
SNPs were predicted using Samtools mpileup and custom scripts which call the variants based on read depth. A minimum read depth of 10 was used to filter heterozygous loci and false positive SNPs. In parallel to this, GATK toolkit was also considered to call variants using haplotype caller command version 3.6-0 for SNP calling keeping default parameters. Misa software was used for filtering SSRs from high quality filtered reads from both the genotypes that were aligned to the contigs and primers were designed using Primer3 software. Minimum primer size of the SSR primers designed using pri-mer3 software was 15 and maximum primer size was 21 with optimal primer size of 18. Estimated product sizes from these SSRs were in the range of 100 to 300 bp.

Variation in wilting and seedling survival
Differences between genotypes under control and drought stress were found significant for wilting and seedling survival ( Fig. 1a and c). The effects of stress were first observed on the leaves after 2 h of air exposure. The sensitive genotype (JL-3) showed wilting much earlier than the tolerant lines. Tolerant and sensitive genotypes severely wilted when exposed to air for 4 h, (Fig. 1b) and when returned into the nutrient solution and kept for 12 h, tolerant genotype (PDL-2) showed much faster recovery. On the other hand, JL-3 showed less recovery when returned to the nutrient solution (Fig. 1c).

Variation in physiological and biochemical traits
Relative water content RWC data showed significant reduction under drought stress in tolerant and sensitive lentil genotypes. However, PDL-2 maintained RWC under drought stress showing a significantly lower reduction (28.6%) in RWC compared to JL-3 (60.1%) under drought stress (Fig. 2a).

Membrane stability index
Membrane stability index was reduced under drought stress, although PDL-2 had lower reduction in MSI with 19.3% reduction over the control as compared to JL-3 which showed quite higher reduction of 57.7% (Fig. 2b).

Proline
Under drought stress, there was observed significant increase in proline concentration of PDL-2 (27.2% over control). On the other hand, sensitive lentil genotype JL-3 didn't show any significant change in proline accumulation (Fig. 3d).

Glycine betaine
In PDL-2, glycine betaine content increased significantly under drought with a 33.3% increase over the control, whereas in JL-3there was no significant change in its content (Fig. 2e).

Lipid peroxidation
Malondialdehyde (MDA) content was increased under drought stress in both the genotypes. However, the magnitude of its increment was prominent in JL-3 (36.5%) (Fig. 2f ).

Antioxidant activities Catalase (CAT)
No significant differences in the activity of catalase enzyme was observed ( Fig. 2j) under drought stress condition as compared to normal environment.

Superoxide dismutase (SOD)
There was higher SOD activity in the shoots of both tolerant and sensitive lentil genotypes under drought stress. SOD activity increased significantly in shoots of both the genotypes as compared to the control. The increase in SOD activity was higher in PDL-2 (32.5%) compared toJL-3 (6.1%) (Fig. 2g).

Ascorbate peroxidase (APX)
The activity of APX increased significantly under drought stress conditions in both tolerant and sensitive genotypes as compared to control. Seedlings grown under drought stress showed higher increase in APX activity in leaves of tolerant, PDL-2 genotype (44.8%) than sensitive genotype, JL-3 (15.9%) (Fig. 2i).

Glutathione peroxidase (GPX)
In drought conditions, GPX activity increased in the leaves of both the genotypes with a greater intensity in the leaves of PDL-2 (57.4%) than JL-3 (14.7%) (Fig. 2h).

Total RNA integrity and cDNA library preparation
Total RNA was extracted from the both control and drought stressed plants and the quality of RNA was tested using Nanodrop and Agilent Bioanalyzer 2100. All the samples were satisfactory for library construction and sequencing with RIN value in between 6.7 and 7.9. Concentration of RNA samples (ng/μl) were 1823, 1493, 3111 and 1708 for 1C, 1 T, 2C and 2 T respectively. The TruSeq RNA fragmentation protocol for transcriptome analysis was performed on RNA after mRNA purification using elevated temperature. The fragmentation resulted in libraries with inserts ranging from 120 to 200 bp with a median size of 150 bp. To keep the selection consistent in subsequent stages, the fragments were eluted according to bead volume and incubation time.

Analysis of differentially expressed genes
Differential gene expression of drought treated genotypes along with their controls revealed a total of 11,435 upregulated and 6,934 downregulated transcripts which were identified in the combinations of. 1C-1 T, 1C-2C, 1C-2 T, 1 T-2 T, 2C-1 T, 2C-2 T, where '1' and '2' represents drought tolerant (PDL-2) and drought sensitive (JL-3) genotypes respectively; 'C' represents control and 'T' as drought treated plants. Graphical representation of number of upregulated and downregulated transcripts in major comparative combinations is presented in Additional file 1: Figure S1. Data obtained from control and treated samples were subjected to combination wise comparison, from which DEGs were identified. A total of 6633 DEGs were identified (including 2919 solely expressed ones). Comparisons of total number of upregulated and downregulated contigs in all the genotypes are presented in Additional file 2: Figure S2 and Additional file 3: Figure S3), respectively in the form of Edward plots, which exhibited sharing of contigs between different combinations. When tolerant and sensitive genotypes were compared with their respective controls, 1514 and 1596 upregulated DEGs and 814 and 1012 downregulated DEGs were identified, respectively. When tolerant and sensitive genotypes were compared, 1417 upregulated and 1001 downregulated DEGs were identified. Upregulated DEGs constituted 65.03%, 61.20% and 58.60% in 1C-1 T, 2C-2 T and 1 T-2 T comparison groups, respectively. The results of significantly Differentially Expressed Genes obtained from DESeq & edgeR were compared and 45-55% of significantly expressing genes were found to be matching using both the software packages (edgeR & DESeq). When assembled transcripts were compared with NCBI nonredundant protein database using BLASTx, around 36% of transcripts were found to have confidence level of atleast 1e -5 , where E-value < = 10 -5 and similarity score > = 40% (Additional file 4: Figure S4a). Around 48% of assembled transcripts had similarity score of more than 60% at protein level with the existing proteins at NCBI database (Additional file 4: Figure S4b). To describe the results, two different analysis methods were applied. Firstly, major DEGs whose expression differed significantly in different combination groups were identified based on the criteria that there p-value < 0.005 and log2_foldchange > 1 (Tables 3, 4, 5, 6, 7 and 8). When compared between sensitive and tolerant genotypes under drought stress condition, some of the DEGs whose expression was significantly upregulated belonged to Cyclin family (c44790_g4_i5, c44790_g4_i4, c44790_g4_i3, c44790_g4_i2), Aldehyde dehydrogenase family (c40657_g1_i1), PTR2/POT transporter family (c40074_g1_i3), Fatty acyl-CoA reductase family (c43159_g1_i4, c43159_g1_i1, c43159_g1_i8), Phosphatase 2A regulatory subunit B56 family (c36460_g1_i2) and the ones with downregulated expression belonged to Adaptor complexes small subunit family (c46344_g1_i1), c46344_g1_i1 (c45522_g2_i4, c45522_g2_i2), SHMT family (c45167_g2_i11), Cystinosin family (c43466_g1_i10), Polyglycohydrolase family (c45115_g1_i3), Cation transport ATPase (P-type) family (c42686_g1_i6), MenG/UbiE family (c44280_g1_i2) ( Tables 5 and 6). In one of the comparison for 1 T-1C v/s 2 T-2C a total of 6720 genes were found significantly expressed with less than 0.05 FDR value (Additional file 5: Table S1).

Functional classification of DEGs
For functional classification of DEGs and to investigate the metabolic pathways in which they are involved for their fight against drought, the GO terms for transcripts were extracted and KEGG pathway annotation analysis  was performed on the identified DEGs. A total of 413 GO annotation terms were extracted where 176 were within molecular process, 128 in cellular and 109 in biological process groups (Additional file 11: Figure S10). When drought tolerant plants were compared to control, the top three significantly enriched GO annotation categories for downregulated DEGs were nucleus (GO:0005634), integral component of membrane (GO:0016021) and plasma membrane (GO:0005886). Apart from last two (GO:0016021, GO:0005886), third significantly enriched GO annotation for upregulated DEGs in tolerant plants was ATP binding (GO:0005524) (Additional file 12: Figure S11). GO annotation in sensitive genotypes revealed that categories similar to tolerant ones were found in sensitive genotypes also. When sensitive genotypes were compared with their controls, highly enriched GO annotation categories were integral component of membrane (GO:001602), plasma membrane (GO:0005886) and nucleus (GO:0005634) in downregulated transcripts along with GO terms GO:0016021, GO:0005886 and GO:0005634 in upregulated transcripts (Additional file 13: Figure S12). Also, when sensitive and tolerant genotypes were compared, 3 most significant GO categories were integral component of membrane (GO:0016021), ATP binding (GO:0005524) and plasma membrane (GO:0005886) in downregulated transcripts and nucleus (GO:0005634), integral component of membrane (GO:0016021) and ATP binding (GO:0005524) in upregulated transcripts (Fig. 6).
The pathway annotation analysis of a total of top 202 genes in the three comparison groups viz. 1C vs 1 T, 2C vs 2 T and 2 T vs 1 T revealed 28 unique drought related DEGs annotated to KEGG pathways including 42 metabolic pathways (Table 10) The most frequently associated pathways were metabolic pathways (15), followed by biosynthesis of secondary metabolites (13), carbon metabolism (6), citrate cycle and protein processing in endoplasmic reticulum (4 each). Some of the DEGs like AT3G55610 for arginine and proline metabolism; AT5G40390 for galactose metabolism were upregulated in both 1C vs 1 T and 2C vs 2 T comparison groups. Similarly, DEGs like AT3G27850 for propanoate metabolism; AT4G37370 for stilbenoid, diarylhetanoid, gingerol biosynthesis and limonene and pinene degradation; AT5G66760 for oxidative phosphorylation; AT3G48000 for fatty acid degradation and beta-Alanine metabolism etc. were upregulated in both 1C vs 1 T and 2 T vs 1 T comparison groups.
Some drought responsive transcription factors that bind to specific DNA sequences and control the rate of transcription were also identified. TFs are the key entities in transduction of stress related signals. When compared in 1C-1 T, 1 T-2 T and 2C-2 T comparison groups most of transcription factors belonged to TF families like AP2/ERF family (subfamily ERF or RAV), WRKY group II b and III families, HD-ZIP homeobox family or BZIP family of transcription factors. Also, several SSRs (9949), SNPs (8260) and INDELs (1248) were identified which could be further developed and used for drought related studies in lentil and other crops using Samtools mpileup toolkit. SNP calling analysis was performed by GATK toolkit, Haplotype caller tool (version 3.6-0) using default parameters (Additional file 14: Table S2).

Plant water stress
Hydroponic is most effective and practical approach for screening large number of genotypes in small area because it is easy in handling and has possibility of better controlled environment. Another advantage is that it is non destructive and plants can be screened at a early stage of growth, and tolerant plants can be selected and transferred to pots or field for further assessment of drought tolerance at subsequent stages of growth [19,20]. In the hydroponic technique, when plants were exposed to air for 4 h, all of them initially get wilted (Fig. 1). But the marked difference occurred when they were returned to the nutrient solution for 12 h, only the tolerant genotype (PDL-2) showed strong recovery, whereas the sensitive genotype (JL-3) showed less recovery. This suggested that even 4 h exposure to stress resulted in completely affected plant metabolism in the sensitive genotype, which was not observed in drought tolerant genotype. The exposure to air did not affect the basic metabolic activities of the plants and they retained the capacity to revive back to normal life when water became available [20]. The visual observations of wilting after relatively short durations of air exposure provided a suitable and reliable ranking of genotypes under longterm and/or more severe drought stress conditions. Therefore, it suggests that visual assessment of plant wilting and seedling survivability may offer suitable parameters for quick characterization of drought tolerance even at seedling stage.

Physiological and biochemical attributes
Plant RWC decreased significantly in both lentil genotypes under drought stress but PDL-2 maintained significantly higher RWC in both the control and water stressed conditions (Fig. 2a). It shows that the higher RWC enabled the tolerant lentil genotype to perform better in terms of physio-biochemical processes under water stress conditions. Tolerant genotype (PDL-2) produced higher stable yield, probably because water retention ability in plant is one of the components of tolerance mechanisms [19]. Higher RWC has been reported to be associated with higher level of photosynthetic pigments, membrane stability index, osmolytes and antioxidant activities in maize [25].
Opening of stomata results in more transpiration and subsequently reduction in RWC. Under these conditions, the genotype loses more water, particularly if drought is prolonged, plant recovery is impossible and ultimately it may die. Tolerant genotypes can maintain higher RWC in their leaves through stomatal closure and consequently reduction in leaf gas exchange [26]. The RWC in leaves of drought stressed cultivars decreased significantly. relative water content and water potential in the leaves under drought stress [27][28][29][30][31]. Drought stress caused a strong loss of photosynthetic pigments although; PDL-2 had higher chlorophyll contents than JL-3 under water stress. Contrasting effects of drought stress on the plant pigments have been reported previously [32,33]. Similarly drought-tolerant genotypes have been reported to maintain higher chlorophyll content than sensitive ones [34].

Many researchers have reported large reductions in
The synthesis of osmolytes including proline and glycine betaine is widely reported in plants to stabilize membranes and maintain the conformation of proteins at low leaf water potential. But variability in synthesis as well as accumulation of osmolytes occurs among intra and inter plant species. Our results showed that lentil genotypes which performed better under water deficit conditions had higher levels of proline and glycine betaine content as compared to those found sensitive to stress conditions ( Fig. 2d and e). Higher proline concentration has been reported to be involved in reducing the photo damage to chloroplast thylakoid membranes by scavenging and/or reducing the production of ROS [35][36][37]. Therefore, higher concentration of proline has been suggested as one of the parameters for selection for stress tolerant plant [38]. Similarly, proline has also been reported to protect and stabilize ROS scavenging enzymes and activate alternative detoxification pathways in plants subjected to various abiotic stresses [39]. Therefore, higher proline may act as a direct antioxidant as well as an activator of antioxidant mechanisms. Accumulation of proline has been associated with drought stress avoidance in maize, wheat and chickpea also [34,40,41]. Glycine betaine (GB) is another effective compatible solute which increases in the chloroplast of plants, when exposed to environmental stresses [42]. An increased accumulation of glycine betaine content was noticed in tolerant genotypes. Higher level of GB under drought stress has been reported in many plants like barley [43,44].
MDA has been suggested as a marker of oxidative stress-induced lipid injury and its concentration varies in response to abiotic stresses [45,46]. Lipid peroxidation, in turn, is an indicator of the prevalence of free radical reactions occurring in tissues and indicates a relationship between drought and oxidative stress [47,48]. The genotypes show better performance under water deficit conditions have been observed to have lower levels of MDA content in the roots (Fig. 2f) thus they protect themselves from lipid peroxidation of membrane systems as compared to the genotypes which had higher levels of MDA content. The similar results were observed in wheat and fababean under drought stress conditions [49,50].
Antioxidant enzymes like CAT and POX are key enzymes in scavenging and detoxification of hydrogen peroxide, a hazardous by product of photorespiration [51]. In this study the activity of these enzymes was induced due to drought stress. The increase in activity was higher in POX across genotypes responding to stress. Similar induced activity of above enzymes under stress conditions has been reported previously [51,52]. Highest activity of CAT in response to drought was observed in PDL-2 while lowest POX activity under stress condition was recorded in JL-3. This is clear that these genotypes experienced oxidative stress and their antioxidant enzymes triggered to detoxify cells. These results are in consistence earlier reports showing higher antioxidant enzyme activity in wheat cultivars [53,54].

Transcriptome analysis and stress responsive genes
Transcriptome sequences data can be a valuable resources especially for the species without a completely sequenced genome like lentil. In lentil drought stress mediated gene expression has been sparsely studied. In transcriptome analysis, N50 is an important criterion to decide the quality of assembly and with a N50 value of approximately 2000, quality of assembly was found to be very high [55][56][57]. A total of 11,435 upregulated and 6,934 downregulated transcripts were identified through differential gene expression of drought treated genotypes and their controls. Gene expression analysis revealed that genes involved in oxidation-reduction process, correct folding of protein, TCA cycle, electron transport chain, organ senescence and reduction of stomatal conductance are more severely upregulated in drought tolerant genotypes than the sensitive ones, whereas genes for transcription binding, GABA synthesis, synthesis of cell wall protein, those involved in negative regulation of absicisic acid etc. are downregulated in tolerant genotype as compared to sensitive ones. In tolerant genotype, the activity of unnamed protein product (NCBI accession no. 291047692, patented by Journal No. WO 2010020654-A2 25-FEB-2010) which belonged to aldehyde dehydrogenase family was most significantly upregulated with a log2 fold change of 7.9. Water stress often results in concentration of reactive toxic molecules like aldehydes, which can cause lipid peroxidation and alteration in proteins and nucleic acids. The aldehyde dehydrogenase family is a large family of enzymes which are regarded as "aldehyde scavengers" and irreversibly convert these aldehydes into acids which will result in less damage caused by several abiotic stresses including drought [58].
Other significantly upregulated DEGs within a log2 fold change of 3 to 6 belonged to 50S ribosomal protein L12-3, citrate synthase 3, succinate dehydrogenase [ubiquinone] flavoprotein subunit 1, Metallothionein-like protein and NADP-dependent malic enzyme 1. 50S ribosomal protein L12-3 is a chloroplast precursor and upregulation of its protein may represent the tolerant cultivar attempt to protect and sustain the correct folding of other protein in addition to accelerated degradation of unfolded/incorrectly folded or stress damaged protein [59]. Citrate synthase 3 is peroxisomal citrate synthase which is required for fatty acid respiration in seedlings where citrate is exported from peroxisomes into mitochondria during respiration of triacylglycerol [60]. Succinate dehydrogenase [ubiquinone] flavoprotein subunit 1 is located in mitochondrial respiratory chain complex II and is involved in mitochondrial electron transport. It was upregulated in Ilex paraguariensis leaves in response to water deficit and abscisic acid [61]. Metallothionein-like proteins are low molecular weight, cysteine-rich, soluble, and metal-binding proteins which are found in both plant and animal tissues. Involvement of these proteins in drought stress has also been confirmed in peanut [62,63]. Gorantla et al. found that metallothionein-like proteins represented the most abundant group of drought stressed transcripts in rice cultivar (Nagina 22) which helped in metal detoxification [64]. Similarly, two genes encoding metallothionein-like proteins were identified which were induced under drought stress in B. napus suggesting function of metallothionein-like proteins in drought stress [65]. NADP-dependent malic enzyme 1 plays an important role as anti-drought. The majority of water loss from plants occurs through stomata. When stomata are open, the concentration of potassium chloride and/or malate is high in guard cells, which enhance their turgor pressure and results in increased pore size. These are widely distributed in plant, which mainly appear in mitochondria, chloroplast as well as cytoplasm and catalyze the oxidative decarboxylation of malate to produce pyruvate, CO 2 and NADPH under metallic ions [66]. When leaf stoma are closed under drought stress, the malate concentration in cell decreases and NADP-ME activity increases [67].
In this study, it was found that under drought stressed condition for tolerant as well as sensitive genotypes there is higher percentage of upregulated DEGs. Further, some of the DEGs involved in TCA cycle, respiratory electron transport chain, ion channel transport, ABC family protein mediated transport, HSFs activation, metabolism of glucose are upregulated when tolerant and sensitive genotypes are compared with their controls whereas it has been reported that expression of genes involved in photosynthesis, photorespiration and carbohydrate metabolism were more drastically downregulated in drought tolerant genotypes than the sensitive ones [68].

Pathway analysis for drought tolerance
The detailed classification of KEGG pathways for DEGs with 3-7 log2fold change fall under 42 pathways, where a significant number of DEGs belonged to secondary metabolism category including stilbenoid, diaryl heptanoid and gingerol biosynthesis, phenylpropanoid biosynthesis, diterpenoid biosynthesis along with that of carbon metabolism and citrate cycle; protein processing in endoplasmic reticulum, biosynthesis of amino acids, galactose metabolism, plant hormone signal transduction, nitrogen metabolism including alanine, aspartate and glutamate metabolism, fatty acid degradation etc. Similar results have been reported in Ammopiptanthus mongolicus leaves under drought stress [69].
Antioxidant enzymes constitute "first line of defence" against ROS generated during environmental stresses like drought [70]. Superoxide dismutase is one of the major classes of antioxidants that catalyzes the first step in ROS scavenging system and removes superoxide radicals by converting them into H 2 O 2 and O 2 [71] . Higher SOD activity in drought tolerant lentil genotype, PDL-2 may be one of the mechanisms for drought tolerance as AT5G51100 gene was highly upregulated in tolerant genotype under drought stress, which is involved in hydrogen peroxide metabolism in peroxisomes. The production of H 2 O 2 can directly be countered by the activities of catalase and nonspecific peroxidise. Phenylpropanoids constitute a secondary antioxidant system, and activated upon due to depletion of primary antioxidant defences and control cellular H 2 O 2 within a sub-lethal concentration range [70]. Peroxidase 52 upregulated in this study are involved in biosynthesis of phenylpropanoids like p-Hydroxy phenyl lignin, Guaiacyl lignin and Syringyl lignin, which must have played a major role in ROS scavenging. Phenyl ammonia lyase activity which is a key enzyme involved in biosynthesis of isoprenoid antioxidant compounds was found to increase sharply in tolerant genotype under drought stress condition. Similar results have been reported in maize inbreds [72]. Biosynthesis of other secondary metabolities like stilbenoid, diaryl hepatanoid and gingerol were also found to increase due to upregulation of gene encoding cytochrome P450, polypeptide 8, which is involved in conversion of resveratrol to piceatanol, a crucial step involved in synthesis of these metabolites.
Stress-inducible galactinol-synthase, a member of glycosyl transferases family plays a key role in enhancing level of galactinol and raffinose under abiotic stress conditions, which is important osmoprotectants for drought stress tolerance [73]. Galactinol synthase is involved to transfer of UDP-D galactose to myo-inositol and is considered the main regulator of this biosynthetic pathway [74]. In this study, galactinol synthase 1 gene has been upregulated many folds in both 1C vs 1 T and 2C vs 2 T comparison groups suggesting its role in regulation of drought tolerance in lentil. Several studies have shown that the expression of galactose synthase gene is involved in response to several abiotic stress tolerance mechanisms [74][75][76][77].
Under drought stress conditions, the endoplasmic reticulum (ER) protein folding machinery reaches a limit and the demand for protein folding exceeds its capacity. Therefore unfolded or miss-folded proteins increase in the ER, and trigger an unfolded protein response. This results in up regulating the expression of genes encoding components of protein folding machinery or the ER-associated degradation system [78]. The same has been visualized in lentil as the expression of genes encoding heat shock protein 70, which is part of ubiquitin ligase complex and ER-associated degradation system and 70B, which is characterized as molecular chaperone and that of heat shock protein 22 (shsp) were significantly upregulated.
Gene encoding Glutathione S-transferase TAU 20 (AT1G78370), which is involved in glutathione metabolism, was upregulated in tolerant genotype under drought condition. Glutathione S-transferases are involved in protection under various stress conditions by detoxifying endogenous plant toxins which increases under oxidative stress [79]. GSTs play a major role under drought conditions by conjugation of glutathione with electrophilic substrates to increase their solubility and facilitating further metabolic processing [80]. There are reports, which suggest involvement of GSTs in drought stress [81,82]. AT1G78370, the gene upregulated in this study, has been shown to physically interact with Far-Red Insensitive 219 (FIN219) in response to light and play a crucial role in cell elongation and plant development [83].
Xyloglucanendo transglucosylase/hydrolase is important enzyme of cell wall and involved in modifying various physiological process for stress tolerance. Choa et al. have reported that constitutive expression of abiotic stress inducible hot pepper CaXTH3, which encodes Xyloglucanendo transglucosylase/hydrolase homolog and showed improved drought tolerance in transgenic Arabidopsis plants [84]. The same protein has been found to be increased in PDL-2 (tolerant) plant, through brassinosteroid induced plant hormone signal transduction pathway under drought stress in this study. Under drought conditions, plant needs to reduce shoot growth while maintaining root growth. This process requires differential cell wall synthesis and remodelling. Formation of reactive oxygen species and peroxidises are key players in this process, which initially cross link phenolic compounds and glycoproteins of the cell walls causing stiffening. Xyloglucan modifying enzymes results in cell wall loosening which allows further growth of stressed organs [85].
Absicisic acid dependent pathway can be considered important for drought tolerance in lentil as, gene encoding 9-cis-epoxycarotenoid dioxygenase which is a well known key enzyme for carotenoid biosynthesis, a precursor for ABA biosynthesis is highly upregulated in drought tolerant genotype under drought condition [86]. Increase in cellular ABA triggers the activation of several stress responsive genes and the closure of stomata to restrict transpiration [87,88]. The same has been confirmed in several plant species including Arabidopsis, cowpea, beans, etc. [88][89][90]. Genes involved in lysine degradation were also upregulated in drought tolerant genotype. In plants, lysine is catabolised to glutamic acid and acetyl CoA with the help of two enzymes viz. lysine ketoglutarate reductase and saccharopine dehydrogenase, whose gene expression was also involved in response to ABA and drought stress in soybean [91]. WRKY33 is involved in drought stress regulation; this gene was found to be upregulated in tolerant genotype under drought condition. WRKY33 was reported to be directly associated with drought tolerance through transcriptional regulation of Ces8A, an Arabidopsis drought tolerant gene [92]. Drought tolerant genes like MYB, ZFP, Metallothiorenin and others identified in our findings were compared to the list of genes identified and were found matching to the Rice database as well [93].