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

Jha, Uday Chand; Nayyar, Harsh; Palakurthi, Ramesh; Jha, Rintu; Valluri, Vinod; Bajaj, Prasad; Chitikineni, Annapurna; Singh, Narendra P.; Varshney, Rajeev K.; Thudi, Mahendar

doi:10.3389/fpls.2021.655103

ORIGINAL RESEARCH article

Front. Plant Sci., 26 July 2021

Sec. Plant Breeding

Volume 12 - 2021 | https://doi.org/10.3389/fpls.2021.655103

This article is part of the Research Topic Accelerating Genetic Gains in Pulses View all 15 articles

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

$\nUday Chand Jha$ Uday Chand Jha¹

Rintu Jha⁴

Annapurna Chitikineni³

Narendra P. Singh¹

Rajeev K. Varshney^3,5

Mahendar Thudi^3,6^*

¹Indian Council of Agricultural Research (ICAR)-Indian Institute of Pulses Research (IIPR), Kanpur, India
²Department of Botany, Panjab University, Chandigarh, India
³Center of Excellence in Genomics and Systems Biology (CEGSB), International Crops Research Institute for the Semi-Arid Tropics (ICRISAT), Hyderabad, India
⁴Institute of Crop Science, Chinese Academy of Agricultural Science (CAAS), Beijing, China
⁵State Agricultural Biotechnology Centre, Centre for Crop and Food Innovation, Food Futures Institute, Murdoch University, Murdoch, WA, Australia
⁶University of Southern Queensland, Toowoomba, QLD, Australia

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 (Jukanti et al., 2012). 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-temperature-tolerant 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-by-sequencing (GBS) approach has been extensively used for single-nucleotide polymorphism (SNP) discovery and mapping traits in several crops for genetic research and breeding applications (Chung et al., 2017), including chickpea (Jaganathan et al., 2015; Thudi 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 (Thudi et al., 2014; 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 yield-related 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.

Materials and Methods

Plant Material

A biparental mapping population, comprising 184 F₇ 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 (Gaur et al., 2012; 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 reproductive-stage 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₇ 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^® optical leafclip meter), NBI^® combines chlorophyll and flavonols (related to nitrogen/carbon allocation) measured by using DUALEX^® optical leafclip meter, chlorophyll content (CHL, using DUALEX^® optical leafclip meter ng/mm²) 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² 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.

Statistical Analyses

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² = Vg/(Vg + Ve/nr), as suggested by Falconer et al. (1996), and pooled broad-sense heritability was estimated as H² = Vg/{(Vg) + (Vge/ne +Ve/(ne × nr))}, as suggested by Hill et al. (2012), where H² 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^® 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 fine-tuned 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 (χ²) 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).

Results

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). The Combined ANOVA indicated the presence of significant genetic variability in the evaluated RILs under both NS and HTS. High to moderate heritability (98.2–61.3%) under NS for both years and 73.3–98.4% heritability under HTS for both years were recorded. Only low heritability of 38.2 and 47.9% for HI was observed under HTS during 2017–2018 and 2018–2019. However, high heritability (77.6–84.7%) was noted under NS conditions.

TABLE 1

Table 1. Phenotypic performance of heat sensitive (DCP 92-3), tolerant (ICCV 92944), RILs, and heritability of traits evaluated under normal and heat stress environments.

FIGURE 1

Figure 1. Frequency distribution of days to flower initiation (DFI, d), total filled pods per plant (FP), and seed yield per plant (SYPP, g) in RIL population derived from DCP 92-3 × ICC 92944 and evaluated, under heat stress environment 2017–18 (A) and under heat stress environment 2018–19 (B).

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 (Figure 2; Supplementary Tables 7, 8; Supplementary Figure 2). The CaLG06 had the highest proportion of the mapped SNPs (23.4%; 185 SNPs), whereas CaLG08 had the lowest proportion of the mapped SNPs (7.6%; 60 SNPs) and the largest linkage group CaLG01 spanned 191 cM, whereas the smallest linkage group CaLG08 spanned 68 cM.

FIGURE 2

Figure 2. Genetic map and QTLs for heat tolerance-related traits. A genetic map comprising 788 SNPs distributed in eight linkage groups (CaLG01–CaLG08) spanning 1,125 cM was constructed using RIL population derived from DCP 92-3 × ICCV 92944. The blue and red-colored bars indicate the QTLs identified under normal and late sown conditions, respectively. The major QTLs are indicated in bold font and the minor QTLs in italics.

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.

TABLE 2

Table 2. Summary of QTLs identified for phenological, physiological, yield, and yield-related traits using RIL population derived from DCP 92-3 × ICCV 92944.

QTLs for Phenological Traits

Under the HTS environment, in the case of DFI, two QTLs each were identified during 2018–2019 and pooled data of 2017–2018 and 2018–2019 on CaLG06 and CaLG08. The PVE ranged from 7.48 to 8.96%. While in the case of DM, all three QTLs (PVE 8.96–18.13%) identified were in 2017–2018 and under HTS environment on CaLG01. In the case of DFI, an additive effect for QTLs on CaLG06 ranged from −2.84 to −2.94 (Table 2; Figure 3).

FIGURE 3

Figure 3. Genomic regions with major QTLs for physiological traits. (A) Under heat stress environments of 2017–18 and 2018–19, two QTLs explaining 8.84% and 9.85% PVE, respectively, for normalized difference vegetation index (NDVI) were identified on CaLG02; (B) similarly, in addition to QTLs under two heat stress environments, QTLs for nitrogen balance index (NBI) were also identified based on pooled data on CaLG08; (C) a major QTL explaining 15.04% PVE for chlorophyll content on CaLG05 based on pooled data of under heat stress environments of 2017–18 and 2018–19.

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. Similarly, under HTS pooled data, one major QTL (PVE 15.04%) and one minor QTL (PVE 6.78%) were identified for CHL on CaLG05 flanked by SCA5_30627756–SCA5_1154130 and SCA5_1154130–SCA5_11665932 markers, respectively (Table 2; Figure 3C). Furthermore, under NS 2017, one major QTL (PVE 33.52%) and one minor QTL (PVE 8.26%) were identified for CHL on CaLG04 flanked by SCA4_48715028–SCA4_14861717 and SCA4_48720330–SCA4_48714912 markers, respectively. In addition, one major QTL (PVE 19.71%) and one minor QTL (PVE 9.92%) were identified on CaLG04 based on pooled data under NS.

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).

FIGURE 4

Figure 4. Genomic regions with major QTLs for yield and yield-related traits. (A) Under heat stress environments of 2017–18 and 2018–19 as well as based on pooled data of these two environments QTLs for days to pod initiation (DPI, d) and biological yield per plant (BYPP, g) were mapped on CaLG06. In addition, days to flower initiation (DPI, d) of two QTLs for normalized difference vegetation index (NDVI) were identified explaining 8.84% and 9.85% PVE, respectively, on CaLG02; (B) QTLs for days to pod initiation (DPI, d) and 100 seed weight (100SDW, g) under stress environment of 2017–18 as well as based on pooled data of both heat stress environments of 2017–18 and 2018–19 were co-localized on CaLG07.

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.

100 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%.

Candidate Genes in QTL Regions

Mining of the candidate genes for heat tolerance revealed 1,498 genes in 24.82 Mb (8.68–33.50 Mb) region on CaLG01, 1,162 genes in 23.49 Mb (8.48–31.98 Mb) region on CaLG02, 1,408 genes in 25.71 Mb (14.86–40.57 Mb) region on CaLG04, 140 genes in 4.88 Mb (39.03–43.91 Mb) region on CaLG06, and 2,074 genes in 38.35 Mb (9.56–47.91 Mb) region on CaLG07 (Supplementary 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).

TABLE 3

Table 3. Key heat stress responsive genes in the QTL regions.

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 (Gaur et al., 2014). 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 (Gaur et al., 2015), 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 (Thudi et al., 2014; 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 (Varshney et al., 2014). 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 (Bajaj et al., 2015; 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 (Singh et al., 2017). 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 (Zhang et al., 2015; 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 extensin-like 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.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

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 Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.655103/full#supplementary-material

Supplementary Figure 1. Mean weekly minimum and maximum day temperature (°C) recorded during the crop season of 2017–2018 and 2018–2019. 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.

Supplementary Table 1. Correlations among different traits under normal sown environments of 2017–2018.

Supplementary Table 2. Correlations among different traits under late sown or heat stress environments of 2017–2018.

Supplementary Table 3. Correlations among different traits under normal sown environments of 2018–2019.

Supplementary Table 4. Correlations among different traits under late sown or heat stress environments of 2018–2019.

Supplementary Table 5. Correlations between different traits based on pooled data of 2 years under normal sown environments.

Supplementary Table 6. Correlations between different traits based on pooled data of 2 years under late sown or heat stress environments.

Supplementary Table 7. Genotyping-by-sequencing–single-nucleotide polymorphism (GBS–SNP) marker data used for construction of genetic map for recombinant inbred line population derived from DCP 92-3 (heat sensitive) × ICCV 92944 (heat tolerant).

Supplementary Table 8. Distribution of SNPs on the genetic map derived from DCP 92-3 × ICCV 92944.

Supplementary Table 9. Genes identified in the genomic regions harboring major QTLs on CaLG01, CaLG02, CaLG04, CaLG06, and CaLG07.

References

Agarwal, G., Garg, V., Kudapa, H., Doddamani, D., Pazhamala, L. T., Khan, A. W., et al. (2016). Genome wide dissection of AP2/ERF and HSP90 gene families in five legume sand expression profiles in chickpea and pigeonpea. Plant Biotechnol. J. 14, 1563–1577. doi: 10.1111/pbi.12520

PubMed Abstract | CrossRef Full Text | Google Scholar

Bajaj, D., Saxena, M. S., Kujur, A., Das, S., Badoni, S., Tripathi, S., et al. (2015). Genome-wide conserved non-coding microsatellite (CNMS) marker-based integrative genetical genomics for quantitative dissection of seed weight in chickpea. J. Exp. Bot. 66,1271–1290. doi: 10.1093/jxb/eru478

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhandari, K., Sharma, K. D., Rao, B. H., Siddique, K. H., Gaur, P., Agrawal, S. K., et al. (2017). Temperature sensitivity of food legumes: a physiological insight. Acta Physiol. Plant. 39:68. doi: 10.1007/s11738-017-2361-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhandari, K., Sita, K., Sehgal, A., Bhardwaj, A., Gaur, P., Kumar, S., et al. (2020). Differential heat sensitivity of two cool-season legumes, chickpea and lentil, at the reproductive stage, is associated with responses in pollen function, photosynthetic ability and oxidative damage. J. Agron Crop Sci. 206, 734–758. doi: 10.1111/jac.12433