Transcriptome landscape of perennial wild Cicer microphyllum uncovers functionally relevant molecular tags regulating agronomic traits in chickpea

The RNA-sequencing followed by de-novo transcriptome assembly identified 11621 genes differentially xpressed in roots vs. shoots of a wild perennial Cicer microphyllum. Comparative analysis of transcriptomes between microphyllum and cultivated desi cv. ICC4958 detected 12772 including 3242 root- and 1639 shoot-specific microphyllum genes with 85% expression validation success rate. Transcriptional reprogramming of microphyllum root-specific genes implicates their possible role in regulating differential natural adaptive characteristics between wild and cultivated chickpea. The transcript-derived 5698 including 282 in-silico polymorphic SSR and 127038 SNP markers annotated at a genome-wide scale exhibited high amplification and polymorphic potential among cultivated (desi and kabuli) and wild accessions suggesting their utility in chickpea genomics-assisted breeding applications. The functional significance of markers was assessed based on their localization in non-synonymous coding and regulatory regions of microphyllum root-specific genes differentially expressed predominantly in ICC 4958 roots under drought stress. A high-density 490 genic SSR- and SNP markers-anchored genetic linkage map identified six major QTLs regulating drought tolerance-related traits, yield per plant and harvest-index in chickpea. The integration of high-resolution QTL mapping with comparative transcriptome profiling delineated five microphyllum root-specific genes with non-synonymous and regulatory SNPs governing drought-responsive yield traits. Multiple potential key regulators and functionally relevant molecular tags delineated can drive translational research and drought tolerance-mediated chickpea genetic enhancement.

Scientific RepoRts | 6:33616 | DOI: 10.1038/srep33616 C. microphyllum, a member of tertiary gene pool is presumed to have originated from India and Turkey. The recent marker-based genetic diversity and population genetic structure studies documented a higher admixture and close phylogenetic relationship of one contrasting C. microphyllum accession domesticated primarily in the natural habitat alongside glaciers of the Hindu Kush Himalayas with cultivated desi chickpea of Indian origin 2,11 . The agro-ecological distribution of this C. microphyllum accession characterized by a peculiar feature of deep tap root system in the diverse rocky soil topography of Himalayan region, infers its strong adaptation towards drought and cold climatic conditions. Furthermore, this perennial wild accession exhibiting a vital characteristic of dense pod dehiscence is believed to serve as a potential source of pod borer resistance in chickpea 18 . The existing diversity among the wild C. microphyllum accessions for pod borer resistance can be thus exploited to broaden the genetic basis of resistance with an objective of driving marker-assisted disease resistance breeding and crop improvement in cultivated chickpea 18 . Collectively, a wild relative C. microphyllum accession revealing close evolutionary relationship with cultivated desi accessions is a potential natural source of novel genes and alleles especially related to drought, salinity and cold tolerance as well as pod borer resistance in chickpea.
In view of aforesaid possibilities, the present study employed a RNA-seq strategy to optimize de novo high-quality transcriptome assembly and generate comprehensive transcript sequences of an Indian origin perennial wild C. microphyllum accession for the first time at a high-resolution scale. A comparative analysis of root and shoot transcriptomes between C. microphyllum and a widely studied drought tolerant cultivated desi chickpea accession ICC 4958 39,40 was performed to identify molecular tags like candidate genes, transcription factors (TFs) and key regulators/metabolic pathways especially governing inherent root-mediated natural adaptation characteristics of C. microphyllum in adverse agro-climatic condition. The informative transcript-derived genic SSR and SNP markers developed at a genome-wide scale were further utilized for molecular mapping of major QTLs governing two drought-responsive important yield traits (yield per plant and harvest-index) on a constructed high-density inter-specific genetic linkage map (transcript map) of chickpea. The high-resolution QTL mapping was further integrated with comparative transcriptome profiling to delineate candidate genes underlying major drought yield QTLs in chickpea.

Results and Discussion
Transcriptome sequencing of C. microphyllum. To generate the global transcriptome of a wild perennial C. microphyllum accession, the transcriptome sequencing of root and shoot tissue samples-derived RNAseq libraries was performed using an Illumina NGS (next-generation sequencing) platform. These libraries were constructed by using the high-quality RNA isolated from three independent biological replicates of root and shoot tissue samples, each comprising pooled tissues from three randomly selected seedlings of a C. microphyllum accession. The natural adaptation of C. microphyllum accession towards abiotic stress (drought, salinity and cold) tolerance in adverse agro-climatic conditions is primarily because of its inherent characteristic feature of deep tap root system required for survival of this perennial wild chickpea in the diverse rocky soil topography of Himalayan region 18 . Considering the above facts, the root and shoot tissues of C. microphyllum were utilized specifically for transcriptome profiling to obtain the maximum representation of genes/transcripts differentially expressed/regulated in this hardy wild perennial accession as compared to a cultivated desi (ICC 4958) chickpea during natural adaptive climatic conditions. The RNA sequencing of both root (mean 54.7 million reads) and shoot (64.1) tissues of C. microphyllum accession generated on an average of ~118.7 million paired-end sequence reads across all three replicates of studied samples (Table 1). Of these, on an average ~105.9 million high-quality filtered sequence reads (root: 48.6 and shoot: 57.3 mean million reads) with a mean Phred quality score ≥ 30 at an individual base position were obtained ( Table 1). The sequencing data generated from three independent biological replicates of each root and shoot tissue sample was highly correlated (91%) based on average Pearson's correlation coefficient estimated among samples. These observations assured high reproducibility among biological replicates despite generating variable number of high-quality transcript sequence reads by transcriptome sequencing of C. microphyllum roots and shoots. All the C. microphyllum transcriptome sequence resource Scientific RepoRts | 6:33616 | DOI: 10.1038/srep33616 generated in the present study was submitted to NCBI-SRA (sequence read archive) database under BioProject ID PRJNA31486 with accession IDs SRR3340324-3340327 for unrestricted use.
Assembly and mapping of C. microphyllum transcriptome. With an effort to constitute a reference-guided assembly of the C. microphyllum transcriptome, the filtered high-quality transcript sequence reads were mapped primarily on the reference draft genome sequences of desi (ICC 4958) chickpea. On an average 55% (ranging 49-61% in shoots and roots) of transcript sequence reads derived from root and shoot samples of C. microphyllum were mapped uniquely across the chickpea genome. Subsequently, to increase the mapping efficiency of uniquely mapped sequence-reads and generate the optimal assembly of C. microphyllum transcriptome, the de novo assembly using all high-quality non-redundant sequence reads (by removing the duplicate reads) obtained from both root and shoot samples was performed. The summary statistics defining the best de novo transcriptome assembly outputs of C. microphyllum optimized in our study, are mentioned in the Table  S1. The mapping of on an average about 94 million de novo assembled high-quality transcript sequences uniquely suggests the use of at least 89% of transcript sequence reads in de novo assembly of C. microphyllum. The optimized assembly resulted in generation of 154352 unique transcript sequences with the average transcript and N50 length of ~896 and ~1641 bp, respectively (Table S1). The mapping efficiency of de novo assembled high-quality (longer average transcript and N50 length) transcript sequence reads constituted by us on the C. microphyllum transcriptome is higher/comparable to that documented earlier in reference-, de novo-and hybrid (reference and de novo) assembly-based transcriptome studies in crop plants including chickpea [19][20][21][22][23][24][25][26][27][28][32][33][34][35]38,41 . Interestingly, about 85% of these improved de novo assembled transcripts of C. microphyllum were mapped to reference desi chickpea genome. The use of 154352 unique high-quality transcript sequences in Cufflink followed by Cuffmerge assembly generated 38149 putative gene loci. The comparison of these candidate gene loci with the genes annotated from desi chickpea genome based on BLAST homology searches (with 90% sequence identity and 70% query coverage) exhibited orthology of 11219 (29.4%) unique C. microphyllum genes with 8465 (27.8%) candidate protein-coding genes (total 30257) annotated from desi genome.

Differential transcriptome responses of C. microphyllum.
To determine the transcriptome responses based on differential gene expression profiling, the high-quality transcript sequence reads generated from each root and shoot tissue sample were mapped on our optimized de novo assembled C. microphyllum transcriptome. This led to the mapping of ~89% transcript sequence reads on the C. microphyllum transcriptome. The abundance of transcripts (referred as genes) from each sample was estimated in accordance with FPKM (fragments per kilobase of exon per million fragments mapped) followed by differential expression profiling through Cuffdiff. This analysis identified 11621 genes exhibiting significant (≥ two-fold at P ≤ 0.05) differential up (7592, 65.3%)-and down (4029, 37%)-regulated expression specifically in roots as compared to shoots of C. microphyllum. Remarkably, 2514 (21.6%) including 1597 (63.5%) and 917 (36.5%) genes showed a pronounced higher (≥ 10-fold at P ≤ 0.05) differential up-and down-regulation, respectively in roots than that of C. microphyllum shoots. Higher accumulation of differentially regulated transcripts especially in the roots of C. microphyllum may underlie its deep tap root system-mediated inherent characteristics of natural adaptation towards drought, salinity and cold tolerance in adverse agro-climatic condition. Therefore, these differentially regulated molecular tags especially expressed in the roots of C. microphyllum have functional relevance to decipher the regulatory gene (TFs) networks and/or metabolic pathways controlling abiotic stress tolerance traits in wild as well as cultivated chickpea.
The comparative analysis of transcriptomes between C. microphyllum and ICC 4958 by correlating the differential expression profiles of genes in the roots and shoots of these accessions was performed. The detail statistics regarding mapping of transcript sequence reads obtained from control unstressed root and shoot tissues of ICC 4958 on the C. microphyllum transcriptome, are mentioned in the Table 1. This identified a total of 12772 genes  38 ). Average sequence reads were obtained by estimating the mean of reads generated from three independent biological replicates of each sample.
differentially up-and down-regulated (≥ two-fold at P ≤ 0.05) under at least one root and shoot tissue sample of C. microphyllum and ICC 4958 (Fig. 1A, Table S2). This included 463 genes differentially expressed commonly between C. microphyllum and ICC 4958. Higher number of root-specific genes (3685) as compared to shoot-specific genes (2100) were differentially expressed in C. microphyllum and ICC 4958. A total of 3242 rootand 1639 shoot-specific genes exhibiting pronounced higher differential expression specifically in C. microphyllum were identified (Fig. 1A, Table S2). Notably, 6987 genes were differentially expressed commonly between roots and shoots of both C. microphyllum and ICC 4958. This included 16 genes commonly expressed between C. microphyllum and ICC 4958, whereas 6740 genes specific to C. microphyllum (Fig. 1A, Table S2). Interestingly, we observed about 7.3 times higher differential up-and down-regulated expression of genes in the roots of C. microphyllum (3242 genes, 88%) as compared to that of ICC 4958 (443 genes, 12%). Therefore, these abundant root-specific differentially expressed genes especially scanned from the C. microphyllum could be vital to decipher the possible key regulatory pathway/mechanism governing superior root-responsive natural adaptive characteristic of this hardy wild chickpea accession in adverse climatic condition (diverse rocky soil topography of Himalayan region) as compared to a known cultivated drought tolerant desi accession ICC 4958.

Expression dynamics of genes and transcription factors.
The KOG-based determination of putative function for 12772 differentially regulated genes identified from comparative transcriptome profiles of roots and shoots between C. microphyllum and ICC 4958 indicated primary roles of 3229 (25.3%) genes in multiple cellular, biological and molecular processes in crop plants. This revealed enrichment of differentially expressed genes (especially root-and shoot-specific genes) involved in translation, ribosomal structure and biogenesis (J, 709 genes, 22%), post-translational modification, protein turnover and chaperones (O, 411, 12.7%), and energy production and conversion (C, 275, 8.5%) ( Figure S1). Among 441 (  Despite higher root-specific gene expression, more than three times differential expression of TF-encoding genes specifically belonging to abundant class of TF gene families like MYB, bHLH, NAC and WRKY in C. microphyllum (136, 76.8%) than that of ICC 4958 (41, 23.2%) was observed. In contrast, we observed more than 2.5 times higher expression of root-specific genes especially belonging to genetic information processing (transcription, translation and post-translational modifications) in C. microphyllum than that of ICC 4958. In this context, remarkable differential regulation/alteration of diverse class of gene transcripts especially in C. microphyllum roots might impart enhanced adaptation to this hardy wild chickpea in adverse climatic condition (Himalayan region) than that of a known cultivated drought tolerant desi accession ICC 4958. The differential transcript accumulation and regulation of TF-encoding genes (especially MYB, NAC and WRKY TFs in this study) representing diverse TF families have implicated their common and/or specific role in multiple adaptive trait mechanism functioning in many crop plants (including chickpea) during acute climatic adversities 38,[42][43][44] . Molecular characterization and functional validation of these TFs and dissection of their regulatory networks would further assist us to understand the global transcriptional reprogramming and thereby useful in defining the precise role of TF genes underlying adaptive trait mechanism in chickpea.

Regulation of molecular pathways.
To determine the possible gene regulatory pathways/interactions involved in natural adaptation responses, the pathway analysis of 12772 genes/TFs differentially expressed in roots and shoots of C. microphyllum and ICC 4958 was performed using KEGG database. Of these, 2212 (17.3%) genes were found to be significantly enriched with diverse pathways related to metabolism [primary (carbohydrate, lipid, nucleotide, amino acid, fatty acid, photosynthesis and abscisic acid) and secondary (flavonoids and terpenoids) metabolism 656, 29.7%], genetic information processing (transcription, translation and post-translational modifications, 595, 26.9%) and environmental information processing (signal transduction and membrane transport, 183, 8.3%) ( Figure S4). Remarkably, more than two times higher altered expression of above-mentioned multiple pathway-related genes in the roots of C. microphyllum than that of ICC 4958 was observed suggesting potential involvement of these genes in imparting root-mediated natural adaptation to this wild chickpea. The pronounced differential regulation of genes involved in various primary and secondary metabolic pathways and post-translational modifications (phosphorylation and ubiquitination) are majorly known to elicit diverse cell signalling and transcriptional events, and thereby play crucial role in various plant adaptive abiotic stress responses 38,42,[45][46][47][48][49] . Collectively, the differentially expressed root-and shoot-specific genes and TFs scanned from the comparative transcriptome profiles between known stress tolerant wild C. microphyllum and cultivated desi ICC 4958 accessions could serve as an immense transcriptomic resource for elucidating the key gene regulatory mechanism involved in differential abiotic stress tolerance adaptation traits especially between cultivated and wild chickpea. This will assist us to select best suitable differentially induced/regulated potential genes (especially root-specific genes) from C. microphyllum for developing climate resilient cultivars and thereby accelerate translational genomics in chickpea.

Experimental validation of gene expression. For experimental validation of differentially up-and
down-regulated genes identified though transcriptome profiling, the genes exhibiting high root-specific expression in C. microphyllum and ICC 4958 were used for differential expression profiling. For this, the RNA isolated from the root and shoot tissues of six drought tolerant (C. microphyllum, ICC 4958 and ICC 296131) and sensitive (ICC 4951, ICCX-810800 and ICCV 93954) desi and kabuli chickpea accessions including a C. microphyllum accession (used for transcriptome sequencing) was amplified with the gene-specific and endogenous internal control (EF1α) primers using the semi-quantitative and quantitative RT-PCR assays. The up-and down-regulated differential expression profiles of genes observed in roots and shoots of both C. microphyllum and ICC 4958 by RT-PCR assay was strongly correlated (85%) with that of our RNA seq-led transcriptome profiling study (Fig. 3, Table S2).
To enrich the informativeness of designed genic SSR markers in C. microphyllum, these markers were compared with that derived from a desi chickpea accession ICC 4958 based on variation of the repeat-units, which led to develop 282 in silico polymorphic genic SSR markers. The fragment length polymorphism detected by in silico polymorphic genic SSR markers between C. microphyllum and ICC 4958 varied from 2 to 106 bp with an average of 6.6 bp. The di-nucleotide class I (132 markers, 87.4%) in silico polymorphic genic SSR markers were found most abundant than that of tri-nucleotide class I (102, 78.5%) and class II (28, 21.5%) SSR markers (Fig. 4F). The in silico polymorphic 257 (91.1%) and 25 (8.9%) SSR markers were physically mapped on the eight chromosomes and unanchored scaffolds of the desi chickpea genome, respectively (Fig. 4G). Highest and lowest number of markers were mapped on the chromosomes 6 (50 markers) and 8 (16 markers), respectively. The physically mapped 282 polymorphic SSR markers were structurally and functionally annotated in different coding and non-coding regulatory sequence components of 282 desi genes. Maximum proportion of polymorphic SSR markers were derived  Table S2).
Validation and polymorphic potential of genic SSR and SNP markers. To experimentally validate the genic SSR markers, 192 C. microphyllum transcript-derived SSR markers and 96 in silico polymorphic (between C. microphyllum and ICC 4958) SSR markers physically mapped on chromosomes were genotyped in 96 wild and cultivated chickpea accessions using agarose gel-based assay (Table S3). This included six known drought tolerant and sensitive accessions selected in our study for differential expression profiling. Two hundred sixty-three of 288 markers (used for experimental validation) produced single reproducible PCR amplicons with an average amplification success rate of 91.3%. One hundred eighty seven of 263 amplified SSR markers exhibited showing in silico fragment length polymorphism between C. microphyllum and ICC 4958 based on repeat-unit variation were experimentally validated successfully using gel-based assay. One hundred eighty-seven polymorphic SSR markers validated in our study were localized in diverse coding and regulatory sequence components of differentially expressed root-specific genes/TFs of C. microphyllum and ICC 4958 (Tables S2 and S3). This implicates the functional relevance of genic SSR markers (including in silico polymorphic markers) developed from the C. microphyllum transcriptome at a genome-wide scale for large-scale genotyping applications including QTL/eQTL (expression QTL) mapping and association analysis to identify potential genes/QTLs controlling important agronomic traits in chickpea.
The validation of 127038 SNPs mined in silico between C. microphyllum and ICC 4958 exhibited the presence of 1124 (0.9%) SNPs common between our present and past studies 11,25,26,50,51 based on their nature/types of SNP alleles and physical positions (bp) on the desi (ICC 4958) genome (Table S4). Interestingly, 30 SNPs in the 26 genes of these were differentially expressed in roots of C. microphyllum and ICC 4958. The SNPs discovered from the genes as well as differentially expressed genes exhibiting differentiation between C. microphyllum and ICC 4958 could serve as an informative marker resource to be deployed in genomics-assisted breeding applications of chickpea. To experimentally validate the C. microphyllum transcriptome-derived genic SNPs, 192 genome-wide well-distributed (physically mapped on chromosomes) SNPs were genotyped in 96 wild and cultivated chickpea accessions using the MALDI-TOF MassARRAY genotyping assay. This included six known drought tolerant and sensitive accessions selected in our study for differential expression profiling. One hundred eighty-three (95.3%, mean PIC: 0.41) of 192 non-synonymous and regulatory SNP loci were successfully validated (polymorphic between C. microphyllum and ICC 4958) and exhibited polymorphism among 96 wild and cultivated (desi and kabuli) chickpea accessions with 95.3% validation success rate. Notably, these experimentally validated informative SNPs were derived from diverse coding and regulatory sequence components of differentially expressed root-specific genes/TFs of C. microphyllum and ICC 4958 (Tables S2 and S4). This suggests the functional significance of genic SNPs discovered from C. microphyllum transcriptome at a genome-wide scale for high-throughput genetic analysis in chickpea.

Functional significance of genic SNPs.
To determine the functional significance of genic SNPs more elaborately, comparative analysis of transcriptomes by correlating the genes differentially expressed in the roots and shoots of C. microphyllum with the control and stress (drought, salinity and cold)-imposed root and shoot tissue samples of ICC 4958 was performed. The detail statistics regarding mapping of transcript sequence reads obtained from ICC 4958 on the C. microphyllum transcriptome are mentioned in the Table 1. This identified a total of 13077 including 3580, 3532 and 3524 C. microphyllum genes that were differentially up-and down-regulated (≥ two-fold at P ≤ 0.05) under at least one root and shoot tissue sample/stress condition of ICC 4958 during drought, salinity and cold stress, respectively (Table S5). The C. microphyllum genes with coding and regulatory SNPs differentially regulated in control and stress-imposed root and shoot tissues of ICC 4958 were correlated with the TFs, KOG and KEGG-based functional annotation information of these genes ( Table 2). The TFs-encoding C. microphyllum genes belonging to different abundant class of TF families such as MYB, NAC and bHLH that were differentially expressed in roots and shoots of ICC 4958 during drought, cold and salinity stress responses contained more than 1.5 times higher frequency of regulatory SNPs than non-synonymous SNPs. Likewise, higher frequency (2.8 to 4-times) of regulatory SNPs as compared to non-synonymous SNPs in the abiotic stress-responsive differentially expressed genes representing abundant class of KOG and KEGG-based functional categories was observed (Table 3). Henceforth, it will be interesting to correlate the differential expression profiles of genes with the SNPs mined from their gene regulatory regions for understanding the molecular basis of altered gene transcriptional regulation underlying abiotic stress responses in wild and cultivated chickpea.
Our global transcriptome sequencing-based comparative differential gene expression profiles in shoots and roots of C. microphyllum and/or ICC 4958 accessions under drought, salinity and cold stress responses was correlated with the occurrence of diverse non-synonymous/large-effect and regulatory SSR and SNP markers in the coding and non-coding sequence components of these genes (Table S6). This analysis identified the presence of 457 SSR and 5687 SNP (including 706 non-synonymous, 9 large-effect and 1948 regulatory SNPs) markers in the 381 and 579 C. microphyllum genes, respectively that were differentially expressed in shoots and roots of ICC 4958 during drought, salinity and cold stress responses. Especially, 67 and 141 C. microphyllum genes with 86 and 1388 SSR and SNP markers, respectively were differentially expressed in roots of ICC 4958 during drought, salinity and cold stress responses. Of these, 21 and 22 genes with 27 SSR and 295 SNP markers, respectively were differentially regulated especially in roots of C. microphyllum (Table S5). Notably, 86 and 217 C. microphyllum genes with 125 SSR and 1930 SNP markers, respectively were differentially expressed in roots and shoots of ICC 4958 commonly across drought, salinity and cold stress responses (Table S6). Of these, 25 and 53 C. microphyllum genes with 31 SSR and 479 SNP markers, respectively were differentially expressed in roots of ICC 4958 commonly across drought, salinity and cold stress responses whereas 10 genes with markers differentially regulated in roots of C. microphyllum (Table S6). The details regarding SSR and SNP markers-containing C. microphyllum genes responsive to at least one specific stress condition (drought, cold and salinity) in both the roots and shoots of ICC 4958 and especially in the roots of ICC 4958, are provided in the Table S6.
Collectively, these analyses delineated a diverse set of root-and shoot-specific differentially up-and down-regulated abiotic stress-responsive genes and TFs with regulatory and non-synonymous SNPs exhibiting differentiation between C. microphyllum and ICC 4958 (Table S6). Henceforth, these gene-based informative SNP markers can be utilized in rapid targeted mapping of drought, cold and salinity stress-responsive differentially expressed root-and shoot-specific genes on the chromosomes by high-throughput marker genotyping and high-resolution QTL mapping in the advanced generation mapping populations of chickpea. Consequently, this will accelerate the delineation of functionally relevant molecular tags like QTLs, eQTLs (expressed QTLs) and potential genes especially regulating abiotic stress tolerance traits in chickpea. The SNPs being derived from the   Table 3.

Characteristics of an inter-specific chickpea genetic linkage (transcript) map comprising eight chromosomes constructed using a RIL mapping population (ICC 4958 × ICC 17163).
Scientific RepoRts | 6:33616 | DOI: 10.1038/srep33616 diverse coding (non-synonymous and large-effect) and regulatory sequence components of drought, salinity and cold stress-responsive differentially expressed genes can also have potential in rapidly establishing marker-trait linkages for identification of potential genes/QTLs and natural allelic variants governing abiotic stress tolerance by high-throughput marker genotyping and high-resolution genetic association mapping in natural germplasm lines of chickpea.

Construction of a high-resolution inter-specific chickpea genetic linkage (transcript) map.
To construct a saturated inter-specific genetic linkage (transcript) map, 500 SSR and SNP markers derived from the drought-responsive differentially expressed root-specific genes of C. microphyllum and ICC 4958 exhibiting polymorphism between parental accessions (ICC 4958 and ICC 17163) were genotyped among 190 individuals of a RIL mapping population (ICC 4958 × ICC 17163). The linkage analysis utilizing 500 marker genotyping information mapped 490 genic marker loci onto eight LGs (linkage groups) of an inter-specific transcript map of chickpea (Table 3, Fig. 6). This eight LGs-based transcript map covered a total map length of 825.4 cM with an average inter-marker distance of 1.68 cM (Table 3, Fig. 6). Longest and shortest transcript map length spanning 112.5 and 72.8 cM were observed in LG5 and LG8, respectively. Maximum (75 markers) and minimum (41) numbers of genic markers were mapped on LG6 and LG8, respectively (Table 3).
LG6 contained a most saturated transcript map (an average inter-marker distance of 1.48 cM) whereas LG2 had the least saturated map (1.99 cM) ( Table 3).
Collectively, an inter-specific transcript map constructed in the present study revealed a comparable/higher mean map-density (an average inter-marker distance of 1.68 cM) than that previously reported in multiple SSR and SNP markers-anchored intra-and inter-specific integrated genetic linkage maps of chickpea 11,17,25,26,32,40,[51][52][53][54] . Therefore, the integrated SSR and SNP markers-based transcript map constructed in our study has adequate map-density to be utilized as a reference genetic linkage map for high-resolution molecular mapping of QTLs governing diverse agronomic traits including abiotic stress tolerance traits in chickpea.  (Table S7). A highly significant positive correlation between YP and HI traits (r: 98% at a P < 10 −4 ) based on estimated Pearson's correlation coefficient (r) was observed. The QTL mapping using the genotyping information of 490 gene-derived SSR and SNP markers mapped on an inter-specific transcript map and field phenotyping data of 190 RIL mapping population identified six major genomic regions harbouring six significant (6.5-11.5 LOD) QTLs associated with YP and HI traits which were mapped on six chromosomes (except chromosomes 1 and 8) of chickpea (Table 4, Fig. 6). These identified QTLs exhibited consistent and stable phenotypic expression with more than 10% PVE (phenotypic variation explained) across rainfed and irrigated environments of two geographical locations and thus were considered as robust QTLs governing YP and HI traits in chickpea. The six major genomic regions underlying robust QTLs covered (1.5 cM on chromosome 5 to 5.0 cM on chromosome 4) with 20 SNP and SSR markers were genetically mapped on LGs/chromosomes and showed positive additive gene effect for increasing yield and harvest-index with large effective allelic contribution from ICC 4958. The proportion of PVE by individual robust QTL varied from 11.2-23.7%. The combined PVE estimated for all six robust QTLs was 39.8% (Table 4, Fig. 6). Six QTLs associated with multiple traits (YP and HI) were mapped on six different genomic regions at the similar marker intervals of chromosomes. The high-resolution QTL mapping identified SSR and SNP markers in the diverse drought-responsive differentially expressed root-specific genes tightly linked to the major drought yield QTLs (YP and HI) which are mentioned in the Table 4 with corresponding marker genetic positions (cM) on the chromosomes. To ascertain the validity of identified robust QTLs, the major genomic regions harbouring six YP and HI QTLs were compared with that documented by previous inter/intra-specific mapping population-led QTL mapping studies 39,40 . None of the YP and HI QTLs identified by us exhibited correspondence with previously reported known drought-responsive yield QTLs based on their congruent physical positions on chickpea chromosomes. This implicates the novelty and population-specific characteristics of six major QTLs identified in our study for dissecting the drought tolerance-related yield trait regulation in chickpea. Therefore, six markers-containing drought-responsive root-specific genes tightly linked with the major YP and HI QTLs mapped on chromosomes could have functional relevance to be deployed in marker-assisted genetic enhancement for developing high-yielding drought tolerant chickpea cultivars.

Molecular mapping of drought yield QTLs in chickpea.
An integrated genomic approach to delineate candidate genes underlying drought yield QTLs in chickpea. Six robust YP (CaqYP2.1, CaqYP3.1, CaqYP4.1, CaqYP5.1, CaqYP6.1, CaqYP7.1 and CaqHI7.1) and HI (CaqHI2.1, CaqHI3.1, CaqHI4.1, CaqHI5.1, CaqHI6.1 and CaqHI7.1)-associated major QTLs genetically mapped on six chromosomes identified by high-resolution QTL mapping were selected to delineate candidate gene(s) regulating drought tolerance-related yield traits in chickpea (Table 4, Fig. 6). The integration of genetic linkage map information of markers flanking the six major QTLs with that of physical maps of desi chickpea genome defined genomic regions (spanning 518931-4888587 bp) harbouring such robust QTLs on six chromosomes ( Table 4). The structural and functional annotation of these physical QTL intervals on six chromosomes identified multiple (19 to 230 genes) candidate protein-coding desi chickpea genes ( Table 4). The differential expression profiles of genes annotated in the major genomic regions harbouring six robust QTLs were determined using the aforementioned comparative transcriptome profiling data generated by us in the roots and shoots of C. microphyllum as compared to that of control and drought imposed root and shoot tissues of ICC 4958. Interestingly, non-synonymous and URR SNPs-containing six C. microphyllum genes/TFs revealing tight linkage with six major YP and HI QTLs (based on single marker analysis-led high-resolution QTL mapping) Scientific RepoRts | 6:33616 | DOI: 10.1038/srep33616 exhibited drought-responsive root-specific expression as well as pronounced up-regulated expression (6.4 to 46.8-fold) especially in ICC 4958 based on our comparative transcriptome profiling study. Our quantitative RT PCR-based gene expression profile validation data also revealed up-regulated (> 10-fold) differential expression of these six SNPs-carrying genes especially in the roots of known drought tolerant (C. microphyllum, ICC 4958 and ICC 296131) chickpea accessions as compared to that of sensitive (ICC 4951, ICCX-810800 and ICCV 93954)  Table 4. Orange and blue boxes indicate the QTLs governing yield per plant and harvest-index, respectively mapped on the chromosomes of a high-density transcript map. The SNPs used as anchor markers to construct a transcript map are illustrated with black colour lines. The non-synonymous/regulatory and synonymous SNPs derived from the differentially expressed genes are indicated with red and blue colour lines, respectively.
Scientific RepoRts | 6:33616 | DOI: 10.1038/srep33616 accessions (Fig. 3). Interestingly, more than three-fold up-regulation of six genes particularly in the roots of a known drought tolerant hardy wild C. microphyllum accession in contrast to other two known drought tolerant cultivated accessions (ICC 4958 and ICC 296131) was observed implicating the efficacy of delineated genes in regulating drought tolerance of chickpea. The five desi genes (Serine threonine protein kinase, ABC transporter G-family, Cytochrome P450, Chalcone synthase and Glutathione-S-transferase) and one TF (WRKY72) with non-synonymous and URR SNPs regulating drought tolerance-related yield traits scaled-down in this study, by integrating high-resolution transcript-derived SSR and SNP markers-based QTL mapping with comparative transcriptome profiling are functionally relevant. Based on previous documentation, these six genes/TFs are known to play crucial role in various abiotic (drought) stress responses by regulating diverse DNA binding and transcriptional activity, hormonal cell signalling, transport, metabolism and accumulation of transcripts/proteins in crop plants [55][56][57][58][59][60][61][62][63][64][65][66][67][68][69] . One of these gene coding for serine threonine protein kinase governing root-mediated drought tolerance-related yield trait has been validated by combining high-density QTL mapping with marker-trait association and gene enrichment analysis in chickpea 70 . The informative non-synonymous/large-effect and regulatory SSR and SNP markers as well as potential differentially expressed root-specific genes underlying major QTLs regulating two most vital drought-responsive yield traits delineated in this study from the comparative C. microphyllum and ICC 4958 transcriptome profiles can be deployed for various functional and translational genomic (marker-assisted breeding and transgenics) research. This will provide multiple essential leads, which could accelerate genomics-assisted crop improvement to develop drought tolerant superior high-yielding cultivars of chickpea.
The current study utilized a RNA-seq strategy to sequence transcriptome of an Indian origin perennial wild C. microphyllum accession for the first time at a high-resolution scale. The optimized de novo high-quality transcriptome assembly and comprehensive transcript sequence data for root and shoot tissues of this wild accession were generated to provide a global view of its gene content. The gene-encoding transcripts differentially expressed between roots and shoots of C. microphyllum accession were correlated with that of an available whole genome transcript sequence data of a desi chickpea accession (ICC 4958) based on comparative transcriptome profiling to identify candidate genes, transcription factors (TFs) and key regulators/metabolic pathways especially governing inherent root-mediated natural adaptation characteristics of C. microphyllum in adverse agro-climatic condition. Based on these, numerous experimentally well-validated C. microphyllum transcript/gene-based SSR and SNP markers derived from the coding and regulatory sequence components of genome and/or root-and shoot-specific abiotic (drought, cold and salinity) stress-responsive differentially expressed genes of ICC 4958 were discovered to accelerate high-throughput genetic analysis in chickpea. These informative genic SSR and SNP

*QTLs
LGs/ chromosomes  markers were genotyped successfully in an advanced generation mapping population to construct a high-density inter-specific genetic linkage map (transcript map) and for molecular mapping of major QTLs governing two drought-responsive important yield traits (yield per plant and harvest-index) in chickpea. An integrated genomic approach by combining high-resolution drought yield QTL mapping with comparative transcriptome profiling was employed to delineate candidate gene(s) harbouring major QTLs regulating yield per plant and harvest-index in chickpea. Overall, this study provides a deep transcriptomic insight into the gene regulatory networks governing differential root-mediated natural adaptation responses in a hardy wild perennial C. microphyllum accession as compared to that of an already existing cultivated high-yielding drought stress tolerant desi chickpea accession (ICC 4958). The comprehensive transcriptome, informative genic SSR and SNP markers, differentially expressed genes, a high-resolution inter-specific genetic linkage/transcript map and major drought yield QTLs and potential genes/natural allelic variants especially generated in our study at a global scale can serve as a useful resource to drive translational genomics for genetic enhancement in chickpea. High-quality transcriptome assembly. For de novo transcriptome assembly, the high-quality filtered sequence reads obtained from three independent biological replicates of root and shoot tissue samples of a C. microphyllum accession were analysed with Trinity (vr2012-05-18) and CLC Genomics Workbench (v4.7.2) tools. To obtain best transcriptome assembly, most vital parameters including k-mer and insert length, and expected coverage were optimized as per Garg et al. 19,42 . Reference-based transcriptome assembly of C. microphyllum using the available draft genome sequence of its evolutionary close desi chickpea accession (ICC 4958) (CGAP v2.0 50 ) as a reference was performed following the criteria as defined by Garg et al. 42 . The diverse quality parameters of developed transcriptome assemblies were evaluated using the perl scripts employed in NGS QC Toolkit. The high-quality transcripts generated for C. microphyllum were functionally annotated by their BLASTX search against desi chickpea proteome 50 , NCBI non-redundant (nr) protein sequence database (http://www. ncbi.nlm.nih.gov/refseq/about/nonredundantproteins) and PFAM database v27.0 (http://pfam.sanger.ac.uk). The KOG (eukaryotic orthologous groups of proteins, ftp://ftp.ncbi.nih.gov/pub/COG/KOG)-based functional annotation of the C. microphyllum transcripts and identification of TF-encoding genes from these transcripts were performed following the methods of Kujur et al. 51 with ≥ 50 and 70% query coverage and percent identity, respectively.

Methods
Differential gene expression analysis. For gene-encoding transcript expression analysis, all the high-quality sequence reads obtained from each of the root and shoot tissue samples of C. microphyllum accession were mapped/aligned individually on the constituted de novo transcriptome assembly using Tophat (v2.0.0) 73 and CLC Genomics Workbench tools (http://www.clcbio.com) with default parameters. The Cufflinks (v2.0.2) 74 , Cuffmerge and Cuffdiff were used to estimate the differential expression of genes/transcripts in roots and shoots of C. microphyllum accession. The genes/transcripts revealing at least two-fold expression change with P value ≤ 0.05, were considered significant to be differentially expressed. For comparative transcriptome profiling, the raw transcriptome sequence data of control (unstressed) as well as drought, cold and salinity stress-imposed root and shoot tissue samples of ICC 4958 that are freely accessible at NCBI Gene Expression Omnibus (GEO) database (GSM1299261-GSM1299268) were retrieved 38 . The processing of these raw sequence reads to obtain high-quality reads, reads mapping to C. microphyllum transcript assembly and differential expression profiling of root and shoot (control/unstressed)-derived transcripts/genes between C. microphyllum and ICC 4958 were performed following the above-mentioned strategy. The heatmaps depicting the differential expression profiles of genes/transcripts (log2 fold-change) in unstressed control roots and shoots of C. microphyllum and ICC 4958 were constructed based on hierarchical clustering employing MultiExperiment Viewer (Mev v4.9). The KOG-based functional annotation of the C. microphyllum and ICC 4958 transcripts and TF-encoding gene identification from the differentially expressed transcripts were performed following the aforesaid strategy. The molecular interactions based on biochemical pathway mapping of differentially expressed genes/transcripts were performed using KEGG pathway database (http://www.genome.jp/kegg/pathway.html).

Experimental validation of gene expression. For experimental validation of gene expression through
semi-quantitative and quantitative real-time (RT) PCR assays, the genes/transcripts exhibiting high root-specific differential expression especially in C. microphyllum and ICC 4958 were selected to design primers using Primer Scientific RepoRts | 6:33616 | DOI: 10.1038/srep33616 Express (v3.0) (Applied Biosystems, USA). The RNA isolated from the root and shoot tissue samples of six drought tolerant and sensitive desi and kabuli chickpea accessions including a C. microphyllum accession (used for transcriptome sequencing) following aforesaid methods. In the RT-PCR assay, RNA isolated from three independent biological replicates of each root and shoot sample and two technical replicates of each biological replicate with no template and primer as control was used for expression study. For gene/transcript expression profiling, the isolated high-quality RNA was amplified with gene-specific primers using ABI 7500 Real-Time PCR System (Applied Biosystems, USA), following the methods as previously described by Kujur et al. 52 and Bajaj et al. 75 . An internal control gene elongation factor 1-alpha (EF1α) was utilized in RT-PCR assay for normalization of gene/transcript expression value. The correlation between transcriptome sequencing and RT-PCR profiles assayed in roots and shoots of C. microphyllum and ICC 4958 was analysed through R program. The differential expression (fold change) profiles of genes/transcripts assayed in roots and shoots of drought tolerant and sensitive chickpea accessions was measured. Significant difference in gene/transcript expression among diverse tissues and/or accessions was determined by LSD-ANOVA significance test using SPSS 17.0 (http://www.spss.com/statistics). The significant difference in gene/transcript expression profile was visualized by a heat map using MeV.

Development and annotation of genic SSR and SNP markers.
To discover genic SSRs, the C. microphyllum transcripts were analysed with MISA (MIcroSAtellite; http://pgrc.ipk-gatersleben.de/misa/) as per Kujur et al. 52 . The in silico polymorphic genic SSRs based on variation of repeat-units length in the transcript/gene sequences between C. microphyllum and ICC 4958 (desi reference genome) were identified following the methods as described previously 27,76,77 . To develop genic SSR markers, the primer-pairs were designed from the sequences flanking the SSR repeat-motifs including the in silico polymorphic SSR repeats using the Primer3 interface module of MISA as per Kujur et al. 52 .
The genic SNPs in the transcript/gene sequences between C. microphyllum and ICC 4958 were mined using Tophat and samtool. The high-quality genic SNPs differentiating the C. microphyllum and ICC 4958 accessions were screened following the criteria as defined by Jain et al. 78 and Kujur et al. 51 . These mined SSR (in silico polymorphic SSR) and SNP markers were structurally and functionally annotated in the diverse coding and non-coding regulatory sequence components of desi chickpea (ICC 4958) genome as per Jain et al. 78 and Kujur et al. 42,51 . The functional significance of SSRs and SNP markers was determined based on their presence particularly in the coding (including non-synonymous and large-effect substitutions) and regulatory regions of genes/TF-encoding genes of C. microphyllum that were differentially expressed in the roots and shoots of ICC 4958 under drought, cold and salinity stresses.

Validation of genic SSR markers. For experimental validation of genic SSR markers, 192 C. microphyllum
transcript-derived SSR markers and 96 in silico polymorphic SSR markers (physically mapped on chromosomes) were selected based on their localization especially in the diverse coding and regulatory sequence components of root-specific genes (TFs)/transcripts of C. microphyllum that were differentially expressed especially in the roots of ICC 4958 under drought stress. These SSR markers were genotyped in 96 wild and cultivated chickpea accessions 2,11 (including known drought tolerant and sensitive accessions selected for differential expression profiling) following the methods as described previously 27,76,77 . Validation of genic SNPs. To validate genic SNPs identified between C. microphyllum and ICC 4958 in silico, the genotyping information of SNPs differentiating these two chickpea accessions was compared/correlated with the available in silico SNP database according to their nature/types of SNP alleles and physical positions (bp) on the desi (ICC 4958) genome 11,25,26,50,51 . For experimental validation and understanding the functional significance of mined SNPs, 192 coding (including non-synonymous and large-effect substitutions) and regulatory SNPs (physically mapped on chromosomes) were selected from root-specific genes/TF-encoding genes of C. microphyllum based on their drought-responsive differential expression especially in the roots of ICC 4958. These SNPs were further genotyped in 96 wild and cultivated chickpea accessions (including known drought tolerant and sensitive accessions selected for differential expression profiling) using Sequenom MALDI-TOF (matrix-assisted laser desorption ionization-time of flight) MassARRAY as per Saxena et al. 2 and Bajaj et al. 11 .
Construction of an inter-specific genetic linkage map and QTL mapping. A 190 F 8 RIL (recombinant inbred line) mapping population derived from the inter-specific crosses between two parental chickpea accessions (C. arietinum desi accession ICC 4958 × C. reticulatum wild accession ICC 17163) was developed to construct a genetic linkage map (transcript map) and for identification/mapping of major drought-responsive yield QTLs in chickpea. The genotyping data of drought-responsive root-specific differentially expressed gene-derived SSR (genotyped by gel-based assay) and SNP (Sequenom MALDI-TOF MassARRAY) markers revealing polymorphism between mapping parental accessions (ICC 4958 and ICC 17163) was utilized to generate an inter-specific genetic linkage map (transcript map) at a higher LOD (logarithm of odds) threshold (> 4.0) with Kosambi mapping function following Saxena et al. 17 and Bajaj et al. 79 . Accordingly, the informative polymorphic genic markers were integrated into eight LGs of a transcript map based on their centiMorgan (cM) genetic distance. The LGs were designated (LG1 to LG8) as per their corresponding marker physical positions (bp) on the chromosomes of desi chickpea genome. The previously documented genic SNP markers mapped on eight LGs exhibiting polymorphism between two mapping parental accessions (ICC 4958 and ICC 17163) were used as anchor markers to define the identity of LGs/chromosomes of an inter-specific genetic map constructed in the present study. The identity of LGs was further ascertained by positions (cM) and groupings shared by anchor markers among eight LGs/chromosomes between current and past studies 25,26 .
The RIL mapping individuals and two parental accessions (ICC 4958 and ICC 17163) were grown in the field [as per randomized block design (RBD) with at least two replications] during crop growing season (minimum 13 ± 2.2 °C and maximum 21 ± 2.5 °C) for two consecutive years (2012 and 2013) under rainfed and irrigated environments at two diverse geographical locations (New Delhi: latitude 28.6 °C and longitude 77.2 °C, Palampur: 32.1 °C and 76.5 °C) of India. Accordingly, these were phenotyped for two vital drought tolerance-related traits, namely yield (g) per plant (YP) and harvest-index (HI%) especially under rainfed and irrigated environments following Varshney et al. 40 . The mean, standard deviation, coefficient of variation (CV), frequency distribution, Pearson's correlation coefficient and broad-sense heritability (H 2 ) of YP and HI traits phenotyped in a RIL mapping population were measured as per Saxena et al. 17 and Bajaj et al. 79 . The QTL mapping was performed by integrating the genotyping data of genic SSR and SNP markers genetically mapped on eight LGs/chromosomes with the field phenotypic data (YP and HI) of 190 RIL mapping individuals and parental accessions based on composite interval mapping function of MapQTL v6.0. The LOD threshold score of more than 4.0 at 1000 permutation was considered significant (p < 0.05) to identify and map the major QTLs on LGs/chromosomes governing YP and HI traits in chickpea. The positional genetic effect and PVE% defined by QTLs were evaluated at a significant LOD as per Saxena et al. 17 and Bajaj et al. 79 .