Initially, four cluster bean cultivars, namely, RGC-1025, RGC-1038, RGC-1055, and RGC-1066, were screened for their drought tolerance based on various parameters, such as germination, seedling growth, biomass, relative water content (Barrs and Weatherley, 1962), cell membrane injury (Leopold et al., 1981), malondialdehyde (MDA) (Hodges et al., 1999), total chlorophylls (Arnon, 1949), and total proline content (Bates et al., 1973).

Seeds of cluster bean (C. tetragonoloba L.) cultivar RGC-1025 were sterilized in 0.5% (W/V) sodium hypochlorite solution for 5 min, then rinsed thoroughly, and soaked in distilled water for 30 min. Seeds were sown in earthen pots containing soil and farmyard manure in a 3:1 proportion maintained in the departmental botanical garden. After 20 days post-sowing, drought stress was induced by withholding water to one set of pots, and respective fully watered controls were maintained in another set of pots. Ten days after stress imposition, fresh leaf samples from five plants were collected, pooled, flash-frozen in liquid nitrogen, and transported immediately to the sequencing facility.

For total RNA-seq, total RNA was extracted using the Qiagen RNeasy Plant Mini Kit with DNAse treatment (Thermo Fisher Scientific, United States) as per the manufacturer’s instructions. The quality and quantity of the RNA were estimated using a NanoDrop Spectrophotometer (Thermo Fisher Scientific, United States) and Qubit Fluorometer (Thermo Fisher Scientific, United States). The integrity of the RNA samples was analyzed using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). RNA-seq libraries were prepared with Illumina-compatible NEB Next^® Ultra™ II Directional RNA Library Prep Kit (New England BioLabs, MA, United States). In total, 500 ng of total RNA was taken for mRNA isolation, fragmentation, and priming. Fragmented and primed mRNA was subjected to first-strand synthesis followed by second-strand synthesis. The double-stranded cDNA was purified using JetSeq Clean Beads (Bioline Meridian Bioscience, Australia). Purified cDNA was end-repaired, adenylated, and Illumina multiplex barcode adapters were ligated as per NEBNext^® Ultra™ II Directional RNA Library Prep protocol, followed by second-strand excision using USER enzyme at 37°C for 15 min. Adapter-ligated cDNA was purified using JetSeq Beads and was subjected to 10 cycles for indexing (98°C for 30 s, cycling (98°C for 10 s, 65°C for 75 s) and 65°C for 5 min) to enrich the adapter-ligated fragments. The final PCR product (sequencing library) was purified with JetSeq Beads, followed by a library-quality control check using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). A total of 7,629,816 short reads were obtained, 150 bp length paired-end (P.E.) reads and average fragment size of ∼400 bp were used, and three biological replicates were processed for sequencing analysis with Illumina HiSeq™ 4000, which was out-sourced at Genotypic Technologies, Bengaluru, India.

The reads were processed for quality assessment and low-quality filtering before the FastQC tool assembly. The reads were then processed by removing the adapter sequences and low-quality bases (<q30) using the Cutadapt tool. Processed reads were assembled using a graph-based approach by the rnaSPAdes program. The characteristic properties, such as N50 length, average length, maximum length, and a minimum length of the assembled contigs, were calculated. De novo transcriptome assembly of the processed reads from all the libraries was done using Bowtie2 with end-to-end parameters. In the second step of the assembly procedure, clustering of the assembled transcripts based on sequence similarity is performed using the Cluster Database at High Identity with Tolerance (CD-HIT)-EST program¹ with 95% similarity between the sequences. This reduces the redundancy without excluding sequence diversity used for further transcript annotation and the DEG analysis.

All unigenes were annotated using the BLASTX search tool on Viridiplantae transcripts from the UniProt database containing 8,058,045 protein sequences and the NCBI non-redundant database (N.R.). The cutoff e-value was 10^––5, and the minimum similarity was more significant than 40%. Gene ontology annotation was carried out using the Blast2go program and visualized using Web Gene Ontology Annotation Plot (WEGO).²

DESeq, an R package, was used for differential expression analysis. Sequencing (variable library size/depth) bias among the samples was removed by library normalization using size factor calculation in DESeq. DESeq normalized expression values were used to calculate fold change for a given transcript. The regulation for each transcript was assigned based on log2-fold change. The transcripts that show a log2-fold change less than −1 are represented as downregulated. The values greater than one are upregulated and between −1 and 1 are termed neutrally regulated. Gene Ontology (GO) enrichment analysis and pathway analysis for DEG were done against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. KAAS server was used to analyze and characterize associated pathways. To obtain the highly significant differential expression genes, the criterion of the absolute value of reads per kilobase of transcript per million reads (RPKM) ratio > 1,000 was used.

Simple sequence repeats (SSRs) were identified using the MISA Perl script in each transcript (MIcroSAtellite identification tool).³ A simple repetition of motif length ranging from 1 to 6 bp was identified with recommended default parameters of MISA.

To evaluate the gene expression pattern from the DEG analysis, the total RNA extracted, as mentioned earlier, from control and drought-stressed cluster bean variety RGC-1025 was used for the SYBR Green qRT-PCR assay. The qRT-PCR assay was performed for 16 different stress-responsive genes selected from DEG analysis. These genes include aldo-keto reductase 1 (AKR1), late embryogenesis Abundant14 (LEA14), non-specific lipid transfer protein, TFs MYB30, NAC4, scare crow-like protein1 (GRAS TF’s), BHLH, GATA, malate dehydrogenase (MDH), aquaporin, DNA helicase, nitrate reductase, proline dehydrogenase (PRODH), serine hydroxy methyltransferase (SHMT), and thaumatin, trehalose 6-phosphate phosphatase (TRE6PH) with actin and tubulin genes as an internal control. For cDNA synthesis, 1 μg of total RNA from control and drought-stressed cluster bean RGC-1025 samples was treated with a Turbo DNASE treatment kit (Thermo Fisher Scientific, United States) as per the manufacturer’s protocol to remove any DNA traces. cDNA was synthesized using Revert Aid M-MuLV Reverse Transcriptase (Thermo Fisher Scientific) as per the manufacturer’s instructions. qRT-PCR mix was comprised of 1× using Power SYBR Green Master Mix (Ambion, United States), 20 ng of cDNA, and 0.2 μM of forward and reverse primers. Supplementary Table 1 shows the primers used for the investigated genes. The RT-PCR analysis was done on Applied Biosystems Step One Real-Time PCR machine with standard cycling comprising 95°C for 30 s, 40 cycles of 95°C for 1 s, 60°C for 20 s, and a melt curve analysis. Relative quantification was studied using 2^{–Δ Δ} CT method (Livak and Schmittgen, 2001). Each gene was analyzed in three biological samples, and three reaction replicates were performed for each biological sample.

Leaf surface waxes exteriorly deposited were extracted and quantified by a colorimetric assay reported by Mamrutha et al. (2010). Carnauba wax was used as the standard for the wax quantification assay. The wax content is represented as μg/gm fresh weight.

An scanning electron microscope (SEM) examined the epicuticular wax crystals. The third and fourth leaves of drought-stressed and control plants were cut to 0.5 cm, mounted onto standard stubs, and coated with gold particles using a fully automated vacuum spotter smart coater (DII-29030SCTR, JOEL, United States). The surfaces of the coated samples were observed through an SEM (JOEL-JSM-IT500, Japan).

Cluster bean cultivars, namely, RGC-1025, RGC-1038, RGC-1055, and RGC-1066, were screened for drought tolerance and we found significant differences at the cultivar level in morphological and biochemical traits under stress treatments. Results revealed that among four cultivars evaluated; cultivar RGC-1025 showed a lesser decrease in seedling growth, better biomass, and relative water content; lesser extent of cell membrane injury, lesser MDA, total chlorophylls content, and significantly higher levels of osmoprotectant, proline when compared to other cultivars (Supplementary Tables), suggesting the relative tolerance of cultivar RGC-1025 over other cultivars to drought stress. Therefore, we further extended transcriptome studies to understand the molecular mechanisms conferring the drought tolerance of cultivar RGC-1025.

To construct the transcriptome of cluster bean cultivar RGC-1025, high-quality RNAs from three replicates of drought-stressed and unstressed conditions (control) were sequenced. An average of 18 million paired-end reads were used for the downstream analysis after pre-processing. In total, 76,129,816 short reads were obtained using Illumina HiSeq 4000 Technology (Table 1). Most of the reads had > 99.9% score. Around 97.5% of the reads were retained in both the samples post-filtering. The cleaned-up reads were assembled using the Bowtie2 tool. CD-HIT was used to cluster redundant and similar isoforms. Finally, 66,838 transcripts were clustered with an average length of 955 bp. The non-redundant transcripts were considered as unigenes and were further analyzed.

All the obtained unigenes were annotated against the Uniprot Viridiplantae sequence database, NCBI non-redundant database, with a cutoff E-value of 10^–5. Around 55.98% of the unigenes were found to have hits in the public databases. The transcripts with more than 30% identity were considered during the analysis. These unigenes were classified as 31 functional categories. The significant fall under the category of DNA templated transcriptional regulation (2.34%) among the biological processes, the category of integral component of the membrane (25.4%) under cellular components, and the category of adenosine 5′-triphosphate (ATP) binding (13.14%) under molecular function. The highest transcript matches during the functional annotation with members of the family Fabaceae, such as Glycine (10,131), Mucuna (5,826), Cajanus (4,533), Cicer (2,632), Vigna (2,183), Phaseolus (2,111), Arachis (1,806), Trifolium (1,597), and Medicago (1,524) (Figure 1).

Differentially expressed genes were subjected to GO analysis to achieve functional classification. As a result, 37,418 DEGs fall into (i) molecular function, (ii) biological process, and iii) cellular components. In total, 30,900 (50.2%) DEGs were associated with molecular function terms, such as ATP binding encoding transcripts, followed by metal ion binding transcripts, DNA binding transcripts, zinc ion binding transcripts, nucleic acid binding transcripts, protein kinase activity transcripts, etc., and 16,410 (26.6%) DEGs were annotated with cellular component terms, represented by the integral component of the membrane encoding transcripts, followed by nucleus components transcripts, cytoplasm components transcripts, ribosome transcripts, plasma membrane transcripts, retrotransposon nucleocapsid transcripts, etc., and 14,191 (23%) DEGs were associated with biological process terms. This study also observed that identical DEG sequences could exist in more than one category. The most represented “biological process” subcategories identified were DNA-templated transcriptional regulation encoding transcripts followed by translation components, carbohydrate metabolic process transcripts, DNA integration transcripts, signal transduction transcripts, intracellular protein transcripts, etc. (Figure 2 and Supplementary Tables 2–4). In general, 12,866 genes were found to be upregulated, 16,177 genes were downregulated, and 27,782 genes had shown no change in their expression levels in cluster bean cultivar RGC-1025 due to drought stress. More interestingly, 3,745 transcripts were expressed only in the stressed sample.

Kyoto Encyclopedia of Genes and Genomes is an online database that deals with genomes and enzymatic pathways, and its identifiers were looked to predict biochemical pathways related to DEGs. Among 66,838 transcripts in unigene pathways, 17,211 transcripts against the KEGG pathways were identified. Of the 203 pathways identified, the top forty pathways are shown in Supplementary Table 5. DESeq analysis of the transcripts revealed that the enzymes with the most frequency of expression in the category of upregulated genes were alpha dioxygenase 2 (9.4-fold), low temperature-induced 65 kDa protein (LTI65; 9.2-fold), putative vacuolar amino acid transporter (9.05-fold), hexosyl transferase (EC 2.4.1.-; 7.95 fold), late embryogenesis abundant protein3 (LEA 3; 7.79-fold), Putative anthocyanidin 3-O-glucoside 2″-O-glucosyltransferase (EC 2.4.1.297; 7.44-fold), Glucosyltransferases, Rab-like GTPase Activators and Myotubularins (GRAM) domain protein/abscisic acid (ABA)-responsive-like protein (putative GRAM domain, P.H. domain-containing protein) (7.30-fold), and cytochrome P450 monooxygenase (EC:1.14.14.80; 7.14-fold). In the category of downregulated genes, the genes encoding the following proteins were found to be downregulated, such as putative CDP-alcohol phosphatidyl transferase class-I family protein 3 (EC 2.7.8.1; 0.5-fold), NEDD4-binding protein 2 (0.49-fold), dihydrolipoamide acetyltransferase component of pyruvate dehydrogenase complex (EC 2.3.1.-; 0.49-fold), putative cyclic nucleotide-gated ion channel 14 (0.498), ATP binding cassette (ABC) transporter G family member 28 (0.49-fold), N-(5-phosphoribosyl) anthranilate isomerase (putative phosphoribosyl anthranilate isomerase) (EC 5.3.1.24; 0.49-fold), F-box/FBD/LRR-repeat protein (0.49-fold), and adenylate kinase (EC:2.7.4.3; 0.49-fold). The heat map of DEGs showed the genes expressed (Figure 3).

Among the drought stress-responsive upregulated TF families, NAC family TF was the most abundant (26%), followed by MYB TFs (12%) and WRKY TFs (9%) (Figure 4A). Figure 4B depicts the top 40 upregulated TFs under drought treatments in cluster bean cultivar RGC-1025. Most renowned drought stress-responsive TFs upregulated include NAC4, NAC3, NAC29, NAC 104, and NAC18 from the NAC family, followed by WRKY12, WRKY50, WRKY6, WRKY33, WRKY30, WRKY24, WRKY42, WRKY53, WRKY70, and WRKY7 from WRKY family, further followed by other TFs, such as homeobox domain TFs, scarecrow/GRAS TFs, and ethylene-responsive TFs. Among the downregulated TF families, the hemophagocytic lymphohistiocytosis (HLH) TF family was the most abundant (20%), followed by the cellular TF family (10%) and homeobox domains (∼8%) (Figures 4C,D).

A total of 21,494 SSRs were identified in the cluster bean data, 14,434 (67.15%) are mono-nucleotide repeats, 3,072 (14.29%) are di-nucleotide repeats, 3,516 (16.35%) are tri-nucleotide repeats, 341 (1.58%) tetra-nucleotide repeats, 69 (0.32%) pentanucleotide repeats, and 63 (0.25%) hexanucleotide repeats (Figure 5). The 21,494 potential SSRs identified from de novo transcriptome sequencing data represent a significant addition to the limited set of genic-SSR markers available in cluster bean cultivar RGC-1025.

To verify the reliability of the expression profiles from RNAseq data, qRT-PCR analysis was performed for DEG (Figure 6). Results of the qRT-PCR assay for 16 stress-responsive genes revealed consistency in gene expression patterns as compared to that of the DEG analysis of the transcriptome of RGC-1025. The relative fold change of both the qRT-PCR and NGS-DEG is represented (Figure 6). Upregulation of stress-responsive TFs, such as NAC4, MYB30, scarecrow-like protein (SCL-1), primary helix-loop-helix TF (SlbHLH22), and TF bHLH 22, and the overexpression of candidate stress-responsive functional genes, such as DNA helicase, MDH, AKR1, LEA14, PDH, and SHMT, during drought stress by improving the ROS scavenging system, increasing osmotic potential, stomatal regulation, pH stability, respiration, and β-oxidation of fatty acids could support further the drought tolerance of cluster bean cultivar RGC-1025 (Ahmad et al., 2017; Waseem et al., 2019). The correlation among the expression patterns of the genes in NGS and qRT-PCR represents the consistency of the data in the current study.

The epicuticular wax content of the cluster bean cultivar RGC-1025 during the drought stress was 610.36 ± 0.53 μg/dm² as compared to the control 432.41 ± 0.4 μg/dm², which is 41.15% higher as compared to its control. These epicuticular wax data were supported by the SEM imaging of the leaf surfaces of control and drought-stressed plants. The SEM image of the leaf surface of RGC-1025 under drought stress showed enhanced wax crystal deposits (Figure 7).

Among the genes differentially expressed in wax biosynthesis pathway, there were upregulated genes that include KCS1 that encodes β-ketoacyl-CoA synthase 1 that is involved in elongation of 24C fatty acids, WSD1 that encodes wax ester synthase/diacylglycerol acyl transferase, which involves in wax ester biosynthesis, KCR1 that encodes β-Ketoacyl-CoA reductase, which is involved in very long-chain fatty acid elongation (VLFCA elongation), FATB that encodes acyl-acyl carrier protein thioesterase, which is engaged in supply of saturated fatty acids for wax biosynthesis, CER4/FAR3 that encodes alcohol forming fatty acyl CoA reductase, which is involved in formation of C24:0 and C26:0 primary alcohols, protein WAX2 encoding gene, CER17 also called Eceriferum1, which encodes for acyl-CoA desaturase-like 4 protein that is involved in n-6 desaturation of very long-chain acyl-CoAs, ABC transporter G family member 11, which encodes ABC transporter proteins that is involved in secretion of surface waxes in interaction with CER5, which is an another ABC transporter protein, and lipid transfer protein gene that encodes a lipid transport protein, which has role in cuticular wax export or accumulation. Finally, the upregulation of these wax genes in the present study through transcriptome DEG data reveals that these gene products are responsible for accumulating or producing epicuticular wax in cluster bean cultivar RGC-1025 (Figure 8). qRT-PCR analysis of selective wax genes showed significant changes in the expression patterns and an increase in the expression of the KCS1 gene was 1.85-fold, WSD1 gene was 1.81-fold, KCR1gene was 3.73-fold, FATB gene was 1.8-fold, CER4/FAR3 gene was 0.4-fold, protein WAX2 gene was 2.36-fold, CER17 gene was 2.89-fold, ABC transporter G family member11 gene was 1.95-fold, and lipid transfer protein gene was 1.49-fold.