Molecular Mapping of QTLs for Heat Tolerance in Chickpea

Chickpea (Cicer arietinum L.), a cool-season legume, is increasingly affected by heat-stress at reproductive stage due to changes in global climatic conditions and cropping systems. Identifying quantitative trait loci (QTLs) for heat tolerance may facilitate breeding for heat tolerant varieties. The present study was aimed at identifying QTLs associated with heat tolerance in chickpea using 292 F8-9 recombinant inbred lines (RILs) developed from the cross ICC 4567 (heat sensitive) × ICC 15614 (heat tolerant). Phenotyping of RILs was undertaken for two heat-stress (late sown) and one non-stress (normal sown) environments. A genetic map spanning 529.11 cM and comprising 271 genotyping by sequencing (GBS) based single nucleotide polymorphism (SNP) markers was constructed. Composite interval mapping (CIM) analysis revealed two consistent genomic regions harbouring four QTLs each on CaLG05 and CaLG06. Four major QTLs for number of filled pods per plot (FPod), total number of seeds per plot (TS), grain yield per plot (GY) and % pod setting (%PodSet), located in the CaLG05 genomic region, were found to have cumulative phenotypic variation of above 50%. Nineteen pairs of epistatic QTLs showed significant epistatic effect, and non-significant QTL × environment interaction effect, except for harvest index (HI) and biomass (BM). A total of 25 putative candidate genes for heat-stress were identified in the two major genomic regions. This is the first report on QTLs for heat-stress response in chickpea. The markers linked to the above mentioned four major QTLs can facilitate marker-assisted breeding for heat tolerance in chickpea.


Introduction
In recent years, the adverse impact of climate change on agriculture is well recognized all over the globe. The ever-increasing day and night temperature is going to affect the production of crops, especially those grown in the winter [1]. In this context, heat-stress due to rise in temperatures remains a challenge in developing crop varieties that are adaptive to changing climatic conditions.
Chickpea is a nutrient-rich grain legume crop cultivated in arid and semi-arid regions. The chickpea grain is an excellent source of proteins along with a wide range of essential amino acids and vitamins. In the fight against hidden hunger all over the globe, the role of legumes such as chickpea is indispensable. Grown in over 60 countries and traded in over 190 countries, chickpea is the second most consumed pulse crop in the world after common bean [2]. Due to global warming, several noticeable changes occurred in the cropping system and intensity in the recent past. These are delaying the cultivation of chickpea to relatively hot conditions [1]. Generally, the crop faces heat-stress during reproductive phase under late sown condition in the tropical and semi-arid regions [3]. Reports state that the exposure to temperature, 35 • C and above, even for a few days, during reproductive phase has a negative impact on optimum yield in chickpea [4,5]. Unlike drought and other abiotic stresses, until recently, the importance of breeding for heat-stress conditions in chickpea has not been realized [1].
Grain yield under heat-stress is considered to be one of the important criteria for assessing heat tolerance in chickpea [3][4][5]. However, chickpea yield is known to be highly influenced by environments [6]. Due to genotype by environment (G × E) interaction, breeding for heat tolerance through conventional breeding approaches based on yield parameter sometimes limits selection for heat-stress tolerance in chickpea.
In recent years, progress has been made in genomics-enabled trait dissection in several crop plants, including chickpea. Several studies have been carried out earlier to identify the quantitative trait loci (QTLs) for tolerance to various biotic stresses [7,8], and abiotic stresses like drought tolerance [9], and salinity tolerance [10][11][12] in chickpea. Moreover, genomic regions associated with heat tolerance have been reported in several crops, including wheat, rice, maize, barley, potato, tomato, cowpea, azuki bean, brassica [13]. Pod setting (seed set) and grain yield have been used as proxy traits to detect QTLs for heat tolerance in different crops [14][15][16][17][18]. Similarly, in chickpea, the number of filled pods, total number of seeds, biomass, and harvest index were found to be significantly associated with heat tolerance [3,19]. However, to date, QTLs for heat tolerance have not been reported in chickpea.
In this study, genotyping by sequencing (GBS)-based single nucleotide polymorphism markers were used to identify key genomic regions responsible for heat tolerance. In addition, putative candidate genes for heat tolerance in these genomic regions were identified using the available chickpea genome sequence information [20].

Response of Parents and Recombinant Inbred Lines (RILs) under Heat-Stress and Non-Stress Environments
The descriptive analysis of parents and RILs are presented in Table 1. Predicted means for all the traits in parents differed significantly in both heat-stress environments, except biomass per plot (BM). In the non-stress environment, predicted means for grain yield per plot (GY), BM, harvest index (HI) and %PodSet were non-significant between parents, while filled pods per plot (FPod) and total number of seeds per plot (TS) were significant. The range of variation in all the traits was high in stress environments ( Table 1). The combined analysis of variance (ANOVA) for both the stress environments revealed that significant variation existed in RILs for all the traits measured, except BM, whereas under non-stress environment relatively low genetic variability was observed. Transgressive segregants in both directions were observed for several traits in the RIL population (Figure 1a,b).
The potential use of a trait in a breeding program relies on the heritability of that trait. Under both the heat-stress conditions, the heritability of all the traits was high (72.0-90.7%), except BM in summer 2014 (49.8%). Whereas, under non-stress environment the heritability of the traits was moderate (47.6-66.0%) ( Table 1).

Relationship between Yield and Yield Determining Traits
Heat tolerance is a complex trait and can be estimated indirectly through yield and yield contributing traits under heat-stress. All the traits-visual score (VS), FPod, TS, BM and %PodSet were positively correlated with yield (r = 0.51 **-0.90 **) under both the heat-stress environments and pooled over analysis except HI (r = 0.32 **) under heat-stress environment of 2013 (Table 2). In addition, VS had positive association with FPod (r = 0.68 **-0.80 **) and TS (r = 0.67 **-0.79 **). Likewise, %PodSet was found to have a strong positive correlation with FPod (r = 0.59 **-0.77 **) and TS (r = 0.60 **-0.78 **) under both the heat-stress environments as well as in pooled analysis ( Table 2). In contrast, under non-stress environment, the correlation with yield was low for %PodSet (r = 0.17 **) and HI (r = 0.33 **), and high for other traits (r = 0.63 **-0.91 **) ( Table 2). Regression analysis between the traits and yield revealed that all the traits exhibited medium to high variation for yield (25% to 81%) in both stress environments as well as pooled over years ( Figure S3a-c). In non-stress environment, %PodSet had low contribution (3%) whereas BM was found to have high variation for yield (82%) ( Figure S3d). A significant correlation between the yield and yield contributing traits under heat-stress environment indicated that these traits can be used in direct or indirect selection for improving heat tolerance in chickpea.

Sequencing Data and SNP Discovery
The parents of the mapping population (ICC 4567 × ICC 15614) were sequenced at higher depth (5× coverage), and a total of 19.63 million reads containing 1.70 Gb for ICC 4567, and 15.79 million reads containing 1.37 Gb for ICC 15614, were generated. In addition, 3333.41 million reads containing 289.70 Gb were generated from 292 RILs. The number of reads generated varied from 6.86 million (RIL099) to 20.66 million (RIL112) with an average of 11.42 million per line. The single nucleotide polymorphisms (SNPs), identified using the software SOAP, were analyzed to remove heterozygous SNPs in the parents, and a set of 396 SNPs were identified across 292 RILs. The sequence details of all SNPs have been provided in Table S1a,b.

Genetic Linkage Map and Marker Distribution
The 396 polymorphic SNPs obtained from GBS were used for genetic map construction. The genetic linkage map covered 529.11 cM of the chickpea genome with an average interval of 1.95 cM between markers (Table S2 and Figure S1). The highest number of markers was in CaLG04 (57), while the lowest number of markers was in CaLG08 (10) ( Figure S1). CaLG08 showed the highest marker density with 1.78 markers per cM on average. The lowest marker density was observed for CaLG02, which had 0.29 markers per cM on average. Overall, the map had on average 0.51 markers per cM (Table S2).
All the major QTLs present in the genomic region of CaLG05 were found to exist in the pooled analysis for the two stress environments (Table 3). In CaLG05, two major QTLs for VS and HI were found explaining 15.1% (LOD 11.1) and 18.5% (LOD 13.0) of phenotypic variation, under the heat-stress environment (2014), respectively (Table S3). In contrast, during the stress environment in 2013, one major QTL for VS was found close to the genomic region on CaLG05 with a phenotypic variation of 13.88% (LOD 12.05) (Table S3). Through single marker analysis (SMA), Ca5_44667768 was co-segregated with the four major QTLs in this genomic region.

Genomic Region on CaLG06
A second genomic region, harbouring QTLs for four important traits in this study, was identified having the marker interval Ca6_14353624-Ca6_7846335 (Table 3 and Figure 2b). The QTLs for FPod (qfpod03_6), GY (qgy03_6), %PodSet (q% podset08_6), and VS (qvs05_6) spanned a genetic length of 19.14 cM (~6.50 Mb on physical map) in CaLG06. The range of phenotypic variation shown by various traits in this genomic region was from 3.92 to 11.07% (Table 3).

QTLs Identified on Other LGs
In the present work, a total of 13 QTLs were identified consistently across two heat-stress environments showing both major and minor effects for various traits measured. Apart from the QTLs identified in CaLG05 and CaLG06, a QTL for GY (qgy01_1) was found in the same position (40.0 cM) demonstrating 7.33% and 10% of phenotypic variation in the first and second year, respectively, on CaLG01 (Table S4).

Mapping of Epistatic QTLs (E-QTLs)
Epistatic interaction analysis revealed that 19 QTL pairs were involved in the epistatic interactions covering seven LGs (Table 4). A significant effect was observed for all the epistatic interactions. However, no significant interaction between epistasis and environment was observed, except for the trait biomass (BM).
Two epistatic QTL pairs for VS were found to have loci distributed on four different LGs accounting for 3.43% phenotypic variation. In the case of FPod, two QTLs were found to be interacting in the same LG, CaLG02. Another QTL pair was found for FPod to interact with each other in two different LGs (Table 4). These two epistatic QTL pairs for FPod together explained a phenotypic variation of 2.94%.
The highest number of epistatic QTL pairs (nine pairs) were detected for TS in this population and have contributed up to 12.38%. The epistatic interaction for TS was found in all the linkage groups, except CaLG03 and CaLG07. One QTL interaction pair was detected for GY interacting from CaLG01 to the locus on CaLG02 with a phenotypic variation 0.83% (Table 4 and Figure S2). Similarly, in the case of %PodSet, four epistatic QTL pairs were found to interact with each other in three linkage groups CaLG01, CaLG03, and CaLG04 showing a phenotypic variation of 5.79%.
In addition, an interaction between non-QTL, and additive and additive × environment-QTL was found in the case of BM, which showed 1.22% phenotypic variation. Concurrently, five loci (loci located at 10.1 cM and 26.4 cM in CaLG01, 2.2 cM and 75.6 cM in CaLG04, and at 44.5 cM in CaLG05) were observed to have interaction simultaneously with several other loci affecting the expression of the particular trait. Two loci (eqts2_1/eqpodset2_1 in CaLG01 and neqfpod4_5/neqts9_5 in CaLG05) controlling two or three different traits were also interacted with other loci (Table 4).

Phenotypic Evaluation of RILs and Parents in Field Condition
Sowing during the month of February proved to be an ideal condition to expose chickpea crop to heat-stress and selecting heat tolerance lines in earlier studies under field conditions at ICRISAT, Patancheru, India [19,21]. A recent study on chickpea reported 34 • C as the threshold temperature for pod setting and also observed that at 35 • C, pod set was reduced by 50% in chickpea genotypes [19]. The average maximum temperatures (37.5 • C and 36.7 • C in summer 2013 and summer 2014, respectively) in both the heat-stress environments found were ideal for phenotyping RIL population. An average maximum temperature of 29.4 • C was recorded in non-stress environment, which was considered as control for this study. This temperature was ideal for sowing in the non-stress environment for the timely sown crop [22].
The frequency distribution of measured traits showed the characteristics of continuous variation (Figure 1a,b). Paliwal et al. (2012) [23] in RILs of wheat and Buu et al. (2014) [24] in BC 2 F 2 population in rice, reported several transgressive segregants for heat tolerance. Similarly, in this present study, transgressive segregants in both directions were observed, indicating that both parents have contributed alleles for heat tolerance in the RILs (Figure 1a,b). A significant variation found among the RILs for all the traits indicate the presence of genetic diversity in the selected parents for the selected traits under heat-stress condition. Parents differed significantly for all the traits in both the heat-stress environments, except biomass (BM).
High heritability (H 2 ) values were observed for all the traits measured under both the heat-stress environments, except for biomass in summer 2014, which indicates that there is a high probability of achieving the same kind of results if the trial is repeated under similar growing conditions. Yield under high temperatures is the key objective for heat tolerance breeding in chickpea. Traits such as FPod, %PodSet and TS contributing to increased yield under high-temperature stress can be treated as a proxy for heat tolerance. The presence of significant correlations between yield and other traits in heat-stress environments indicated that these traits can be used as selection criteria for heat tolerance.
FPod and TS had a strong correlation with yield (88 to 90%) under both the stress environments. Such high correlation of these traits toward yield was reported earlier in chickpea under abiotic stress [10,11]. In addition, VS and %PodSet was also found to have good correlation (50 to 79%) with yield. However, BM and HI showed large difference in correlation with yield in both the heat-stress conditions. Positive and strong association of the four traits-FPod, TS, VS and %PodSet with grain yield revealed the importance of these traits in determining yield under heat-stress environment. Hence, detecting QTLs of these traits under stress would be helpful in heat tolerance programme.

QTL Mapping for Heat Tolerance
The genomic region in CaLG05 harbours QTLs for FPod, TS, GY, and %PodSet, which were reportedly associated with heat tolerance in chickpea [3,19]. Interestingly, the positions of the QTLs (qts02_5, qgy02_5, q% podset06_5) for TS, GY, and %PodSet were identified in the same position over the years, which strongly confirm the QTLs in these positions.
The presence of four major co-localized QTLs (qfpod02_5, qts02_5, qgy02_5, and q% podset06_5) suggests tight linkage or the phenomenon of pleiotropy and the phenotypic correlations between these traits were highly significant in both the stress environments. Moreover, the tolerant parent ICC 15614 is contributing the desirable alleles for all the QTLs found in the two genomic regions in CaLG05 and CaLG06.
Identification of QTLs at the same positions in both the heat-stress environments indicate their possible practical utility in breeding for heat-stress tolerance in subsequent studies [25]. Several co-localized QTLs for various traits were found which could possibly due to pleiotropy or tightly linked QTLs. Fine mapping of the target genomic region will further help in resolving the issues of pleiotropy and tight linkage. The incorporation of a higher number of markers into the existing genetic map can further narrow down the genomic regions identified.
QTLs for traits such as FPod, TS, and GY were not expressed under non-stress condition, confirming the fact that these QTLs were only expressed under high-temperature condition. Two major QTLs for HI were identified in CaLG01 and CaLG04 explaining the phenotypic variation of 12.03% (LOD 8.8) and 12.53% (LOD 7.9), respectively. In addition, three minor QTLs including one for HI and two for %PodSet were found in different LGs. The fewer number of detected QTLs and their unique positions in the non-stress environment is a strong evidence that there is no correspondence between QTLs found in non-stress with the QTLs found in heat-stress environment. This phenomenon proves the fact that those QTLs identified in heat-stress condition were independent and exclusive for heat tolerance.

Epistatic QTLs for Heat Tolerance
Epistatic interaction is one of the key factors controlling the expression of a complex trait. The epistatic interaction analysis of QTLs provides a more comprehensive knowledge of the QTLs and their genetic behaviour underlying the trait [26,27].
In the current study, 19 pairs of digenic epistatic QTLs were found to be associated with the six traits: VS, FPod, TS, GY, BM, and %PodSet. Maximum number epistatic QTLs loci were observed for TS (nine), followed by %PodSet (four). In this study, some loci such as eqts2_1/eqpodset2_1, eqts2_1/eqpodset2_1, eqpodset2_1/eqts2_1, neqts9_5/neqfpod4_5, neqfpod4_5/neqts9_5 were simultaneously controlling more than one trait indicating the pleiotropy nature of the traits.
Four categories of epistatic interaction were found in this study such as, additive × additive, additive × non-QTL, non-QTL × non-QTL, and additive × (additive-environment) × non-QTL interaction. FPod and VS showed two epistatic interactions each. Out of two epistatic interactions, one additive × additive epistatic interaction was found for both FPod and VS.
For GY, one additive × additive QTL epistatic interaction was found. For TS, five additive × additive QTL epistatic interactions, three non-QTL × non-QTL interaction and one additive × non-QTL interactions were observed. Similarly, two additive × additive QTL interactions, one non-QTL × non-QTL interaction and one additive × non-QTL interaction were observed for %PodSet. All the epistatic interactions were found to be significant.
The additive effects were found in both directions for all the traits. Nine interactions had negative additive effects, meaning that recombinant allele combinations could increase the particular trait value. Similarly, ten epistatic QTL interactions having positive additive effects, indicating parental allele combinations, would help to improve the trait [28].
Presence of epistatic interactions for a given trait will make the selection difficult. Interestingly, all major QTLs had no epistatic interaction and this will increase the heritability of the trait and make the selection easy.

Putative Candidate Genes for Heat Tolerance
Recent progress in functional genomics facilitates the elucidation of the important role of candidate genes for expression of tolerance against abiotic stress in plants [29][30][31]. In the present study, mining of the candidate genes for heat tolerance revealed 236 genes in 2.28 Mb (44.6-46.9 Mb) region in CaLG05 and 550 genes in 6.50 Mb (7.85-14.35 Mb) in CaLG06 (Tables S5 and S6). Based on functional categorization, many genes were found to be associated with biological processes (168 genes in CaLG05 and 365 genes in CaLG06) in the two genomic regions.
Gene ontology classification revealed a total of 25 putative candidate genes (11 in CaLG05 and 14 in CaLG06) known to function, directly or indirectly, as heat-stress response genes in several plant species (Table S9a,b). Of the 25 candidate genes, five genes encode protein like farnesylated protein 6 (AtFP6), ethylene-responsive transcription factor ERF114, ethylene-responsive transcription factor CRF4, F-box protein SKP2B, and ethylene-responsive transcription factor RAP2-11. These genes were identified to have key roles in heat acclimation and growth of plants under severe heat-stress condition. Many transcription factors, enzyme, and stress responsive element binding factors responsible for heat tolerance in various plant species were reported earlier [32]. Furthermore, various heat shock proteins (HSPs), ethylene forming enzymes (EFEs), and ethylene-responsive element factors (ERFs) were found to be candidate genes for heat tolerance in soybean and cowpea, two of the plant species closest to chickpea [32].
The role of various heat shock proteins and heat-stress transcription factors has been widely accepted and reported in different crops [33]. The role of HSP90 transcription factors under heat-stress conditions was also reported in chickpea [34]. Five putative genes were identified in the two examined genomic regions, encoding for either heat shock proteins or heat shock transcription factors contributing for thermo-tolerance.
Oxidative stress can occur in parallel with heat-stress through the formation of reactive oxygen species (ROS) [35]. Three putative candidate genes were also observed in this study to have a role in defying oxidative stress and recovering plants from heat-stress damage. These genes encode different types of proteins like protein tansparent testa glabra 1, peroxidase 52, and zinc finger protein CONSTANS-LIKE 5. In addition, certain signalling molecules like ethylene, abscisic acid (ABA), and salicylic acid are among a few to have a significant role in the development of heat tolerance [36]. In this study, a few genes-MYB44, AKH3, and RAN1-were found to involve with these signalling molecules through upregulation process to mitigate the heat-stress. Being a preliminary study, evaluation of these putative candidate gene-functions in chickpea through fine mapping and gene expression study is necessary to use them for further study. In all the environments, the field was solarized using polythene mulch during the preceding summer to sanitize the field, especially to avoid incidence of root diseases. Sowing was done on the ridges using ridge and furrow method with inter-and intra-row spacing of 60 × 10 cm. Each plot consisted of a 2 m long row. Need-based insecticide sprays were provided to control pod borer (Helicoverpa armigera) and the experimental plots were kept weed-free through manual weeding. Before sowing, seeds were treated with the mixture of fungicides 0.5% Benlate ® (E.I. DuPont India Ltd., Gurgaon, India) + Thiram ® (Sudhama Chemicals Pvt., Ltd., Gujarat, India).

Plant Material and Treatment Condition
The experimental design was laid out in a 15 × 20 alpha lattice design with three replications. The sowing for the non-stress environment was done on the residual moisture in the last week of November 2013 and provided with essential irrigation. The planting was done in the first week of February for stress environments to expose the reproductive phase of RILs to heat-stress (>35 • C). The stress experiments were provided with irrigation to avoid the confounding effect of moisture stress during the heat screening.
In chickpea, a temperature higher than 35 • C during reproductive phase adversely affects growth, development, and yield [1,19]. The parents used for developing RIL population for this study showed significant variations at this temperature (35 • C and above) in an earlier study [19]

Variables Measured
Number of filled pods per plot (FPod), total number of seeds per plot (TS), grain yield per plot (GY, g), harvest index (HI, %), biomass (BM, g) and percent pod setting (%PodSet), were reportedly found to be associated with heat tolerance in chickpea [3,19]. These six traits along with visual score on podding behaviour (VS) were recorded in the RIL population. The data for FPod, TS, GY, BM, and HI were recorded from a half-meter (0.5 m) long continuous patch out of the 2-m plot. VS at maturity and %PodSet were recorded from the entire plot. For visual scoring, score-1 was considered most sensitive (least number of pod-bearing ability), whereas, score-5 was taken as the most tolerant (maximum number of pod-bearing ability) under heat-stress. In the non-stress environment, all RILs were assumed to behave more or less the same. Hence, no visual score data were recorded in this environment.

DNA Extraction, Genotyping, and SNP Calling
DNA from 292 RILs, along with the parents, was isolated from 15-day old seedlings following the high-throughput mini-DNA extraction method [37]. Genotyping was done using GBS approach [38]. The GBS libraries from the parental lines and RILs were prepared using ApeKI endonuclease (recognition site: G/CWCG) and were sequenced using the Illumina HiSeq 2000 platform (Illumina Inc, San Diego, CA, USA). The detailed procedure of genotyping approach was described by Jaganathan et al. (2015) [25].
For SNP calling the raw reads obtained were first de-bimultiplexed using sample barcodes, and adapter sequences were removed using a custom Perl script ( Figure S5). The reads having more than 50% of low-quality base pairs (Phred < 5%) were discarded and filtered data were used for calling SNPs after due quality check (Q score > 20). The high-quality data from each sample were aligned to the draft genome sequence (CaGAv1.0) of chickpea [20] using SOAP [39]. After SNP calling, the polymorphic loci were determined by following the criteria defined in [25].

Linkage Map Construction, QTL Detection and Mining of Candidate Genes
By adopting a stringent selection criterion including the missing percentage, minor allele frequency, and percent heterozygosity, the final number of SNPs included in the analysis were 396. The selected panel of robust SNPs were used for construction of genetic maps.
A linkage map was constructed with the 396 SNPs using JoinMap 4.1 [40]. Composite interval mapping in QTL Cartographer-V 2.5 [41] was employed to identify the QTLs responsible for heat tolerance with a forward and backward stepwise regression (threshold p-value < 0.05). A window size of 10 cM, along with a walking speed of 1.0 cM, and 1000 permutations for p < 0.05 were chosen for the QTL analysis. QTL × QTL and QTL × E interactions were estimated using the QTL Network version 2.0 (http://ibi.zju.edu.cn/software/qtlnetwork/) which is based on a mixed linear model.
First-dimensional genome scan (with the option to map epistasis) and second-dimensional genome scan (to detect epistatic interactions with or without single-locus effect) were applied. A significance level of 0.05 with 1000 permutations, 1.0 cM walk speed, 10.0 cM testing window and filtration window size were employed for the epistatic QTL analysis. QTL was named with prefix "q" for main-effect QTL, "eq" for epistatic QTL and "neq" for non-QTL epistasis followed by the abbreviated trait name and the identity of the linkage group involved.
The identified markers along with the flanking sequences were mapped on the chickpea reference genome CaGAv1.0 [20]. The genes present within the physical locations of these markers were extracted from the genome features file and were searched against TrEMBL and Swiss-Prot databases. Further functional annotation was done using UniProtKB. The Gene Ontology annotations were categorized into three categories: biological processes (BP), molecular function (MF) and cellular components (CC).

Analysis of Variance, Predicted Means (BLUP), Heritability, and Correlations
The analysis of variance (ANOVA) for the RIL population was performed using GenStat (17th Edition), for individual environments using mixed model analysis. For each trait and environment, the analysis was performed considering entry and block (nested within replication) as random effects and replication as fixed effect.
To pool the data across environments, and to make the error variances homogeneous, individual variances were estimated and modelled for the error distribution using residual maximum likelihood (ReML) procedure. Z value and F value were calculated for random effects and fixed effects, respectively. For single and multi-environment, QTL mapping was performed using predicted means (BLUP-Best Linear Unbiased Prediction) [42]. Broad

Conclusions
The present study identified two potential genomic regions harbouring major QTLs for several heat responsive traits that are directly related to heat tolerance in chickpea. The two regions consistently appeared at the same map position across two years. Epistatic effects were not observed for major QTLs and no QTL × E interaction in the CaLG05 region. The results laid a foundation in understanding heat tolerance and increases the confidence of breeders to proceed with early generation selection for heat tolerance through marker-assisted breeding. In addition, the candidate genes identified in the two genomic regions further help to understand the mechanism of heat tolerance.
Supplementary Materials: The following are available online at http://www.mdpi.com/1422-0067/19/8/2166/ s1. Figure S1: Intra-specific genetic map of chickpea RIL population (ICC 4567 × ICC 15614) with 271 GBS-based SNPs covering 529.11 cM. Genetic distances (cM) were shown on the left side and the markers were shown on the right side of the bars. The map was constructed using JoinMap 4.1 and Kosambi function, Figure S2: The epistatic QTLs on linkage groups detected by QTLNetwork v 2.0 in the RIL population (ICC 4567 × ICC 15614). Lines joining two QTLs represents the epistatic interaction between them, Figure S3 Figure S5: Pipeline of Bioinformatics analysis: GBS data processing and SNP calling, Table S1: (a) Summary sequence data generated genotyping-on 292 RILs and two parents (ICC 4567 and ICC 15614) using GBS approach; (b) Summary of called SNPs on 292 RILs and two parents (ICC 4567 and ICC 15614) using GBS approach, Table S2: Features of intra-specific genetic map developed using 271 SNPs and RIL population ICC 4567 × ICC 15614, Table S3: Summary of QTLs identified in two heat-stress environments, pooled environments and non-stressed environment in RIL population (ICC 4567 × ICC 15614), Table S4: Consistent QTLs found across heat-stress environments (2013 and 2014) in RIL population (ICC 4567 × ICC 15614), Table S5: Gene ontology classification for CaLG05,  Table S6: Gene ontology classification for CaLG06, Table S7: Gene ontology categorization of 236 genes identified on the genomic region flanked by markers Ca5_44667768-Ca5_46955940 on CaLG05, Table S8: Gene ontology categorization of 550 genes identified on the genomic region flanked by markers Ca6_7846335-Ca6_14353624 on CaLG06, Table S9: (a) List of putative candidate genes found to be associated with heat stress on CaLG05 in chickpea, (b) List of putative candidate genes found to be associated with heat stress on CaLG06 in chickpea.
Author Contributions: P.M.G. conceived the idea and coordinated this project. P.M.G. and S.S. were involved in developing the mapping population. P.M.G., S.S., S.K.C., S.B.S. and G.R.L. provided guidance to P.J.P. in conducting field experiments and phenotyping. A.R., R.R.D. and S.S. helped P.J.P. in statistical data analysis. P.J.P. and M.T. were involved in genotyping of the mapping population, construction of linkage maps and QTL analysis. P.J.P., A.W.K. and M.T. were involved in bioinformatics work. P.J.P., P.M.G., MT, S.B.S. and S.S. contributed to writing of the manuscript, and R.K.V. and G.R.L. provided their inputs. All the authors reviewed and approved the final manuscript.