Unraveling the genetics of heat tolerance in chickpea landraces (Cicer arietinum L.) using genome-wide association studies

Chickpea, being an important grain legume crop, is often confronted with the adverse effects of high temperatures at the reproductive stage of crop growth, drastically affecting yield and overall productivity. The current study deals with an extensive evaluation of chickpea genotypes, focusing on the traits associated with yield and their response to heat stress. Notably, we observed significant variations for these traits under both normal and high-temperature conditions, forming a robust basis for genetic research and breeding initiatives. Furthermore, the study revealed that yield-related traits exhibited high heritability, suggesting their potential suitability for marker-assisted selection. We carried out single-nucleotide polymorphism (SNP) genotyping using the genotyping-by-sequencing (GBS) method for a genome-wide association study (GWAS). Overall, 27 marker–trait associations (MTAs) linked to yield-related traits, among which we identified five common MTAs displaying pleiotropic effects after applying a stringent Bonferroni-corrected p-value threshold of <0.05 [−log10(p) > 4.95] using the BLINK (Bayesian-information and linkage-disequilibrium iteratively nested keyway) model. Through an in-depth in silico analysis of these markers against the CDC Frontier v1 reference genome, we discovered that the majority of the SNPs were located at or in proximity to gene-coding regions. We further explored candidate genes situated near these MTAs, shedding light on the molecular mechanisms governing heat stress tolerance and yield enhancement in chickpeas such as indole-3-acetic acid–amido synthetase GH3.1 with GH3 auxin-responsive promoter and pentatricopeptide repeat-containing protein, etc. The harvest index (HI) trait was associated with marker Ca3:37444451 encoding aspartic proteinase ortholog sequence of Oryza sativa subsp. japonica and Medicago truncatula, which is known for contributing to heat stress tolerance. These identified MTAs and associated candidate genes may serve as valuable assets for breeding programs dedicated to tailoring chickpea varieties resilient to heat stress and climate change.


Introduction
Chickpea (Cicer arietinum L.) is an important annual legume plant with a genome size of 738 Mb (Varshney et al., 2013).India stands as the world's largest producer of chickpeas, contributing to 75% of the global production (Gaur et al., 2019).Recent advancements in genetic resource enhancement for chickpeas involve next-generation sequencing (NGS) initiatives, which have revolutionized genomics research.NGS has enabled cost-effective and rapid sequencing of chickpea genomes, unveiling thousands of genetic markers that hold immense potential for enhancing chickpea breeding (Varshney et al., 2019).Additionally, genomics research plays a pivotal role in identifying the specific genes governing vital traits in chickpeas, including drought, salinity, and heat tolerance (Bhat et al., 2020;Jha et al., 2021;Kumar et al., 2021).
The ongoing impact of climate change has brought significant changes in Indian agriculture, resulting in a shift toward chickpea cultivation in warmer regions (Mannur et al., 2019).This change has exposed chickpea crops to drought and heat stress, leading to yield penalties of up to 70%.In India, an increase in seasonal temperature of 1°C leads to a yield reduction of 474 kg/ha, and high temperatures (≥35°C) can lead to a 39% overall yield reduction in chickpeas (Jain et al., 2023).Also, as per the intergovernmental panel on climate change, the current rate of global warming is 0.2°C per decade and is predicted to touch 1.5°C between 2030 and 2050 (Mohanty et al., 2024).Such rising temperatures will lead to high heat stress and a severe threat to global food security by impacting yield losses of up to 10% to 15% in most of the food crops for every 1°C rise in the optimum temperature (Devasirvatham and Tan, 2018).
The intricate nature of heat stress tolerance, which involves multiple genes with pleiotropic and epistatic effects, presents a substantial challenge in breeding.Phenotypic evaluation in earlier studies has identified some heat-tolerant genotypes in chickpeas (Jha et al., 2022).However, there are no reports of screening chickpea landraces and wild germplasm found near the region of crop domestication [West Asia and North Africa (WANA) region] for heat stress tolerance.Most research on chickpeas to date has been centered on varietal evaluation for yield and growth traits at a specific developmental stage for heat stress tolerance.However, the expanding genomic resources have facilitated the identification of heat stress-responsive quantitative trait loci (QTLs) in chickpeas through QTL mapping and genome-wide association studies, using both mapping populations and natural germplasm collections (Kushwah et al., 2021).However, only a limited number of initiatives have successfully identified specific genes/alleles or significant QTLs that confer heat tolerance across a wide range of genetic backgrounds.Given the critical importance of developing heat stress-tolerant varieties in today's context, the pursuit to identify potential and major genomic loci that regulate heat tolerance using the landraces can be a valuable strategy for heat stress breeding programs.
With this background and aim, the current study focused on the identification of QTLs for heat tolerance in 153 chickpea landraces.These landraces were selected to encompass a broad spectrum of genetic diversity and geographical origins.Our investigation focused on an in-depth exploration of growth traits, spanning from the flowering stage to harvestable maturity, employing robust phenotyping methodologies.Genome-wide association studies (GWASs) were conducted to delve into the interrelationships among various agronomic traits to identify the potential markertrait association for heat stress tolerance.

Plant material and experimental setup
The association panel under study comprised 153 diverse landraces along with one tolerant check variety JG14 in the current investigation, predominantly representing the centers of diversity for chickpeas in the WANA region.These landraces were procured from the chickpea molecular breeding laboratory, Division of Genetics, ICAR-Indian Agricultural Research Institute, New Delhi (Supplementary Table 1).Multi-environment trials (METs) were conducted at three distinct geographic locations, i.e., Amlaha, Dharwad, and Delhi, for the two consecutive years, 2021-22 and 2022-23 (Table 1).The trials encompassed two growing conditions, timely sown and late sown (30 days after the former).The experiment was meticulously designed using an alpha lattice design.Each genotype was replicated three times at each location.All the genotypes were randomly distributed and sown in plots with rows measuring 2 m in length, a plant spacing of 30 cm, inter-row spacing of 50 cm, inter-plot spacing of 1.5 meters, and a 2.0-m gap between replications.All the standard recommended package of practices was followed to raise a good chickpea crop in each location.

Phenotypic observations
The chickpea landraces were extensively phenotyped at all the locations for both the conditions for all the two years for recording the following quantitative traits: days to 50% flowering (DTF; number of days), days to maturity (DTM; number of days), plant height (PH; cm), biomass yield (BY; g/1 m of row), hundred-seed weight (HSW; g), plot yield (g/1 m of row), and harvest index (HI), which was calculated by dividing plot yield by biomass yield.

SNP genotyping
For high-throughput genotyping, DNA was isolated from leaves of 7-day-old seedlings grown under controlled conditions.The DNA extraction procedure followed the cetyl trimethylammonium bromide (CTAB) method originally described by Murray and Thompson in 1980 and modified by Kumar et al. (2013).To assess DNA quality, 0.8% agarose gel electrophoresis was used, and DNA samples measuring 30 ng/µL were selected for subsequent singlenucleotide polymorphism (SNP) genotyping.FASTQ files of all the genotypes having 64-bp sequences with barcode adopters and common adopters were sorted, and counted sequence tags were merged to TagCount files used to align to Reference genome CDC Frontier v1 and were used for TagsOnPhysicalMap (TOPM).Parallelly, Master TagCounts + original FASTAQ files determined the distribution of master tags among samples (taxa) to generate the TagsByTaxa (TBT) file.The combined TagsOnPhysicalMap (TOPM) + TagsByTaxa (TBT) file was used for calling and filtering SNP and updating TOPM with variants followed by the production of a ready TOPM file.Finally, the HapMap Format (hmp) was provided and used for the GWAS analysis after filtering and imputation.Genotyping by sequencing (GBS) was employed to obtain SNP marker data in the HapMap format, which comprised 16,892 SNPs.Subsequently, after filtering for monomorphic alleles and considering criteria such as a minor allele frequency (MAF) of less than 0.05, a missing data frequency exceeding 0.18, and a heterozygote frequency surpassing 0.20, a total of 4,530 SNPs were retained for the subsequent GWAS analysis using TASSEL 5.0 (Bradbury et al., 2007).

Data analysis
The phenotypic data generated in multi-location trials (MLTs) under different conditions were analyzed using the R package agricolae version 1.3-6.The comprehensive analysis includes ANOVA and adjusted means for each genotype based on the alpha lattice design.Principal component analysis (PCA) was carried out using the R package FactoMineR version 2.4.The graphical representation of the PCA results was generated using the R package factoextra version 1.0.7 (Kassambara, 2020).Furthermore, Pearson's correlation coefficients were computed to assess the relationships among the studied traits, and visual representations of these correlations were constructed using the R package corrplot (Wei et al., 2017).
For effective analysis of genotypic data, a total of 4,530 SNPs, evenly distributed at approximately 1-Mb intervals across the genome, were filtered and utilized to assess the population structure using the STRUCTURE version 2.3.4 (Pritchard et al., 2010).The analysis was conducted with the following specific parameters: 100,000 burn-in cycles and 100,000 Markov chain Monte Carlo (MCMC) with three iterations performed for each k value ranging from 1 to 10.The optimal number of subpopulations (delta K) was determined using the Evanno (Evanno et al., 2005) method outlined by Structure Harvester (http://taylor0.biology.ucla.edu/structureHarvester/).In the context of cluster analysis, a distance matrix was generated using TASSEL version 5. Subsequently, a neighbor-joining (NJ) tree file in Newick format was exported to iTOL version 6.5.2 (accessible at https:// itol.embl.de/).This allowed us to create a dendrogram using the neighbor-joining method.PCA and kinship analysis based on SNP markers were carried out using the GAPIT (Lipka et al., 2012).Graphical representation of the genetic position of MTA-QTLs was carried out using MapChart 2.3 (https://www.wur.nl/en/show/mapchart.htm).

Association mapping analysis
We computed the r 2 values for 4,530 pairs of SNP markers, then filtered them focusing on pairs within each chromosome, and created a linkage disequilibrium heat map to identify significant linkage disequilibrium (LD) block and its size, which falls diagonally in the heat map at a p-value of 0.001.We generated a whole genome and sorted for individual chromosomes by utilizing TASSEL version 5. Subsequently, we used these files to generate LD decay curves for all eight chromosomes individually and for the entire genome.To estimate the sizes of LD blocks, we plotted the r 2 values against the distance in base pairs (bp) while setting a threshold at r 2 = 0.2 and recorded the distance at which LD decay reached its midpoint.For the genome-wide association analysis, we utilized the 4,530 SNP markers along with the adjusted mean and best linear unbiased prediction (BLUP) values obtained from META-R (Multi Environment Trial Analysis with R for Windows) version 6.0 in linear mixed models.These models estimated random effects for each trait, and we analyzed them using GAPIT version 3 in R with PCA 3 as a default parameter (Wang & Zhang, 2021).We chose the BLINK (Bayesianinformation and linkage-disequilibrium iteratively nested keyway) model for its computational efficiency and ability to avoid false positives in identifying quantitative trait nucleotides (QTNs) (Huang et al., 2019;Zhou et al., 2023), as it considers computationally efficient fixed effect model (FEM) and avoids computationally expensive random effect model (REM).To evaluate the quality of the association model fitting, we employed a quantile-quantile (Q-Q) plot, which compared the expected and observed −log 10 (p) values.To ensure stringency in our selection of marker-trait associations (MTAs), we applied the Bonferroni correction.This entailed establishing a significance threshold at −log 10 (p) value of 0.05 divided by the total number of markers, which was 4,530, in order to minimize the risk of false positives.We visualized these significant MTAs using a Manhattan plot.We also identified stable MTAs across different locations.Moreover, we pinpointed pleiotropic SNPs, which were associated with multiple traits simultaneously.To further explore the potential candidate genes linked to these significant SNPs, we conducted a search for gene-coding regions within a 100-kb flanking region of the MTAs.We performed this search against the NCBI Reference genome ASM33114v1, which has a size of 530.8 Mb.We utilized the pulse database for chickpea (https://www.pulsedb.org/blast/report/)for annotating these significant QTNs.

Phenotypic variability in chickpea landraces
The pooled analysis of variance indicated a significant difference between the studied traits for genotype, treatment, season, location, genotype by treatment, and treatment by location.Furthermore, we noticed the genotype by treatment by location was non-significant for DTF, DTM, and HSW (Table 2).The frequency distribution curve indicated the normal to near-normal distribution for the majority of the traits under normal timely sown and late-sown conditions (Figure 1; Supplementary Figure 1).We performed the Shapiro-Wilk test to prove this statement.p-Value ≥0.05 is considered a normal distribution, ≥0.001 is a near-normal distribution, and ≤0.001 is a non-normal distribution.
The mean values for the traits exhibited significant variations across locations.Under normal sown conditions across locations, environmental variation ranges from 4.73 for DTM_DN2021 to 82,184.86 for the (BY_DeN2022) trait.The genotypic variation ranges from 3.05 for DTM_DeN2021 to 76,340.1 for the BY_DN2021 trait.While phenotypic variation ranges from 18.01 for HSW_DeN2021 to 140,100.4 for BY_AN2022.The broad-sense heritability ranges from 7.75% (HI_DeN2022) to 95% (DTM_DN2021), and genetic advance ranges from 1.21 (DTM_DeN2021) to 501.94 (BY_DN2021).Notably, under latesown conditions, we observed the highest phenotypic variation for traits such as biological yield, plot yield, and harvest index.This was coupled with broad-sense heritability figures ranging from 13.24% (for BY_DL2022) to 93% (for DTM_DL2021).We observed that the genotypic mean square for HI_DeN2022 was non-significant, but all the remaining traits show significant differences at p < 0.001 (Supplementary Table 2).

Phenotypic correlation and PCA
Pearson's correlation coefficient calculated for the traits for both conditions indicated a positive correlation among the traits such as plot yield (PY), BY, and HI across various locations.However, PY showed a negative correlation with DTM (p < 0.05) and exhibited non-significant correlations with PH, DTF, and HSW under normal sown conditions across different locations.In the late-sown condition, we observed negative correlations between DTF and DTM with PY, while PY showed non-significant correlations with HSW (Supplementary Figure 2).The PCA under the normal sown conditions of Amlaha in 2021 showed that the first principal component explained 27.54% of the variation.The primary contributors to this component were BY and PY.Meanwhile, the second dimension, which accounted for 24.7% of the variation, was influenced by traits like PH and DTF, with additional contributions from HSW and HI.In the third dimension, HSW and HI played significant roles, contributing to the variation.DTM had an impact on the fifth dimension.Notably, BY and PY were closely clustered together in the analysis, forming an acute angle, which indicated a positive correlation between these two traits.However, HSW, DTF, DTM, and PH clustered together at an acute angle but had a straight angle with HI, suggesting a negative correlation between HI and HSW, DTF, DTM, and PH traits.Similarly, we observed that DTM and DTF were negatively correlated to HI and PY where an angle of 180°between coordinate lines contributed to PC1 of 39.6%, similar to BY and HSW under the late-sown condition of Amlaha 2021.PCA in Amlaha 2022 indicated that dimension 1 explains 26.6% and 38.1% with dimension 2 explaining 21.8% and 21.1% for normal and late-sown conditions, respectively.In Delhi 2021, under controlled conditions, dimension 1 accounted for 37.9% of the variance, and dimension 2 for 21.3%.For late conditions, dimension 1 explained 34.9%, and dimension 2 explained 17.9%.In Delhi 2022, for timely sowing, dimension 1 explained 29.5%, while dimension 2 explained 23.8%.For late sowing, dimension 1 explained 36.9%, and dimension 2 explained 18.8%.For Dharwad 2021, under controlled conditions, dimension 1 contributed to 27.2% of the variance and dimension 2 to 24.8%.For late conditions, dimension 1 explained 36.5%, and dimension 2 explained 22.4%.In Dharwad 2022, for timely sowing, dimension 1 accounted for 29.9%, and dimension 2 for 21.9% of the variance.For late sowing, dimension 1 explained 37.4%, and dimension 2 explained 19.2% (Supplementary Figure 3).

Analysis of SNP marker distribution, population diversity, and linkage disequilibrium
Population structure analysis revealed the best DK vs. K-value to be 2, indicating the presence of two distinct subpopulations within the GWAS panel (Figures 3A, B).Subpopulation 1 comprised 33.33% (51) of total genotypes, while Subpopulation 2 consisted of 50.32% (77) of total genotypes, and 16.33% (25) of total genotypes were classified as part of the admixture population (Supplementary Table 3).This finding was corroborated by PCA based on SNP marker data, which also depicted two distinct clusters (Figure 3C), reinforcing the presence of two subpopulations.Additionally, kinship and neighbor-joining cluster analyses also supported the presence of these two clusters (Figures 3D, E).LD between marker pairs was determined using r 2 values, and LD decay plots were created by plotting r 2 values against genetic distance in base pairs (bp).A significant LD block size of 0.14 Mb (140 kb) was   4).

Genome-wide association study
For the GWASs, BLUP values were analyzed for individual traits and conditions for both the years separately and combined across locations, years, BLUPs, and the overall BLUPs across treatments of locations.A total of 75 unique MTAs were identified above ≥−log 10 (4.00) for all traits at a cut-off p-value of <0.05 as a Bonferroni correction for stringent selection [−log 10 (p) > 4.95].After this correction, only 27 highly significant MTAs (Table 4; Supplementary Table 4) were found.Out of 27 MTAs, 10 MTAs for BY, one for DTF, six for DTM, two for HI, five for HSW, and three for PY were retained, and they were depicted using Manhattan and Q-Q plots (Figure 5; Supplementary Figures 4A, B  physical location with −log 10 (p) ranges from 5.06 to 6.44 under Dharwad 2022, Amlaha 2021, Dharwad 2022(PY) and by acrosstreatment BLUP (Supplementary Figure 4B).Similarly, pleiotropic effects and cross-confirmations were noted for SNPs Ca7:41673233 and Ca8:10963827 (Supplementary Table 4), while chromosomes 2, 3, and 6 had a maximum number of MTAs (6) to lowest (1) on chromosome 5 (Figure 6; Supplementary Figure 5).

Allelic effects of identified genomic regions on respective phenotypes
Twenty-seven major MTAs were analyzed to determine the range of phenotypic variations for all traits (Figure 7; Supplementary Figure 6).Boxplots display the phenotypic values associated with reliable QTNs that exhibit significant effects (p < 0.01) on their Whole-genome and individual chromosome-wide linkage disequilibrium (LD) decay in the GWAS panel.GWAS, genome-wide association study.Danakumara et al. 10.3389/fpls.2024.1376381Frontiers in Plant Science frontiersin.orgcorresponding traits.The landraces should be categorized into two groups based on whether they carry the superior or inferior allele for each QTN.The x-axis of the plot represents the two allele types for each QTN, while the y-axis depicts the phenotypic values.Association panel genotypes were divided into two classes according to allele types.It was observed that all 27 QTNs demonstrated a significant effect on respective traits (p ≤ 0.01).Among these, marker Ca2:2311917 at Chr2 was pleiotropic for HSW_DN2021, HSW_AN2022, and HSW_AN.These significant associations suggest their plausible role in determining heattolerant traits.

Putative candidate genes associated with MTAs
After rigorous selection criteria were applied, a total of 27 MTAs were identified and retained.To explore potential candidate genes associated with these MTAs, a search was conducted within a 100-kb region flanking each MTA considering LD decay of 140 kb.This search utilized sequence information obtained from the SNP markers and their genome position, which were identified through a BLAST search in PulseDB (Pulse Database) against the respective individual chromosomes of chickpeas.The identified marker Ca1:270126601 for  Manhattan and respective quantile-quantile (Q-Q) plots of significant associations for studied traits using individual BLUPs.BLUP, best linear unbiased prediction.

Discussion
Chickpea is an important dietary legume crop of arid and semiarid regions and holds a prominent position as an economical source of protein-rich food.However, the potential yield of chickpeas faces significant challenges due to high-temperature stress at various growth stages.This stress has adverse effects on factors such as pollen viability, pollen germination on the stigma, pollen tube growth, and poor pod formation, ultimately affecting its yield potential.To combat these problems, it is imperative to gain a thorough understanding of the genomic regions that influence the heat-tolerant abilities of chickpeas.This knowledge is essential for developing specific varieties that are highly resilient to heat stress while exhibiting higher yields.In earlier studies, QTL-hotspot regions associated with drought tolerance and related traits were successfully introgressed to enhance the drought resistance of elite crop varieties through marker-assisted backcross (MABC) breeding (Bharadwaj et al., 2021).For such purposes, the identification of markers linked to the target traits is a fundamental requirement.Therefore, in this study, we sought to identify MTAs, linked to yield and yield-related traits under heat stress conditions, using germplasm collected from the WANA region.
Analysis of variance for all the studied traits (except HI under normal sown conditions in 2022) in each year indicated significant variability among these traits.This variability is a crucial prerequisite for conducting genetic studies and implementing breeding programs.Notably, there was significant genetic variability observed for the tested traits under both normal and heat stress conditions in both years (Jha et al., 2022).Higher trait values were observed in normal sown conditions as compared to late-sown conditions, such as DTF and DTM with the lowest mean values as stable across seasons and locations and higher values for all other traits studied considered under this study; for example, genotype IG5866 was stable under normal timely sown condition, or PY and ILC8666 were stable under late-sown condition (Supplementary Table 6).It is important to consider that PY and BY are complex traits characterized by a quantitative pattern of inheritance and susceptibility to environmental influences, including season and treatment due to soil nutrient conditions.Furthermore, when evaluating phenotypic characteristics, it was observed that the phenotypic coefficient of variation (PCV) and broad-sense heritability were notably higher for PY and BY when compared to the HI.Among the studied traits, DTM consistently exhibited high heritability across all 12 environments, followed by Allelic effects of selected MTAs identified in multiple locations for the studied traits under study.MTAs, marker-trait associations.
HSW (Hussain et al., 2021).It is worth noting that the heritability of PY varied from low to medium, spanning a wide range from 16.66% in Delhi under normal conditions in 2022 to 76% in Delhi under late conditions in 2021 (Paul et al., 2018).A significant positive correlation between grain yield and high heritability along with high genetic variability and fewer yield losses under optimal conditions are essential for a characteristic to be expressed as a heat-tolerant marker (Maqbool et al., 2017;Chandora et al., 2020).Interaction ANOVA showed a significant variation (p < 0.001) for all tested traits by considering genotype, treatment, season, replication, location, replication by block, genotype by location, genotype by treatment, treatment by location, and genotype by treatment by location.However, DTF, DTM, and HSW were non-significant for genotype by treatment by location and DTM for genotype by location, as these traits showed higher broad-sense heritability as near qualitative nature of traits, and traits DTM and HI were nonsignificant for replication within the block due to lower phenotypic variation for these traits influenced less by soil heterogeneity across location and treatments.
It was observed that Pearson's correlation coefficient was found positive between PY, BY, and HI under all the treatments across the locations.Notably, PY showed a non-significant correlation with HSW and negative in Delhi 2022 late condition along with PH p < 0.05.This is because the association mapping (AM) panel predominantly consists of kabuli type, which is bold seeded, hence increasing the weight of seeds, and the costing number of seeds leads to independent contribution toward plot yield.Sandhu et al. (2012) observed a negative correlation that is due to the effect of heat stress hampering seed weight and ultimately yield.PCA shows that the component traits such as PY, BY, and HI predominantly contributed to PC1.These traits exhibited parallel trends and clustered together consistently across various locations and treatment conditions.Traits influenced by additive gene action and showing positive correlations can be collectively and efficiently improved irrespective of environmental influences (Borah et al., 2018).
GWAS panels were grouped into two sub-groups with 25 admixtures, clearly indicating the WANA region's representativeness and relatedness.The LD decay over genetic distance in a population indicates the requirement of a marker density to capture the markers close enough to the causal loci (Flint-Garcia et al., 2003;Danakumara et al., 2021).A large LD block size of 0.14 Mb was found for the whole genome.However, a slow rate of LD decay was observed for chromosome 6 at 0.30 Mb followed by chromosomes 3 and 4 at 0.16 Mb size, whereas faster LD decay on chromosome 5 of 0.01 Mb size was observed (Figure 4) (Sokolkova et al., 2020).A higher LD was observed in chickpeas due to low effective recombination rates as compared to cross-pollinated crops (Niu et al., 2019;Samineni et al., 2022).The extent of LD decay in the association panel of chickpeas was observed at 200-300 kb (Basu et al., 2019) and 5 cM in the chickpea reference set (Thudi et al., 2014).The extent of LD can vary due to the complexity and size of the genome and marker number (Valdisser et al., 2017).The LD may vary in different populations because of population size, genetic drift, admixtures, selection, mutation, nonrandom mating, mode of pollination, gene-rich region, and recombination frequency (Vos et al., 2017).
A genome-wide association study was conducted using the BLINK model within the GAPIT package, which is considered superior for identifying QTNs and minimizing the occurrence of false positive associations as the most reliable method (Huang et al., 2019).A total of 27 significant MTAs were identified with a p-value threshold of <0.05.This rigorous approach was taken to increase the stringency of selection and reduce the likelihood of false positive results, and it involved applying a Bonferroni correction.In this analysis, a total of 22 individual markers were found to be significantly associated with the traits of interest (Kaler and Purcell, 2019).Recent advancements in the sequencing and annotation of the chickpea genome have provided valuable insights into the identification of candidate genes within the genomic regions pinpointed through GWAS.These candidate genes are believed to play a role in modulating the variation observed in heat-responsive traits (Srungarapu et al., 2022).We compiled a table of candidate genes located within 100-kb regions surrounding the identified MTAs, which included information such as their start and stop positions, InterPro IDs, gene ontology (GO) terms, and accessions, as well as the probable proteins they encode.In the genomic regions associated with linked SNP markers for PY, an MTA on chromosome 3 at 15.5 Mb was previously reported by Kalve et al. (2022).However, the remaining MTAs were novel findings specific to our study.For instance, SNP marker Ca1:27012660, linked to BY_AL2022, was located near candidate regions encoding a sucrose non-fermenting 4-like protein ortholog of A. thaliana.Through in silico analysis, we identified coding regions in the vicinity of Ca1:27012752 to Ca1:27025275 on chromosome 1.These regions also encoded a pentatricopeptide repeat-containing protein orthologous to A. thaliana.This protein has been recognized as important for resistance against both biotic and abiotic stresses in Oryza sativa (Qiu et al., 2021).Moreover, mQTL-seq analysis revealed a functionally relevant candidate gene within a major QTL region governing pod number in chickpeas.This gene belongs to the pentatricopeptide repeat (PPR) family and exhibits over 80% sequence conservation with its Arabidopsis ortholog, At1g52620 (Das et al., 2016).Additionally, our investigation identified genes such as pentatricopeptide repeatcontaining proteins and serine threonine protein kinases, both of which have well-established roles in seed development across various crop species (Li et al., 2014).These findings shed light on the potential involvement of these genes in stress management and nutrient content regulation within chickpea grains.
The MTA associated with DTM for the year 2021 is linked to marker Ca2:18671666.This marker likely encodes the indole-3-acetic acid-amido synthetase GH3.1 ortholog, which has similarities to its counterparts in A. thaliana and Solanum lycopersicum.Furthermore, expression profiling studies have indicated the diverse roles of GH3 genes in various aspects of development and responses to abiotic stress in leguminous crops like chickpeas.Notably, genes such as CaGH3-3 in chickpeas;  in Lotus were found to be upregulated in root tissues.This suggests their potential involvement in root development processes.Moreover, certain GH3 genes, specifically CaGH3-1 and CaGH3-7 in chickpeas, as well as MtGH3-7, MtGH3-8, and MtGH3-9 in Medicago, exhibited high levels of induction under conditions of drought and/or salt stress (Singh and Jain, 2014).A common MTA for HSW_HSW_DN2021, HSW_AN2022, and HSW_ANP is a putative disease resistance protein ortholog of Solanum bulbocastanum and M. truncatula with the marker location ranging from Ca2.2322855 to Ca2.2323917, which is also encoded to deoxynucleoside triphosphate triphosphohydrolase, SAMHD1 homolog, ortholog of Dictyostelium discoideum, and G. max effects on the control of cell proliferation and apoptosis (Felip et al., 2022).Trait DTM_AL2021 with marker Ca2:31252891 association encodes cation/H( + ) antiporter 15 ortholog of A. thaliana and M. truncatula role as universal stress protein domain containing droughtresponsive genes in pigeon pea (Cajanus cajan L.) (Sinha et al., 2015).Trait BY_AN2021 with marker Ca3:10159944 encodes probably inactive leucine-rich repeat receptor-like protein kinase role (LRR-RLK) as represents the largest group of RLKs in plants and plays vital roles in plant growth, development, and the responses to environmental stress (Song et al., 2022).Trait BY_AL2022 of SNP Ca3:171579 encodes to tetraketide alpha-pyrone reductase 1 ortholog of A. thaliana and dihydroflavonol-4-reductase as ortholog of M. truncatula and has a role in fatty acyl-CoA ester synthesis.Trait BY_AL is associated with marker Ca3:23273262 encoded by ATPdependent helicase and HI_N marker Ca3:37444451 that encodes for aspartic proteinase ortholog sequence of O. sativa subsp.japonica and M. truncatula as reported recently by Varshney et al. (2019) using a dataset of 3.65 million SNPs derived from the resequencing of 429 chickpea germplasms collected globally.Additionally, several potential candidate genes were highlighted in their investigation, including TIC, REF6, aspartic protease, cc-NBS-LRR, and RGA3.These genes were found to play crucial roles in conferring tolerance to both heat and drought stress conditions.Furthermore, the identified MTAs and genomic regions were recognized as stable controllers of traits related to pods per plant, yield, and phenological traits (Rani et al., 2020).Trait DTM_AL2021 with SNP Ca3:39084979 marker linked to AT-hook motif nuclear-localized protein 14 role as AT-hook motif nuclear localized (AHL) gene family is a highly conserved transcription factor critical for the growth, development, and stress tolerance of plants (Zhang et al., 2023).Marker Ca4:5907421 for trait HI_AN2022 codes likely proteins as dormancy-associated protein homolog 3 ortholog of A. thaliana and M. truncatula, as its over-expression of glucose-6-phosphate/phosphate translocator 2 (GPT2) was investigated in chickpea leaves with roles as heat tolerance.Also, transgenic A. thaliana over-expressing metallothionein 1 (MT1) gene of desi chickpea was subjected to transcriptome analysis, and it drought tolerance was evaluated in 7-day-old plants (Kumar et al., 2022).Trait PY_DN2021 with SNP Ca4:6352125 relates to glycogen synthase kinase-3 homolog and shows genome-wide identification and expression analysis of glycogen synthase kinase encoding genes in foxtail millet (Setaria italica L.) under salinity, dehydration, and oxidative stress (Singh et al., 2023).Marker Ca4:8669498 for HSW_DN2021 and HSW_AN2022 codes for probable disease resistance RPP8-like protein 4. Similarly, trait DTM_AL2021 with marker Ca5:40828566 codes for methionine gamma-lyase ortholog sequence of A. thaliana and cystathionine gamma-lyase of ortholog M. truncatula; methionine gamma lyase maintains the equilibrium of isoleucine in a variety of plants under different environmental conditions such as drought (Joshi and Jander, 2009;Khan et al., 2019).Pleiotropic marker Ca6:10230657 for traits BY_N, BY_DN2022, BY_AN2021, and PY_DN2021, which codes serine/ threonine-protein phosphatase PP2A catalytic subunit, is involved in several physiological responses in plants, playing important roles in developmental programs, stress responses, and hormone signaling.Six PP2A catalytic subunits (StPP2Ac) were identified in cultivated potato (Muñiz Garcıá et al., 2022).Trait BY_AN2021 associated with marker Ca6:9109096 codes for magnesium transporter MRS2-4A rootexpressed magnesium transporter of the MRS2/MGT gene family in A. thaliana and allows for growth in low-Mg 2+ environments (Gebert et al., 2009).Also, trait BY_AN2021 with marker Ca6:9109096 codes for 2Fe-2S ferredoxin-type domain participates in numerous biological processes, including carbon fixation, nitrogen assimilation, chlorophyll metabolism, and fatty acid synthesis, and it contributes to plant resilience against heat stress (Meng et al., 2022).Traits PY_DN2022, BY_N, and BY_DN2022 with marker Ca7:41673233 codes for bgalactosidase-like protein ortholog of M. truncatula-related proteins, including b-galactosidase, glucanase, sucrose synthase, cystathionine gamma-synthase, 1-aminocyclopropane-1-carboxylic acid oxidase, abscisic acid b-glucosyltransferase, and late embryogenesis abundant proteins, which all impart heat stress tolerance in chickpeas (Parankusam et al., 2017).Ca7:41673233 also codes for phosphatidylinositol-4-phosphate 5-kinase, and its role as heat shock triggers phospholipid-based signaling pathways (Yadav et al., 2017).Traits DTM_DL2021 and DTM_AL2022 of marker Ca8:10963827 code for pentatricopeptide repeat-containing protein ortholog sequence of A. thaliana and M. truncatula function as exposure to heat stress; the gene expression pattern was significantly altered in relation to various key factors such as heat shock proteins (HSPs), ubiquitin-protein ligases, transcription factors, and pentatricopeptide repeat-containing proteins.Notably, several genes were found to exhibit significant changes in their expression levels.Specifically, genes encoding heat shock proteins (CL2311.Contig3 and CL6612.Contig2), cytochrome P450 enzymes (CL4517.Contig4 and CL683.Contig7), and basic helix-loop-helix transcription factors (bHLH TFs) (CL914.Contig2 and CL8321.Contig1) showed distinctive induction patterns following 4 days of exposure to heat stress (Wang et al., 2019).The candidate regions we investigated contained a significant number of genes, most of which have been recognized for their vital contributions to plant development.Many of these genes have closely related counterparts in other species, and extensive research on these counterparts has provided valuable insights into their functions.Our analysis, which included assessing gene expression and applying gene ontology, pinpointed specific genes with elevated expression levels, shedding light on their roles within cellular organelles.To gain a deeper knowledge of this genomic region and its potential significance, future research could delve into detailed investigations of each identified region.

Conclusion
Past efforts have successfully utilized drought-related QTLhotspot regions to enhance drought tolerance in elite chickpea varieties through marker-assisted backcross breeding.To achieve this, the identification of markers linked to the trait of interest is crucial.In this study, we employed a diverse mapping panel to identify MTAs related to yield and yield-related traits under heat stress.Our analyses revealed significant genetic variability for the tested traits under both normal and heat stress conditions, laying the foundation for genetic studies and breeding programs.Phenotypic coefficients of variation and broad-sense heritability were relatively high for yield-related traits, indicating their suitability as potential heat tolerance markers.Additionally, we observed significant interactions between genotypes, treatments, seasons, locations, and other factors, reflecting the complex nature of these traits.Pearson's correlation coefficients indicated positive relationships between yield traits and harvest index across most treatments and locations.However, we observed some negative correlations due to specific seed characteristics and the impact of heat stress on seed development.Principal component analysis highlighted the importance of yield-related traits in explaining variation across locations and treatments.The LD analysis demonstrated a slow LD decay across the chickpea genome, indicating extended LD blocks and low effective recombination rates, which is a characteristic of self-pollinated crops.Through GWASs using the BLINK model, we identified 27 significant MTAs linked to various traits associated with heat stress tolerance.These MTAs can serve as valuable targets for further research and breeding efforts.Notably, we identified candidate genes near these MTAs, and we further explored candidate genes situated near these MTAs, shedding light on the molecular mechanisms governing heat stress tolerance and yield enhancement in chickpeas such as indole-3-acetic acid-amido synthetase GH3.1 with GH3 auxin-responsive promoter and pentatricopeptide repeat-containing protein.These identified MTAs and associated candidate genes serve as valuable assets for breeding programs dedicated to crafting resilient chickpea varieties in the face of climate change.Trait HI associated with marker Ca3:37444451 encodes aspartic proteinase ortholog sequence of O. sativa subsp.japonica and M. truncatula, contributing to heat and drought tolerance as stable significant MTAs/genomic regions controlling yield trait.These genes have the potential to be integrated into the well-performing but heat-and drought-sensitive popular chickpea cultivars, which ultimately contribute to improved chickpea crop performance under adverse environmental conditions.The insights gained pave the way for future breeding endeavors, ultimately contributing to global food security.

FIGURE 1
FIGURE 1Frequency distribution plot yield (g) (PY) conducted under timely and late conditions in consecutive years 2021 and 2022.
FIGURE 3 Population in the GWAS panel from model.(A) Population structure-based grouping of genotype from STRUCTURE analysis.(B) K vs. DK of structure harvest.(C) The 2D plot of the principal component-based grouping.(D) Neighbor-joining tree-based diversity.(E) Heat map of pairwise kinship matrix.GWAS, genome-wide association study.

FIGURE 2 SNP
FIGURE 2 SNP density plot indicating the distribution of filtered SNPs across chromosomes.SNP, single-nucleotide polymorphism.

FIGURE 6
FIGURE 6Distribution and position (in Mb) of identified marker-trait associations (MTAs) at their respective chromosome with associated traits over the seasons on chickpea chromosomes for plot yield (black), BY (green), DTF (brown), DTM (red), HI (purple), and HSW (cyan).BY, biomass yield; DTF, days to flowering; DTM, days to maturity; HI, harvest index; HSW, hundred-seed weight.

TABLE 1
Details of sowing conditions, locations, and year of experiment along with abbreviations.

TABLE 2
ANOVA mean sum of squares of each trait tested under study.