Major QTLs and Potential Candidate Genes for Heat Stress Tolerance Identified in Chickpea (Cicer arietinum L.)

In the context of climate change, heat stress during the reproductive stages of chickpea (Cicer arietinum L.) leads to significant yield losses. In order to identify the genomic regions responsible for heat stress tolerance, a recombinant inbred line population derived from DCP 92-3 (heat sensitive) and ICCV 92944 (heat tolerant) was genotyped using the genotyping-by-sequencing approach and evaluated for two consecutive years (2017 and 2018) under normal and late sown or heat stress environments. A high-density genetic map comprising 788 single-nucleotide polymorphism markers spanning 1,125 cM was constructed. Using composite interval mapping, a total of 77 QTLs (37 major and 40 minor) were identified for 12 of 13 traits. A genomic region on CaLG07 harbors quantitative trait loci (QTLs) explaining >30% phenotypic variation for days to pod initiation, 100 seed weight, and for nitrogen balance index explaining >10% PVE. In addition, we also reported for the first time major QTLs for proxy traits (physiological traits such as chlorophyll content, nitrogen balance index, normalized difference vegetative index, and cell membrane stability). Furthermore, 32 candidate genes in the QTL regions that encode the heat shock protein genes, heat shock transcription factors, are involved in flowering time regulation as well as pollen-specific genes. The major QTLs reported in this study, after validation, may be useful in molecular breeding for developing heat-tolerant superior lines or varieties.


INTRODUCTION
Given the global climate changes, heat stress is becoming a major challenge to crop production and food safety. As per Intergovernmental Panel on Climate Change, the current rate of global warming is 0.2 • C per decade and is predicted to reach 1.5 • C between 2,030 and 2,052 (https://www.bbc. com/news/newsbeat-48947573). Such an increase in temperatures leads to heat stress and costs the global economy US$2.4 trillion a year (https://news.un.org/en/story/2019/07/1041652). More than 15% of the global land area becomes exposed to high levels of heat stress with an additional 0.5 • C increase to the 2 • C (Sun et al., 2019). Heat stress, besides affecting producers directly, also reduces labor productivity (Kjellstrom, 2016), further compounding the effects of increasing temperature on crop yields. In recent years, shifts toward more sustainable and healthy diets, which are typically characterized by high consumption of vegetables and legumes, have been evidenced (Scheelbeek et al., 2018).
Chickpea (Cicer arietinum L.) is an important cool season grain legume crop cultivated in the arid and semi-arid regions across the globe. It is an excellent source of proteins, essential amino acids, vitamins, and minerals . Major chickpea producing countries are India, Australia, Pakistan, Turkey, Russia, Myanmar, Iran, Mexico, Canada, and USA. In India, Madhya Pradesh, Maharashtra, Rajasthan, Karnataka, Uttar Pradesh, and Andhra Pradesh are the major chickpea growing states. Although India is the largest producer of chickpea, in order to attain self-sufficiency by 2050, the chickpea production in the country needs to reach 16-17.5 Mt from an area of about 10.5 Mha with an average productivity of 15-17 q/ha (Dixit et al., 2019). Drought and heat, the two most important environmental factors, can cause more than 70% yield loss in chickpea (Varshney et al., 2019). Traditionally, chickpea requires prolonged winter for better growth and cultivation in the northern states of India. However, in the northern states, the pulse area especially chickpea cultivation was reduced due to the green revolution. Southern and Central India, where significant chickpea area increased, are exposed to drought and heat stresses. The rise in ambient air temperature (≥35 • C) that coincides with the reproductive processes leads to various anomalies in reproductive events, especially during fertilization, pod formation, and pod filling in chickpea (Devasirvatham et al., 2013;Kaushal et al., 2013;Gaur et al., 2014).
The genetic mechanism of heat stress in different crop plants has been reviewed extensively (see Janni et al., 2020). In general, the impact of heat stress depends on the intensity, duration of exposure, and degree of the elevated temperature. In the case of legumes like chickpea, heat stress has deleterious effects on the morphology, physiology, and reproductive growth (Sita et al., 2017). The effects of heat stress on the development of various male and female tissues in different legume species have been reviewed recently (Liu et al., 2019). In the case of legume crops, heat shock proteins (HSP), HSP gene families, and various metabolites were reported to control heat stress response (see Janni et al., 2020). Heat stress adversely affects pollen viability, fertilization, and seed development, which leads to a reduced harvest index consequently, and these events greatly impact chickpea yield. In the cool season, legumes such as chickpea, lentil, faba bean, and field peas, the temperature above 30 • C lead to yield losses (Jiang et al., 2015;Bishop et al., 2016;Bhandari et al., 2017). As heat stress is a complex trait governed by many genes/QTLs, breeding for heat stress tolerance in chickpea remains challenging (Krishnamurthy et al., 2011;Devasirvatham et al., 2013). Therefore, the effects of heat stress on chickpea growth, development, and yield are important to understand by observing agronomic traits to develop high-temperaturetolerant cultivars.
Genomic revolution, during the last two decades, simplified understanding of the complex responses to biotic and abiotic stress in several crop plants (Roorkiwal et al., 2020;Thudi et al., 2020). Chickpea research community has access to genome sequence (Varshney et al., 2013), genome-wide variations among diverse germplasm lines at the sequence level (Thudi et al., 2016a,b;Varshney et al., 2019) for trait dissection, and the development of climate-resilient chickpea varieties (Mannur et al., 2019;Bharadwaj et al., 2020). The genotyping-bysequencing (GBS) approach has been extensively used for singlenucleotide polymorphism (SNP) discovery and mapping traits in several crops for genetic research and breeding applications (Chung et al., 2017), including chickpea (Jaganathan et al., 2015Thudi et al., 2020). Besides proteomic and metabolomic approaches to understanding the molecular mechanism of heat tolerance (Parankusam et al., 2017;Salvi et al., 2018), efforts were made to map QTLs and markers associated with heat tolerance in chickpea Paul et al., 2018a;Varshney et al., 2019;Roorkiwal et al., 2020).
In this study, we reported the construction of a high-density genetic map using SNPs derived from the GBS approach and major QTLs for phenological, physiological, yield, and yieldrelated traits based on phenotyping of recombinant inbred line (RIL) population (DCP 92-3 × ICCV 92944) under two environments (normal and late sown) for 2 years (2017-2018 and 2018-2019). In addition, we also reported the potential candidate genes implicated for heat tolerance in the QTL regions.

Plant Material
A biparental mapping population, comprising 184 F 7 RIL lines, derived from the cross DCP 92-3 × ICCV 92944 segregates for heat tolerance was used for identifying genomic regions and candidate genes for heat tolerance. DCP 92-3 is a logging and Fusarium wilt-resistant variety released by the Indian Council of Agricultural Research (ICAR)-Indian Institute of Pulses Research (IIPR), Kanpur, Uttar Pradesh, India for cultivation in Punjab, Haryana, Delhi, Northern Rajasthan, and Western Uttar Pradesh. Pollen viability at a critical temperature of 35 • C differentiates the heat-sensitive and heat-tolerant genotypes. Based on physiological, biochemical, yield, and yield-related trait studies conducted earlier Kumar et al., 2012;Kaushal et al., 2013;Bhandari et al., 2020), the chickpea genotype ICCV 92944 was reported as a heat-tolerant genotype and was released as BARI Chola-10 in Bangladesh, as Yezin 6 in Myanmar, and as JG 14 in India and is performing well under late sown conditions.

Phenotyping of Recombinant Inbred Line Population
In the case of chickpea, the optimal temperature for its growth ranges between 10 and 30 • C. Chickpeas are sensitive to heat stress particularly at the reproductive phase (flowering and seed development). A few days of exposure to high temperatures (35 • C or above) during the reproductive phase can cause heavy yield losses through flower and pod abortion. Late sowing, a simple and effective field screening technique for reproductivestage heat tolerance in chickpea developed at the International Crops Research Institute for the Semi-Arid Tropics (ICRISAT), Patancheru, India, was adopted for phenotyping the RILs for heat stress tolerance. The F 7 RILs (184 individuals) and parents DCP 92-3 × ICCV 92944 were evaluated for two consecutive years 2017-2018 and 2018-2019 under normal sown environment (NS; second week of November) and late sown or heat stress environment (HTS; third week of December) at ICAR-IIPR, Kanpur, Uttar Pradesh, India (26 • 26 ′ 59.7228 ′′ N and 80 • 19 ′ 54.7356 ′′ E). The experiments were conducted under field conditions in a plot admeasuring 3 × 0.6 m, and the distance between plants is 10 cm. The RIL population was evaluated in augmented block design along with the parents DCP 92-3 and ICCV 92944 and two elite chickpea genotypes JG 11 and ICC 4958. All the individuals of the population were apportioned into a total of 10 blocks along with the four checks replicated in each block. The maximum, minimum, and mean temperatures were recorded weekly during the entire cropping season for both years (Supplementary Figure 1). The mapping population was phenotyped for physiological traits like normalized difference vegetation index (NDVI;using GreenSeeker, Optical Sensor Unit, 2002 114 NTech Industries, Ukiah, USA), nitrogen balance index (NBI, using DUALEX R optical leafclip meter), NBI R combines chlorophyll and flavonols (related to nitrogen/carbon allocation) measured by using DUALEX R optical leafclip meter, chlorophyll content (CHL, using DUALEX R optical leafclip meter ng/mm 2 ) and cell membrane stability (CMS, %), yield, and yield-related traits [(total filled pods per plant (FP), biological yield per plant (BYPP, g), seed yield per plant (SYPP, g), harvest index (HI, %), and 100 seed weight (100SDW, g)]. To avoid the biasness, the mean of 10 individual plants was sampled for seed yield/plant taken from each planted genotype instead of seed yield/m 2 per plot. Furthermore, the mean of 10 plants randomly chosen from each line was used for recording the abovementioned traits for all the individuals of mapping population under NS and HTS for both years. Two irrigation and same agronomic package of practice were followed for both NS and HTS sown genotypes for both years. NDVI was measured as per the following formula NDVI = NIR-RED\NIR + RED (Myneni et al., 1995), and CMS was measured as per the formula used by Blum and Ebercon (1981). CMS = 100-membrane injury index (MII), where MII is calculated as a ratio of C1 and C2, with C1 and C2 denoting the electrolytes measured at 40 and 80 • C, respectively.

Analysis of Variance, Best Linear Unbiased Prediction (BLUP), and Heritability
The ANOVA for the RIL population was performed using GenStat (17th Edition), for individual environments using the 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 effects. To pool the data across environments and to make the error variances homogeneous, the individual variances were estimated and modeled for the error distribution using the residual maximum likelihood (ReML) procedure. The Z-value and F-value were calculated for random effects and fixed effects, respectively. Broad-sense heritability was calculated as H 2 = Vg/(Vg + Ve/nr), as suggested by Falconer et al. (1996), and pooled broad-sense heritability was estimated as H 2 = Vg/{(Vg) + (Vge/ne +Ve/(ne × nr))}, as suggested by Hill et al. (2012), where H 2 is the broad-sense heritability, Vg is the genotypic variance, Vge is the G × E interaction variance, Ve is the residual variance, ne is the number of environments, and nr is the number of replications.
DNA Extraction, Genotyping, and Single-Nucleotide Polymorphism Calling DNA from 184 RILs, along with the parents, was isolated from 2-week-old seedlings following the high-throughput mini-DNA extraction method (Cuc et al., 2008). The quality of DNA was checked by using 0.8% agarose gel electrophoresis, and the quantity was assessed by Qubit R 2.0 Fluorometer (Thermo Fisher Scientific Inc., USA). The GBS approach was used for SNP calling between the parents and genotyping the RILs as described by Elshire et al. (2011). GBS libraries from the parental lines and RILs were prepared using ApeKI endonuclease (recognition site: G/CWCG), followed by ligation with uniquely barcoded adapters using T4 DNA ligase enzyme. Such digested ligated products from each sample were mixed in equal proportion to construct the GBS libraries, which were then amplified, purified to remove excess adapters, and used for sequencing on HiSeq 2500 platform (Illumina Inc., San Diego, CA, USA). Sequence reads from raw FASTQ files were used for SNP identification and genotyping using the reference-based GBSv2 analysis pipeline implemented in TASSEL v5.0 (Bradbury et al., 2007). In brief, all reads that begin with one of the matched barcodes immediately followed by the expected four base remnants of the enzyme cut site are sorted, de-multiplexed, and trimmed to first 64 bases starting from the enzyme cut site. Reads containing N within the first 64 bases after the barcode are rejected. The remaining good quality reads (called as tags) were aligned against the draft genome sequence of chickpea using Bowtie2 software. The alignment file was then processed by using the GBSv2 analysis pipeline for SNP calling and genotyping.

Linkage Map Construction and Identification of QTLs
In order to construct the genetic map, all markers were grouped into eight linkage groups with the logarithm of odds (LOD) threshold of 5.0. Marker order within a linkage group was assigned using the regression mapping algorithm with a maximum recombination frequency of 0.4 at a LOD of three and a jump threshold of five. The ripple command was finetuned by adding each marker locus to confirm the final marker order. The Kosambi mapping function was used to calculate the map distance in centimorgan (cM). The segregation distortion and chi-square (χ 2 ) values were detected using JoinMap V4.0, and markers with heterozygosity and significant segregation distortion were excluded (p < 0.001) from the analysis. The linkage map was constructed using ICIMapping 3.2 software (Meng et al., 2015). The QTL analysis was conducted for NS 2017, NS 2018, NS pooled data, HTS 2017, HTS 2018, and HTS pooled data together with the genotyping data and genetic map information using software windows QTL Cartographer version 2.5 (Wang et al., 2012). The composite interval mapping (CIM) analysis was conducted by scanning the intervals of 1.0 cM between markers and putative QTLs with a window size of 10.0 cM and by using the parameters of model six and 1,000 times of permutation with the 0.05 significance level along with the function of "Locate QTLs" option to locate QTLs.

Identification of Candidate Genes Within QTL Confidence Intervals
Based on the physical position of the SNPs/markers flanking the QTL regions, the candidate genes present within the determined QTL intervals were retrieved from the draft genome sequence (CaGAv1.0) of chickpea (Varshney et al., 2013). The identified genes in QTL intervals were searched against NCBI-nr protein database using BLAST program. The gene ontology (GO) terms associated with the genes were searched for GO terms, using BLAST2GO software (Conesa et al., 2005).

Phenotypic Performance and Genetic Variability of the Parents and Mapping Population
A considerable amount of genetic variation for various phenological, yield, and yield-related traits was observed in both the parents and the derived RILs under NS and HTS environments for both years. The descriptive statistics are shown in Table 1. Transgressive segregates in both directions were observed for days to flower initiation (DFI), FP, and SYPP traits in the RIL population (Figure 1)

Relationships Among Different Traits
To investigate the relationship among different traits, we calculated the pairwise correlations among different traits within each environment (NS and HTS). During 2017-2018, under HTS environment, a positive and high significant correlation was observed between DFI with that of days to pod initiation (DPI) (p < 0.01) and days to maturity (DM) (p < 0.01) ( Supplementary Table 1). Similarly, during 2018-2019, under HTS environment, a positive and high significant correlation was also observed between DFI with that of DPI (p < 0.01) and DM (p < 0.05) (Supplementary Table 2). Furthermore, during 2017-2018, under HTS environment, NBI and CHL were found to possess a positive and high significant correlation (p < 0.01). However, no significant correlation was observed during 2018-2019 under HTS environment. A number of filled pods (NFP) and SYPP had a significant and positive correlation under heat stress environments during both years. Furthermore, NPF has a significant positive correlation with BYPP in 2017-2018 and with HI in 2018-2019 (Supplementary Tables 1, 2). Nevertheless, 100SDW possesses a positive and high significant correlation with HI and SYPP during both years under HTS environment (Supplementary Tables 1, 2). Similar positive and high significant correlations were also observed under NS environments in both years as well as pooled data of NS environments for the abovementioned traits (Supplementary Tables 3-6).

Single-Nucleotide Polymorphisms-Based Genetic Map
A total of 49.89 Gb (49 million reads) clean GBS reads were generated using HiSeq2500 on the RIL population derived from DCP 92-3 × ICCV 92944. The number of reads generated per individual ranged from 0.86 to 5.3 million. A total of 3,425,458 genome-wide SNPs were identified on aligning the data to CDC Frontier reference genome (Varshney et al., 2013) using TASSEL-GBS pipeline. After excluding ambiguous SNP calls, SNPs that are monomorphic among the parental genotypes, and SNPs with segregation distortion, a total of 7,947 polymorphic SNPs were used for the linkage map analysis using ICIM. As a result, a genetic map comprising 788 SNPs distributed on eight linkage groups (CaLG01-CaLG08) spanning 1,125 cM was constructed (

QTLs for Heat Stress Tolerance Traits
By using CIM, a total of 77 QTLs (37 major QTLs and 40 minor QTLs) were identified for 12 of 13 traits phenotyped for two seasons (2017-2018 and 2018-2019) and two environments (NS and HTS). Of 77 QTLs, 37 QTLs were major explaining ≥10% phenotypic variation (PVE), and 40 QTLs were minor explaining <10% PVE ( Table 2). A positive value of the additive variance of a given QTL indicates that the female parent (DCP 92-3) has a positive effect on the trait; while a negative value indicates that the male parent (ICCV 92944) having a positive effect on the trait.

QTLs for Physiological Traits
A total of 36 (17 major and 19 minor) QTLs were identified for physiological traits with PVE ranging from 6.69% to 34.02%.

Normalized Difference Vegetation Index (NDVI)
A total of 16 QTLs (seven major with PVE 10.27-34.02% and nine minor with PVE 6.69-9.85%) were identified for NDVI, out of which six were identified based on HTS environment and 10 were identified based on NS environment. Interestingly, for this trait, QTLs were identified on all linkage groups except CaLG07. Furthermore, the majority of QTLs (25%) were present on CaLG02 that explained 8.84-26.31% PVE. However, the QTL on CaLG01 that explained the highest PVE (34.02%) among all QTLs identified for this trait was based on pooled HTS environment ( Figure 3A).

Nitrogen Balance Index (NBI)
A total of 10 QTLs (five major with PVE 10.26-13.93% and five minor with PVE 7.39-9.95%) were identified for NBI. Among these 10 QTLs, five were on CaLG08, three on CaLG07, and two on CaLG06. Of five QTLs identified on CaLG08, four QTLs were flanked by SCA8_6301805 and SCA8_11012719 markers and one QTL was flanked by SCA8_6301805-SCA8_11012719 markers. Furthermore, among these five QTLs, two each were identified based on the pooled data from HTS environments and HTS of 2018-2019 and one based on HTS 2017-2018 ( Figure 3B).

Chlorophyll Content (CHL)
Of the seven QTLs identified for CHL, four QTLs were on CaLG04, two were on CaLG05 and one on CaLG02. Under HTS 2018, one minor QTL (PVE 8.14%) was identified on CaLG02.

Cell Membrane Stability (CMS)
Of three QTLs identified for CMS, one QTL each was on CaLG03, CaLG04, and CaLG06. Among these QTLs, two were identified based on pooled data from the NS environment and one based on pooled data from the HTS environment.

QTLs for Yield and Yield-Related Traits
Eighteen major and 16 minor QTLs were identified for yield and yield-related traits with PVE ranging from 5.88 to 43.49% ( Table 2).

Days to Pod Initiation
Under HTS environments (2017)(2018)(2018)(2019), and pooled data), a total of 6 QTLs (three major with PVE 10.33-43.49% and three minor with PVE 5.88-8.45%) were identified for DPI. Out of these, one QTL was on CaLG01, three on CaLG06, and two on CaLG07. Furthermore, all QTLs identified on CaLG06 were flanked by SCA6_39028647-SCA6_43908965 markers, while QTLs on CaLG07 were flanked by SCA7_9555338-SCA7_47907019 markers. However, the QTL explaining the highest proportion of PVE was present on CaLG07 (Figure 4).

Days to Pod Filling and Number of Filled Pods
Of three QTLs (two major with PVE 11.96-11.97% and one minor with PVE 9.38%), two were on CaLG04 and one was on CaLG08. The major QTL on CaLG08 was identified under HTS 2018, while the remaining two were based on the pooled data of the NS environment. However, all three QTLs were flanked by different markers ( Table 2). One minor QTL (PVE 6.6%) for NFP was identified on CaLG06 based on the pooled data of the NS environment.

Seed Yield per Plant (g)
A total of four QTLs (three major with PVE 11.88-18% and one minor with PVE 8.66%) were identified for SYPP, of which two major QTLs were under the HTS 2018 environments and the remaining two were based on pooled data of the NS environment.

Biological Yield per Plant (g)
A total of eight QTLs (three major with PVE 10.7-11.16% and five minor with PVE 6.92-9.43%) were identified for BYPP. Of eight QTLs, five were identified in the HTS environments of 2017-2018 and 2018-2019 and based on pooled data, and three were identified in the NS environments and based on pooled data. Among these QTLs, seven QTLs were present on CaLG06 and one on CaLG02. Furthermore, all eight QTLs were flanked by different markers (Table 2; Figure 4A).

Harvest Index (%)
Of eight QTLs, one minor QTL each was identified for HI under the HTS environment of 2018 and pooled data of both years, while the remaining six QTLs were on the NS environments of 2017-18 and 2018-19 and pooled data of both years. Among six QTLs under the NS environments, three QTLs were major with PVE 12.38-39.31% and three were minor with PVE 8.08-9.24%. Furthermore, among eight QTLs, three QTLs were located on CaLG06, three on CaLG07, and one each on CaLG05 and CaLG08.

Seed Weight (g)
A total of three major QTLs were identified for 100SDW under HTS 2017 (one QTL) and based on pooled data (three QTLs) under HTS environments for both years. Among four QTLs, two were located on CaLG07 (Figure 4B), while one each was located on CaLG01 and CaLG04. Furthermore, the PVE for these four QTLs ranged from 31.3 to 37.23%.  Table 9). Based on functional categorization, many genes were found to be associated with biological processes in these genomic regions. Using GO classification, we further identified a total of 32 candidate genes (7 on CaLG01, 3 on CaLG02, 14 on CaLG04, and 8 on CaLG07) known to function, directly or indirectly, as heat-stress response genes in chickpea (Table 3). Among seven genes on CaLG01, six genes were encoding heat shock proteins, while one gene was encoding pollen-specific leucine-rich repeat extensin-like protein 1. While in the case of CaLG02, of three selected candidate genes, Ca_16007 encoded pollen-specific leucine-rich repeat extensin-like protein 1, Ca_24649 encoded a truncated transcription factor CAULIFLOWER A-like, and Ca_22033 encoded heat shock protein-binding protein. Among 14 selected candidate genes on CaLG04, six genes were pollen-specific, four were related to heat shock protein, three were DnaJ heat shock amino-terminal domain protein, and one was related to protein PHOTOPERIOD-INDEPENDENT EARLY FLOWERING 1 isoform X1 ( Table 3). The eight genes on CaLG04 encode heat shock protein/heat shock factor protein HSF24-like (Ca_18924, Ca_16239, and Ca_09277), pollen-specific leucine-rich repeat extensin-like protein 1/pollen receptor-like kinase 3 (Ca_16434 and Ca_16155), protein EARLY FLOWERING 3/flowering time control protein FY (Ca_10118 and Ca_17996), and calmodulin-binding heat-shock protein (Ca_13761) ( Table 3).

DISCUSSION
In the context of climate change, every degree increase in aerial temperature has a severe impact on crop production, especially on chickpea that is predominantly cultivated on residual soil moisture in marginal environments . Therefore, understanding the nature, impact, and molecular mechanisms of heat stress tolerance will help in designing strategies to overcome production losses. In chickpea, previously, very few studies were focused on understanding the nature, impact, and existing diversity in germplasm lines as well as identifying the genomic regions responsible to some extent. In this study, we reported major QTLs and novel genes in these genomic regions, which are associated with responsible for heat stress tolerance. Late sowing, a simple and effective field screening technique for reproductive-stage heat tolerance in chickpea developed at the ICRISAT, Patancheru, India, was adopted for phenotyping the RILs for heat stress tolerance. The late sowing approach was adopted earlier in understanding the genetic variability for heat stress among genotypes as well as in identifying the genomic regions responsible for heat stress tolerance; for instance, in cool season crops, namely, chickpea (Paul et al., 2018b), wheat (Sareen et al., 2020), brassica (Branham et al., 2017), and rice (Prasanth et al., 2016). As the selection was based on yield per se results in a slower response because of genotype × environmental interactions, we also phenotyped the mapping population for physiological traits like CMS, NDVI, NBI, and CHL, which could be used as an indirect selection criterion to improve heat tolerance in chickpea as this was used in other crop plants. A sufficient amount of genetic variability for various phenological, physiological, and yield-related traits was noted in both the parents, and the RIL population was studied under NS and HTS environments for both years. The similar results were also reported earlier in chickpea based on evaluating the germplasm lines as well as one RIL population under HTS environment (Krishnamurthy et al., 2011;Devasirvatham et al., 2013;Paul et al., 2018b). High heritability for physiological traits like CMS, NDVI, NBI, and CHL contents under HTS environments indicates that the selection for heat tolerance relying on these traits could be effective. Earlier, heat stress was reported to reduce the total CHL and showed moderate to high heritability for  NDVI, CMS, and CHL content under stress condition in wheat (Bhusal et al., 2018;Condorelli et al., 2018;Pradhan et al., 2020), maize (Naveed et al., 2016), carrot (Nijabat et al., 2020), and pea (Tafesse et al., 2020). Low heritability for HI trait under HTS environments observed in this study indicates that the election for this trait will not enhance yield under stress. As yield traits remain the primary objective for improving the heat tolerance in all crop plants including chickpea, a positive and significant correlation among the yield and yield-related traits especially, FP, SYPP, and BYPP could serve as an important parameter for developing heat-tolerant chickpea genotype. In this study, using a RIL population with a dense genetic map and phenotyping under NS and HTS allowed us to precisely identify major QTLs for heat stress in chickpea. Our genetic map has approximately 3-fold more markers compared to the previous study reporting the QTLs for heat stress tolerance (Paul et al., 2018a). In addition to reporting QTLs for phenological, yield, and yield-related traits under heat stress environments, to the best of our knowledge, this study is the first comprehensive study that reports QTLs for physiological traits like CHL, NBI, NDVI, CMS, and NPP in chickpea. A better understanding of phenology in response to heat stress will enable designing the breeding strategies. Minor QTLs for DFI were identified on CaLG06 and CaLG08, while a major QTL was identified for DM on CaLG01. Physiological traits like CMS, NDVI, and NBI, to date, have been used as proxy for grain yield under stress mostly in cereals (ElBasyoni et al., 2017;Bhusal et al., 2018;Condorelli et al., 2018;Getahun et al., 2020;Khanna-Chopra and Semwal, 2020). In this study, we reported major QTLs for these traits in chickpea, which may be used, after validation, for marker-assisted breeding for heat stress tolerance in chickpea.
In the case of chickpea, four flowering time genes (efl-1 from ICCV 96029, efl-3 from BGD 132, and efl-4 from ICC 16641) and their allelic relationships were reported , and major QTLs corresponding to these genes were mapped on CaLG04, CaLG08, and CaLG06, respectively (Mallikarjuna et al., 2017). Furthermore, marker trait associations for flowering time were reported earlier Upadhyaya et al., 2015;Varshney et al., 2019). However, in this study, we reported QTLs for DFI on CaLG06 and CaLG08 under heat stress environments as well as based on the pooled data of 2 years. Furthermore, interestingly, QTLs for DPI and DFI co-localized or mapped in the same genomic region under HTS environments of both years as well as based on pooled data. These observations indicate that introgression of one of the traits simultaneously improves both the traits, which are key for achieving resilience to heat stress.
Earlier two major genomic regions harboring QTLs for heat tolerance-related traits were mapped on CaLG05 and CaLG06; however, none on the QTLs explained >20% PVE (Paul et al., 2018a). Nevertheless, HI in this study explained >30% PVE. Similarly, CaLG04 also harbored QTLs for five traits (CHL, CMS, NDVI, DPF, and 100SDW), among these QTLs for traits CHL and 100SDW explained >30% PVE. Except HI (PVE 39.13%, NS 2018), none of the QTLs mapped on CaLG06 had PVE >15%. Similarly, QTLs for traits like CHL, NDVI, and HI were mapped on CaLG05 which explained 6.78-18.69% PVE. Earlier, a genomic region refereed as "QTL-hotspot" was reported to harbor several QTLs for different drought tolerance-related traits including 100SDW on CaLG04 . A genomic region on CaLG07 harbors QTLs explaining >30% PVE for DPI and 100SDW as well as QTL for NBI explaining >10% PVE. For 100SWD, a total of four major QTLs were identified on CaLG01, CaLG04, and CaLG07 under HTS, and no QTLs were detected under NS. For SYPP trait, two major QTLs were identified on CaLG02 and CaLG06 under HTS. Only one QTL was identified on CaLG06 under NS. However, yield-related QTLs were not consistently recorded under all the conditions suggesting their environmental specific expression. Likewise, QTLs contributing to pods per plant, seed yield per plant, and seed weight were reported on CaLG01 and CaLG06 Srivastava et al., 2016). In the case of cowpea, four QTLs were identified for pod set number per peduncle under HS and markers, which were utilized in breeding applications (Lucas et al., 2013;Pottorff et al., 2014). Similarly, in lentils, a major QTL controlling the seedling survival and pod setting traits under heat stress was noticed . In addition, QTLs for SYPP and BYPP (full names) were mapped in the same genomic region under NS environments of both years as well as based on pooled data. For yield and yield-related traits like DPI, DPF, and SYPP under HTS major alleles were contributed by ICCV 92944. For the 100SDW trait under HTS major alleles were contributed by DCP 92-3. On the other hand, almost all of the traits DPI, BYPP, and SYPP were contributed by DCP 92-3 under NS. For the trait HI under NS major alleles was contributed by both parents. In addition to these major QTLs, several QTLs were identified that were environmentally specific under NS and HTS, which only appeared in this study in the first year (NS or HTS). A total of nine major QTLs were located on CaLG06, which highlight the importance of this region in the heat tolerance mechanism in chickpea. Some QTLs were largely affected by environmental factors and that could be detected in only one season, and for these QTLs, further verification should be required.
HSP genes play a pivotal role in heat stress tolerance. In this study, 32 genes were identified in the QTL regions of CaLG01, CaLG02, CaLG04, and CaLG07. Similarly, in the case of soybean, 38 Hsfs were identified that were located on 15 soybean chromosomes (Li et al., 2014). HSP gene families were reported to be involved in drought and heat stress responses in soybean seedlings Das et al., 2016;Liu et al., 2019). HSP90 gene families in five legumes and expression profiles in chickpea were reported earlier (Agarwal et al., 2016). Furthermore, based on genome-wide association studies, especially eight flowering time-regulating genes (efl1, FLD, GI, Myb, SFH3, bZIP, bHLH, and SBP) were reported (Upadhyaya et al., 2015). The genes reported in this study can be further explored for haplotypes based on the germplasm sequence information available in the public domain that has the potential for genetic improvement of the trait (Varshney et al., 2020). In addition, pollen-specific leucine-rich repeat extensinlike protein 1 genes identified in the QTL regions were reported to synergistically maintain pollen tube cell wall integrity; thus, they play critical roles in pollen germination and pollen tube growth (Wang et al., 2018). Recently, cloning of SHY in tomato, a pollen-specific gene that encodes a leucine-rich repeat (LRR) protein, demonstrated its role in a signal transduction pathway mediating pollen tube growth (Guyon et al., 2004).

CONCLUSIONS
In this study, we identified a total of 37 major QTLs across the genome for 12 traits. DFI, DPI, and DM are the key traits for escaping the heat stress in chickpea especially reproductive heat stress that hampers chickpea production. In this study, we reported major QTLs explaining >30% PVE for these key traits that contribute to yield under heat stress. In addition, we also reported for the first time major QTLs for proxy traits (physiological traits like CHL, NBI, NDVI, and CMS). Furthermore, 32 candidate genes in the QTL regions that encode the HSP, heat shock transcription factors, genes are involved in flowering time regulation as well as pollen-specific genes. The major QTLs reported in this study may be useful in molecular breeding for developing heat-tolerant superior lines or varieties.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: NCBI repository, BioProject ID PRJNA695065 (http://www.ncbi.nlm.nih.gov/ bioproject/695065).

AUTHOR CONTRIBUTIONS
UJ and MT conceived the idea. UJ, HN, and RJ conducted the experiments. UJ performed the statistical analysis. MT, UJ, RP, and RV prepared the manuscript. VV and PB performed the GBS analysis. NS contributed to the development and phenotyping of mapping population. AC and RV contributed to consumables and the generation of genotyping data. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We are grateful to the Bill and Melinda Gates Foundation and Science and Engineering Research Board (SERB), Government of India for partial support to this study. This work has been carried out as part of the CGIAR Research Program on Grain Legumes and Dryland Cereals. ICRISAT is a member of the CGIAR System Organization.

SUPPLEMENTARY MATERIAL
The Temperature >30 • C during reproductive stages especially at podding stage and pod filling were recorded, which is more than the critical temperatures that hamper the production of chickpea.
Supplementary Figure 2 | High density genetic map comprising 788 single-nucleotide polymorphism markers mapped on eight chromosomes of chickpea. The map distance is indicated as a common scale for all eight linkage groups (CaLG01-CaLG08) on the left side figure and marker names on the right side of each linkage group.