Transcriptomic signature of drought response in pearl millet (Pennisetum glaucum (L.) and development of web-genomic resources

Pearl millet, (Pennisetum glaucum L.), an efficient (C4) crop of arid/semi-arid regions is known for hardiness. Crop is valuable for bio-fortification combating malnutrition and diabetes, higher caloric value and wider climatic resilience. Limited studies are done in pot-based experiments for drought response at gene-expression level, but field-based experiment mimicking drought by withdrawal of irrigation is still warranted. We report de novo assembly-based transcriptomic signature of drought response induced by irrigation withdrawal in pearl millet. We found 19983 differentially expressed genes, 7595 transcription factors, gene regulatory network having 45 hub genes controlling drought response. We report 34652 putative markers (4192 simple sequence repeats, 12111 SNPs and 6249 InDels). Study reveals role of purine and tryptophan metabolism in ABA accumulation mediating abiotic response in which MAPK acts as major intracellular signal sensing drought. Results were validated by qPCR of 13 randomly selected genes. We report the first web-based genomic resource (http://webtom.cabgrid.res.in/pmdtdb/) which can be used for candidate genes-based SNP discovery programs and trait-based association studies. Looking at climatic change, nutritional and pharmaceutical importance of this crop, present investigation has immense value in understanding drought response in field condition. This is important in germplasm management and improvement in endeavour of pearl millet productivity.

SCIeNTIfIC RePoRTS | (2018) 8:3382 | DOI: 10.1038/s41598-018-21560-1 and cookies, for ethanol production in industries, used as poultry and animal feed and biofuel crop 6 . Since it is rich in antioxidants, it is used in nutraceutical industry 7 .
In a study projecting climate change till 2099, it is indicated that there will be extreme increase in drought selectively in Western Hemisphere (major part of Eurasia and Africa) but there will be contrarily high moist regions, from Alaska to Scandinavia 8 . Adversely, drought affected region also covers major pearl millet growing area decreasing the productivity at least by 25% 9 . In such geoclimatic situation, improved millet variety with further drought tolerance would be a preferred crop due to its high photosynthetic efficiency with low transpiration as C4 plant 10 .
Limited number of Quantitative Trait Loci (QTLs) 11 and genomic resources are reported in pearl millet. The reported genomic resources of pearl millet are in the form of simple sequence repeats (SSRs) 12 , expressed sequence tags (EST) 9,13 , EST-SSR [14][15][16][17] , ISSR 18 and SCAR 19 . Gene-based single nucleotide polymorphism (SNP) and Conserved-Intron Scanning Primer (CISP) marker 20 and DArT markers 21 . None of these are available in the form of web-based genomic resources. Looking at genome size, these markers are too less and moreover till now there is no trait specific genic region marker discovery using transcriptomic approach. Approaches like marker-trait associations, genomic selection, genome sequence and genotyping-by-sequencing are still warranted 22 . Genomic resources in pearl millet and management of abiotic and biotic stress are required to have insights of functional genes enabling crop improvement by molecular breeding approaches 23 .
Earlier work were based on heat treatment in phytotran 13 and chemically induced drought using hyperosmotic polyethylene glycol (PEG) 9 for induction of drought. All these were based on subtractive hybridization and cDNA (complementary Deoxyribonucleic Acid) library approach thus yielding very limited number of differentially expressed transcripts (<500) 13 . Moreover, tissues selected in previous experiments were confined to leaves only. For water and nutrient intake, root systems are most critical 24 thus they should also be included for better holistic view at physiological and molecular level. Moreover, root senses and responds to drought first 25 and mediates stress signals for root biomass adjustments, anatomical alterations, and physiological acclimations 26 thus transcriptomic investigation of root tissue is much more required 27 .
Since the effect of drought on root and leaf varies significantly in millets 28 , thus specific tissue approach is imperative for identification of pathway and marker discovery. None of these experiments were based on withdrawal of experiment which truly mimics the actual drought condition of field.
Transcriptomic approach of a single crop genotype with drought treatment has been reported successfully in tea plant deciphering major pathways regulating drought tolerance 29 . Such single genotype-based transcriptomic analysis has been successfully used to delineate major physiological response against drought in tomato 30 and cassava 31 also. In case of field crop like soybean, water deficit response has been found variety/ accession specific 32 . Since single genotypes-based experiment has been found more holistic in understanding the basic physiological mechanism operating at species level, thus such experiments are needed in investigation of drought response in pearl millet species also. Transcriptome database is available for large number of crops but no such database has been developed for orphan crop like pearl millet to be used as genomic resources in crop improvement research. Since whole genome of the pearl millet has recently been available 1 , thus development of such web resources can be done advantageously with SNP discovery by comparing with reference sequence along with genome annotation analysis. Such genic region putative markers (SSRs, SNPs and InDels) discovery can be rapid and cost effective which can be used in drought trait improvement program in future for better productivity of pearl millet.
The present study aims at identification of differentially expressed genes in leaf and root tissue of millet in response to drought induced by withdrawal of irrigation in the field. It also aims at identification of transcription factors (TFs), genic region putative markers viz., SSRs, SNPs and Indels (Insertions and Deletions), gene regulatory network (GRN) having hub genes and development of web-genomic resources.

Material and Methods
Seeds sowing in green house and drought treatment with recording of soil parameters.
Drought tolerant pearl millet J-2454 variety was sown in greenhouse during summer season of the year 2014-2015 at Agriculture Farm of Junagarh Agriculture University, Gujarat, India (21.5222° N, 70.4579° E, 107 meter above mean sea altitude). Seeds were sown in 2 kg polythene plastic bag under small greenhouse and polythene bag fill with equal weight soil mixture of sand, vermicompost and FYM in ratio of (40:40:20) and 25 to 30 seeds sown per polythene bag with three replication of one genotype to comparative study with control and water stress (or water withhold). Soil parameters, namely, average pH, electrical conductivity (EC), Maximum water holding capacity were recorded. Parameters were also recorded for water used in irrigation. During experiment in greenhouse condition day and night temperature and relative humidity (RH) were recorded. Seedlings were maintained by thinning post 10 days after sowing (DAS). Drought conditions were created by withdrawal of irrigation for 6 days after 23rd day of sowing of millet. Leaf and root tissues of pearl millet were collected after 6 days of withdrawal of irrigation (i.e., 29 th days after sowing) from control and drought treated plants. Withdrawal of irrigation at different time point can be an artificially created drought best suited by gene expression profiling in field crop 33 . In control plant, regular irrigation interval of alternate day was maintained. Major physiological information can be seen at transcriptomic level on 29th day of sowing as reported in other millet crop species 34 . Sampling at 29th day of sowing were as per the leaf and root tissues (two sets each) were collected from control and drought induced plants for transcriptomic studies. Tissues were frozen in liquid nitrogen and stored at −80 °C.
Plant tissue collection and RNA extraction. In order to minimize variability across samples, sample pooling approach was followed by taking ten biological replicates of both tissues, viz., root and leaf 35 . Total RNA was isolated from these tissues under control and drought condition using the standard protocol of TRIZOL RNA

Pre-processing, de novo assembly and identification of differentially expressed genes (DEGs).
The single-end Illumina reads were generated using root and leaf RNA of pearl millet genotype J-2454 having approximately 5678956 million and 6463411 million single-ends reads under normal and drought stress condition. The raw reads were pre-processed to remove any adaptor contamination using trimmomatic 37 software with parameters read length ≤36, poor quality ≤3 and HEADCROP:10 bases. These pre-processed data were used for transcriptome assembly along with the identification of differentially expressed genes and other analysis, transcription factor identification and putative genic markers' prediction. The different combinations used were root control vs. root under drought (RC, RT), leaf control vs. leaf under drought (LC, LT), root control vs. leaf control (RC, LC) and root under drought vs. leaf under drought (RT, LT).
The pre-processed high-quality reads were assembled using Trinity platform 38 . The abundance estimation for transcripts obtained was performed using 'RNA-Seq by Expectation-Maximization (RSEM) 39 . For the above mentioned transcriptome datasets, differentially expressed genes were identified using edgeR package 40 of Bioconductor. The significant DEGs were obtained with stringent parameters taking fold change value as two and FDR (False Discovery Rate) <0.05 41 .
Functional Annotation of transcriptome assembly. The sequence similarity search was conducted against the National Center for Biotechnology Information (NCBI) non redundant protein (Nr) database, and Swiss-Prot protein database using the BLASTx algorithm specifying E-values of less than 10 −3 . BLAST2GO 42  Gene Regulatory Network Analysis. The highly differentially expressed genes with fold value ≥8 from up and down regulated gene were considered for the gene interaction network. Network were visualized and carried out for further analysis using Cytoscape (version 3.2.1) 45 which is an open source platform for visualizing complex networks. ARACNE (Algorithm for the Reconstruction of Accurate Cellular Networks), a novel algorithm, specifically designed to scale up to the complexity of regulatory networks operating in the living cells was used. Network Analyzer plug-in was used to analyse the network centrality. The plug-in computes specific node centrality parameters and describing the network topology. Hub genes of complex networks were also obtained according to analysis of degree, betweenness and stress. The genes at the top of degree, betweenness and stress distribution were defined as hub genes.  Variant discovery (SNP and InDels) was done by both approaches namely, comparison of transcripts with de novo transcriptome assembly of genotypes J-2454 and also with available reference sequence of millet (genotype 23D2B1-P1-P5). For reference sequence based SNP discovery, genomic data of pearl millet Cenchrus americanus was downloaded from NCBI (https://www.ncbi.nlm.nih.gov/assembly/GCA_002174835.1/). Circos tool was used for generation of circular map of variants 48 .
Reads were mapped using BWA-0.7.5a (Burrows-Wheeler Aligner) and Samtools 49 . The SNPs were identified at read depth (d) ≥8 and quality depth (Q) ≥20 50 . Using the SAMtools program "vcfutils", SNP sites were further filtered, based on the criteria of 90 bp on both sides of the SNP in the alignment to ensure it in exon 51 . Validation and Expression Analysis by qRT-PCR. The first strand cDNA was synthesized from an aliquot of total RNA for each sample using RevertAid First Strand cDNA Synthesis Kit (Genetix, USA) and served as template for qRT-PCR (Quantitative Real Time Polymerase Chain Reaction). For quantitative PCR, 13 transcripts (6 for leaves and 7 for root) were randomly selected for primer designing using Primer 3 software 46 . The qRT-PCR was performed using QuantiFast SYBR Green PCR Master Mix (Genetix, USA) on ABI-7300 Real-Time PCR detection system, (Applied Biosystem) using standard 40 cycles along with melt curve step. Housekeeping gene, actin was used as endogenous reference for normalization. To obtain linear relationship, PCR conditions were optimized for each set of gene. Finally, differential gene expression were computed in terms of ΔΔCT fold change value 52 .
Web-genomic resource development. Pearl millet transcriptome database (PMDTDb) catalogues the information related to assembled contigs or transcripts, DEGs, the pathways in which these are involved, detailed SSR markers, and variants such SNPs and indels and miRNAs. It has three-tier architecture, i.e., client tier, middle tier and database tier. Web pages are developed for browsing the database along with the queries by user in client tier. All the information regarding contigs, markers, variants etc. are arranged in tables in various tables in MySQL in the database tier. For execution and fetching of user's query, scripting is done in PHP (Hypertext Preprocessor) in the middle tier. The pearl millet web-resources is available at http://webtom.cabgrid.res.in/pmdtdb/.

Standard of Reporting
Resource Identification. Germplasm resource used in the studies are completely disclosed by name of variety.

Results and Discussion
Greenhouse conditions and soil parameters. Soil and weather parameters were recorded. Irrigation withdrawal method was used successfully to induce drought. Average pH and EC of soil mixture were found to be-7.37 ± 0.01 and-1.18 ± 0.01 (ms), respectively. Maximum water holding capacity of soil mixture was found to be 30.15 ± 0.28%. Average pH and EC of water used in irrigation were found to be −7.07 ± 0.14 and −0.36 ± 0.01 (ms), respectively. Maximum day temperature (36-37 °C) and minimum night temperature (24-25 °C) were recorded along with day (84-86%) and night (49-51%) humidity. Other soil parameters and weather information are given in the Table 1.
Pre-processing and de novo assembly of transcriptomic data. Transcriptome data was generated successfully from both set of tissue samples, viz., leaf and root. A total of 12142367 SE reads were generated representing 2272632 (root, control), 3360164 (root, treated), 3406324 (leaf, control) and 3103247 (leaf, treated) reads. After data cleaning and quality assessment at Q ≥25 for control samples and Q ≥35 for treated samples, a final dataset was obtained. A total of 7563927 SE reads were represented by 1980048 (root, control), 2154770 (root, treated), 1216328 (leaf, control) and 2212781 (leaf, treated) reads. De novo transcriptome assembly using different assemblers viz., Trinity, MIRA, Cap3 and CLC revealed trinity to be the best one (N50 = 949, assembly size ~46.26 MB) and was considered for further analysis ( Table 2).  Table 3. Shared and unique DEGs are depicted in the Venn diagram ( Fig. 1) which shows that root is having more unique differential expressed genes (1444) 2). Out of these, 89 pathways were found common to all the four sets. The pathways most represented by contigs were purine and thiamine metabolism, biosynthesis of antibiotics, starch and sucrose metabolism, aminobenzoate degradation, glycine, serine and threonine metabolism and phenylpropanoid biosynthesis. Our analysis of top 30 common pathways reveals that leaf tissue are showing less metabolic activity with respect to root. This reflects that root has more DEGs than leaf for energy production and production of metabolites in response to drought for survival (Fig. 2).

Identification of Differentially Expressed Genes (DEGs
Higher purine metabolism observed in KEGG analysis is due to role of purine metabolism in Abscisic acid (ABA) accumulation mediating abiotic stress pathway of plant defence 53 . KEGG pathway analysis further reveals that, tryptophan metabolism was conspicuously high in both leaf and root tissues in response to drought. Abiotic stress in crop leads to accumulation of proline along with other amino acids. This is because of tryptophan's multifaceted role as osmolyte, ion transport regulation, stomatal control, detoxification reactions and redox-homeostasis 54 . This pathway analysis also reveals higher aminobenzoate degradation in both the tissues. Such higher degradation has been reported in abiotic stress of other crop like soybean 55 . This degradation increases the proline concentration in cell to protect cell from water deficit or abiotic stress. In case of drought     is higher for root. Again, this is due to higher cellular activity of root in response to drought. Such cellular activity involves cell organelles like mitochondria and glycosomes for catabolic metabolism. Higher membrane activity is inevitable for such cellular and catabolic activity 25 .
RAB gene reported to mediate polar root growth in response to drought stress in common bean 57 was found with highest magnitude among the upregulated genes in root. Similarly the role of other observed DEGs like serine acetyltransferase in sulfur assimilation pathway 63 , DHN9 (dehydrin 9) in membrane stabilizion, heat-shock gene (HSP17.8) and dehydrin 3 (DHN3) are involved in drought tolerance by controlling stomatal closure through controlling carbon metabolism 64 . We also observed differential expression of genes involved in energy balance and anti-oxidant activities which are already reported in drought response by various crops like SRP (stress responsive proteins) and DHN5 gene, controlling osmotic stress in wheat 65 , peroxidase controlling Reactive Oxygen Species (ROS) and antioxidant activities detoxifying. In fact SRP gene controls various other stress responsive genes 66 .
ROS concentration is increased in response to drought which may lead to cellular damage but gets detoxified due to antioxidant activities of peroxidase gene 58 .
In drought response, there is change in energy balance and metabolic pattern. Genes controlling these activities were found differentially expressed in our dataset, for example, LEA (Late Embryogenesis Abundant) gene   reported to be associated with salt stress 67 , ATP (Adenosine Triphosphate) citrate synthase reported to control energy balance in drought 68 , aspartate kinase-homoserine dehydrogenase reported to control tricarboxylic acid (TCA) cycle 69 , glycolysis and Krebs cycle controlled by NADP dependent malic enzyme gene 70 . High affinity of nitrate transport uptake has been reported in drought treated roots and the same was observed in differential expression of nitrate transporter gene 59 . Unplaced accessions 1088 Table 7. Chromosome wise distributions of Variants in Pearl millet reference genome. , respectively. We found 7 miRNAs viz., ssp-miR444a; osa-miR169d; ssp-miR444b.2; bdi-miR169l; osa-miR444f; osa-miR414; hvu-miR169 common in all the four sets (Fig. 3). This is due to wide conservation of miRNA having diverse functions in seedling growth and also in response to abiotic stress 77 . Few miRNAs may have multiple targets of different genes in the regulatory networks 78 . Thus, this is obviously expected with at least in some extensively conserved miRNA we may get but still they may have diverse functions. miR444a is reported to mediate drought stress for its adaptation in wheat. It is also known to interact with MADS-box transcription factors, which is associated with stress response 79 . Besides, miR444 responded to drought stress in Dongxiang wild rice 80 . miR169d plays role in fighting drought stress in cotton 81 as well as sorghum 82 . Also, miR169 family responds differently to drought in plants 83 . Gene Regulatory Network Analysis. In root GRN analysis, our dataset reveals response of drought by major intracellular signal transduction mediated by Serine Threonine Protein (STP) kinase which are actually Mitogen Activated Protein Kinases (MAPK) 84 . Drought triggers low expression of this gene to lower MAP Kinase activity. High affinity potential transporter gene was found with lower expression in drought. This genes acts as a sensor for drought in root and also controls uptake of potassium. This gene has been reported to be downregulated in root to balance potassium translocation 60 . Drought leads to high meristematic activity in root mediated by lower expression of O methyltransferase zrp4 85 . Similar low expression of this gene was found in root in response to drought. In drought, ROS gets accumulated in root tip and NADH controls its meristematic activity through ABA pathway by energy balance 86 . This gene was also found differentially regulated in root in response to drought. The details of root GRN hub genes is given in Table 4. Also the leaf subnetwork of important hub DEGs are shown in Fig. 4.

Transcription Factors identification.
In leaf GRN analysis, our dataset reveals Dehydrin gene family is well known as candidate genes associated with drought tolerance in crop. Among these gene family, in our dataset, we found differential expression of genes viz., dehydrin DHN1 and dehydrin COR410-like along with isoform of dehydrin COR410. Dehydrin is known to play role in drought response for crucial protective functions of root tissue. They are produced in response to ABA pathway mediated abiotic stress 87 . This gene family has consensus amino acid seq KIKEKLPG which is reported to be associated with differential expression along with its various isoform in response to drought 88 .
The gene, late embryogenesis abundant protein (LEA) is also reported to be a candidate gene of drought having larger gene family. In our dataset, we found up-regulation of late embryogenesis abundant protein, group 3-like isoform X2, late embryogenesis abundant protein D-34, late embryogenesis abundant protein 3, late embryogenesis abundant protein Lea5. This gene family is reported to play protective role in cells under drought 89 . Two larger gene families, viz., dehydrin and LEA were found differentially expressed in our dataset with its isoform. Similar observation has been reported in durum wheat in desiccation stress 90 . The details of leaf GRN hub genes is given in Table 5. Also the leaf subnetwork of important hub DEGs are shown in Fig. 5.  Table 6. Maximum abundance (%) was of trinucleotides repeats (51.62) followed by mononucleotides (24.3), dinucleotides (19.82), tetranucleotides (3.48), pentanucleotides (0.62) and hexanucleotides (0.14). SSRs were also mined from 4 sets of root and leaf tissues and found 385, 244, 618 and 508 markers in (RC, RT), (LC, LT), (RC, LC) and (RT, LT) respectively (Table 6). These putative SSRs markers can be validated by genotyping them in various cultivars. These putative SSRs markers can be validated by genotyping them in various cultivars. Out of 4192 genic region SSRs, we could successfully design primers for 2828 for de novo full assembly (Additional file 5). Our reported SSR can be further used as functional domain marker in millet variety improvement program especially with respect to drought tolerance. Such genic region based SSR markers from leaf associated with drought has been reported in mulberry plant 91 . Such genic region SSR markers have an advantage over genomic region, viz., transferability, a priory information about gene itself with known functionality. All these are desirable in crop improvement program 92 . SNP and InDel identification. SNPs and INDels were discovered successfully by both approaches.
Transcriptome based approach using a stringent pipeline, we identified 9318 total variants, having 5587 single nucleotide polymorphisms (SNPs) and 3736 InDels using the SAMtools software (Additional Files 6). In reference based SNP discovery, a total of 18360 SNPs and InDels were obtained. Chromosome wise SNP distribution shows highest SNP over chromosome 2 and lowest number of SNP over chromosome 4 (Table 7). Circular map was also generated to depict chromosome-wise SNP distribution among two genotypes of millet (Fig. 6). More SNPs were found by reference based method than transcriptome based. This is obviously expected as transcriptome based SNP represents intra-genotype variation due to presence of heterozygotes only whereas inter-genotype variation based SNP discovery includes SNPs obtained by alignment of sequences between two genotypes from throughout the genome. This has led to discovery of higher number of SNPs by reference based method. Such alternative alignment has been reported as an efficient method of SNP discovery using two genotypes 93 .
Putative SNPs in pearl millet transcriptome could be due to millet being predominantly protogynous crop, having high cross pollination resulting into very high heterozygosity also 94 . Being field crop, pearl millet shows high genetic variability within a single open pollinated variety contributing SNP discovery in transcriptome based approach. If the crop has large effective population size and its wind-pollinated reproductive trategy, then there is further increase in heterozygosity yielding these putative SNPs 95 . Besides this, the potential source of SNP could be variation between each of the two sets of 10 plantlets pooled for RNA extraction for transcriptomic data generation. SNPs detected in silico may be either true SNP showing allelic polymorphisms or it may be gene duplication or paralogous or homologous sequence variations 96 . Further observed variants (SNP and InDels) might be due to various other factors such as alignment ambiguity and undetected paralogs 97 . Transcriptome based putative SNPs can be further validated by Sanger sequencing followed by detection of double peak in chromatogram using SNP discovery by heterozygote approach 98 .
These discovered informative SNPs will be of immense use in designing of various throughput genotyping assay using their respective chromosomal location and genomic co-ordinates over reference sequence. Genic region SNP discovery has been used in various crop trait improvement programs like cold tolerance 99 and dormancy 100 in wheat, RIL for mapping in Brassica 101 , rust resistance in switchgrass 102 , cold tolerance in sorghum 103 and blister rust in white pine 104 . Such genic region SSRs and SNPs from candidate genes and hub genes of GRN are valuable genomic resource for eQTL discovery, high-density mapping and trait improvement in various crops like common bean 105 and chickpea 106 . Validation and Expression Analysis by qRT-PCR. Relative gene expression value obtained by qRT-PCR analysis of all the 13 genes (6 leaf, 7 root) were in correspondence with computed log fold change value of DEGs (Additional File 7, Fig. 7).
Web-genomic resource development. Drought associated web genomic resource of pearl millet, PMDTDb (http://webtom.cabgrid.res.in/pmdtdb/) has been developed using transcriptome data. It catalogues assembled contigs or transcripts, 19983 DEGs, 7596 transcription factors and a total of 34652 genic region putative markers (SSR markers, SNPs and InDels). The flowchart for its usage is shown in Fig. 8. For other species like foxtail millet, similar genomic resources with molecular markers 107 and transcription factors database 108 has been developed. For another species finger millet (Eleusine coracana (L.) Gaertn.), such transcriptome sequence has been reported to provides insights into drought tolerance and nutraceutical properties 109 . Based on the present study, we report the first web-based genomic resources of pearl millet. These resources can be used for further in candidate genes-based SNP discovery programs and trait-based association studies in drought improvement.
Potential application of millet genomic resources. For varietal improvement, traditional breeding program can be supplemented with molecular breeding approach by using genomic resources 110 . Our web-genomic resource is having 17856 SNPs (mined from two millet genotypes) which can be used in variety improvement. Such SNP discovery using even single genotype has been reported in forage crop like Artemisia 95 . Since there is no report of SNP discovery in candidate genes associated with drought in millet thus further research can be done by targeted sequencing of the reported genomic resources in large number of varieties and populations. Such approach can reduce the time and cost required for resequencing of large genotypes without prioritizing the candidate genes involved. Similar use of DEG-based SNP discovery has been reported in other crops like wheat for abiotic stress tolerance 94,95 , Switchgrass for biotic stress tolerance 102 . Candidate genes depicted in GRN can be a preferred source of SNP discovery 111 required for future association studies. Further, such transcriptome-based approach has been reported for eQTL discovery in development of high-density map like Brassica rapa 101 .
Our web-based genomic resources is also having 4192 SSRs along with ready to use primers for genotyping can be used in variety identification and improvement program. Similar transcriptome-based SSR has been used in various crops like barley 112 , S. tuberosum 113 , sugarcane 114 , capsicum 115 , eggplant 116 and differentiation of Basmati rice from non-Basmati 117 and also for DUS (Distinctness, Uniformity and Stability) testing for varietal identification 118,119 . Our enlisted SSR can be further used in MAS programme of millet improvement as reported in sorghum 120 , tagging stem rust resistance gene Sr35 in wheat 121 , Saltol QTLs in rice 122 .
Genomic resources of miRNA and its target in PMDTDb can be further used in research. It can be useful for both knowledge discovery (mechanism/regulation of drought responsive genes) and application oriented research especially for variety improvement. Since in finger millet, it has been reported that drought tolerance can be increased by use of gene silencing of drought associated miRNA thus our enlisted predicted miRNA can be used for similar work 123 . It has been widely reported that miRNA can be used to enhance crop yield along with increased tolerance to biotic and abiotic stress. Further, genome editing tool like CRISPR-cas9 can be used to control expression of miRNA 124 , thus there is greater scope for further research on enlisted miRNAs/genes by genetic modification for crop improvement. Such genome editing approach has been reported to be very promising in reducing anti-nutrients in millet, thus making it further enriched cereal 125 .
Our enlisted TF genes can be used for future SNP discovery for traits improvement of crop. Similar approach has been reported in oil plant for regulation of dwarfism by negative regulation of DELLA protein 126 . In case of Prunus, traits like the flowering, fruit quality, and biotic and abiotic stress resistance have been found regulated with TF 127 .

Conclusion
Present work reports root and leaf transcriptomic signature of drought induced by irrigation withdrawal in pearl millet. We found 19983 DEGs, 7596 TFs and GRN having 45 genes controlling drought response. A total of 34652 genic region putative markers viz., 4192 SSRs, 12111 SNPs and 6249 InDels are reported. Validation of gene expression by 13 randomly selected genes was in correspondence with computed FPKM values. We report major candidate genes such as LEA, Dhn, ATP-citrate synthase family, peroxidase, stress responsive protein and Aspartate Kinase-Homoserine Dehydrogenase genes. Enlisted candidate genes can be used for further SNP discovery programs and association studies. Looking at climatic change and nutritional and pharmaceutical value of this crop as well as its genetic potential of resilience, present investigation is of immense value in understanding drought response of millet in field condition. Such information is important in germplasm management and improvement in endeavour of pearl millet productivity.