Genome-wide association analysis of stress tolerance indices in an interspecific population of chickpea

Chickpea is a cool season crop that is highly vulnerable to abiotic stresses such as heat and drought. High temperature during early flowering and pod development stages significantly reduces the crop yield. The wild relatives of chickpeas can be potential donors for the introgression of heat and drought tolerance into cultivated chickpeas for crop improvement. Initially, 600 interspecific lines were derived from crosses between two elite cultivars, CDC Leader (kabuli chickpea) and CDC Consul (desi chickpea), and 20 accessions of Cicer reticulatum. The F₅ interspecific lines were tested for agronomic and seed quality traits including reaction to ascochyta blight disease under field conditions at two locations in 2018. A subset of 195 lines were selected based on resistance to ascochyta blight and acceptable seed quality. These lines were evaluated for their performance under suboptimal conditions at Lucky Lake (2019 and 2020) and Moose Jaw (2019), Saskatchewan, Canada, and Yuma, Arizona, United States (2019–2020). The lines were grown and evaluated at two seeding dates, normal (SD1) and late (SD2) seeding dates, at each location and year. The same lines were genotyped using Cicer60K Axiom® SNP chip. The population structure was determined based on 35,431 informative SNPs using fastStructure, and the interspecific lines were clustered at a k-value of 15. Significant marker-trait associations were identified for seed yield from SD1 and SD2 seeding dates, and stress tolerance indices (ATI, K₁STI, MP, SSPI, and TOL) using phenotypic values both from individual locations and combined analyses based on BLUP values. SNP marker Ca2_34600347 was significantly associated with yield from both the seeding dates. This and other SNP markers identified in this study may be useful for marker-assisted introgression of abiotic stress tolerance in chickpea.

Introduction

Chickpea (Cicer arietinum L.) is the world’s second most important pulse crop after the common bean (Varshney et al., 2013; Rani et al., 2020). It was one of the earliest domesticated legume crops and is currently grown in 59 countries. In 2019, the world production of chickpeas was around 14.2 million tons (FAO, 2019). Among abiotic stresses, drought and heat are the major environmental constraints limiting chickpea production worldwide in recent years (Devasirvatham and Tan, 2018; Rani et al., 2020). It has been reported that heat and drought can cause more than 70% yield loss in chickpea (Varshney et al., 2019). Chickpeas are mostly grown under rainfed conditions without irrigation. Therefore, soil moisture deficit toward the end of the crop season (terminal drought) affects about two-thirds of the global chickpea area (Gaur et al., 2008, 2019). Moreover, being a cool season food legume, chickpea yield is sensitive to heat stress exposure during the reproductive stage (Devasirvatham et al., 2015; Gaur et al., 2015).

Extreme heat and dry conditions are among the main abiotic stresses that affect crop yield across Canada. The prolonged heat wave and lack of precipitation in recent years in Western Canada¹ has had an adverse impact on chickpea yield. Tolerance to abiotic stresses such as heat and drought is a complex trait that is the result of various morphological, physiological, and biochemical changes in plants (Kaloki et al., 2019). Moreover, these abiotic stresses are quantitatively inherited with a large effect of genotype x environment interaction (Jha et al., 2014). The development of cultivars with abiotic stress tolerance and yield stability is critical in chickpea breeding programs. However, very slow progress has been made in developing tolerant cultivars due to the physiological and genetic complexity of the trait. The variability and unpredictability of stress conditions during trials limit the selection efficiency. Therefore, knowledge of the traits responsible for the adaptation of chickpea to suboptimal environments is important for the development of cultivars with improved abiotic stress tolerance.

Crop wild relatives (CWR) preserve higher levels of genetic diversity as they have been challenged in natural environments for many years in comparison to domesticated cultivars. Hence, these wild relatives are crucial genetic resources used by plant breeders for crop improvement (Hajjar and Hodgkin, 2007). CWR has been used as a source of abiotic stress tolerance in many cultivated species (Hajjar and Hodgkin, 2007; Zhang et al., 2019). Recent advances in genotyping, breeding, and genomics have accelerated the use of CWR for crop improvement by marker-assisted introgression of wild alleles into cultivated germplasm (Bohra et al., 2021). Few examples included enhanced drought tolerance in cultivated germplasm of sunflower (Hussain et al., 2019), and improved drought related traits such as water use efficiency, earliness, and yield of cultivated groundnut by introgression of alleles from the wild groundnut species Arachis duranensis and Arachis batizocoi (Dutra et al., 2018). Another example is enhanced drought resistance and productivity of elite durum (Triticum turgidum ssp. durum) and bread wheat (Triticum aestivum L.) cultivars by the reintroduction of alleles from wild emmer wheat (Merchuk-Ovnat et al., 2017). More studies that have demonstrated wild progenitors as a valuable source for the enrichment of the domesticated gene pool for abiotic stress tolerance included reintroducing wild alleles in lentils (Gorim and Vandenberg, 2017) and wheat (Placido et al., 2013).

Genome-Wide Association Study (GWAS) provides a higher resolution of marker-trait association than classical QTL analysis using bi-parental populations (Korte and Farlow, 2013). Both GWAS and QTL mapping studies were used to identify the genetic loci associated with various abiotic stresses including drought and heat stress (Jha et al., 2021). Nested association mapping (NAM) has been a valuable approach to dissecting the genetic architecture of complex quantitative traits (Gangurde et al., 2020; Lakmes et al., 2022). NAM population consists of multiple families of recombinant inbred lines (RILs) derived from multiple inbred lines crossed to a single reference inbred line (Yu et al., 2008). It utilizes the combined power of QTL mapping and association mapping to identify the trait-associated markers (Buckler et al., 2009). A publicly available collection of wild chickpeas, especially Cicer reticulatum, which survives under suboptimal environment, is an important resource to improve stress tolerance in current chickpea cultivars (Singh et al., 2008; von Wettberg et al., 2018). Linkage drag is a known bottleneck for the introgression of QTLs from wild accessions to cultivated. To overcome the linkage drag and to introgress and expand the genetic basis of cultivated chickpea, introgression of wild alleles from multiple accessions of C. reticulatum into cultivated chickpea germplasm has been initiated at the University of Saskatchewan chickpea breeding program. The main objective of the current study was to examine the performance of the interspecific population derived from C. arietinum x C. reticulatum crosses under suboptimal conditions and to identify the genetic loci associated with the traits crucial for plant performance under suboptimal environments using genome-wide association analysis.

Materials and methods

Plant materials

A chickpea interspecific population consisting of 600 lines derived from crosses between elite cultivars (C. arietinum) and 20 accessions of Cicer reticulatum were developed. The elite cultivars were CDC Leader (kabuli chickpea) and CDC Consul (desi chickpea). The design of this population is like a nested association mapping (NAM) design. Each of the 20 C. reticulatum accessions (tester lines) was crossed with CDC Leader and CDC Consul (founder lines). A single F₁ seed from each cross was grown in the greenhouse in 2014 and selfed to produce the F₂ seed. The F₂ plants of each cross were advanced to F₅ as single seed descents (SSDs). The seeds of each F₅ line were bulked and used for the evaluation of their agronomic performance including their reaction to ascochyta blight disease and seed visual quality at two locations, Limerick and Lucky Lake, SK during the 2018 growing season. One hundred and ninety-five F₇ interspecific lines with improved resistance to ascochyta blight and acceptable visual seed quality were selected and used in the current study. All the selected lines were derived from crosses with CDC Leader as a common parent.

Growing conditions and phenotypic data analysis

The 195 F₇ interspecific lines were grown at Lucky Lake (2019 and 2020) and Moose Jaw (2019), Saskatchewan, Canada, and Yuma, Arizona, United States (2019–2020). The population was planted at two different seeding dates [normal (SD1) and late seeding (SD2)] at each location and year. The purpose of late seeding was to expose the plants to higher temperatures during flowering. The normal seeding date in Saskatchewan was in the second week of May (14 May in Moose Jaw, 2019, 17 May in Lucky Lake, 2019, and 12 May in Lucky Lake, 2020) while the late seeding was planted around 2 weeks after the normal seeding (28 May in all the locations; Supplementary Table 1). At Yuma, Arizona, the first seeding was on 6 November while the second seeding was on 14 January (Supplementary Table 1) to maximize the chance that the population was exposed to a temperature above 27°C during flowering. At each location, the lines were planted in a one square meter plot. In each plot, the seeds were planted in three rows with a density of 60 seeds m⁻². A randomized complete block design (RCBD) with three replications was used at each location and seeding date. Daily maximum and minimum temperatures and other meteorological data were recorded at each location. There were multiple occurrences of temperatures above 26°C–27°C in all the locations during flowering on late seeding dates (Supplementary Table 1).

Days to flowering (DTF) were recorded for each line on the plot basis when 50% of the plants had flowered. Days to maturity (DTM) was also documented for each line when 50% of the plants in each plot were matured. Plant height was measured from the ground level to the tip of the plants when the pods reached physiological maturity. Flower color was recorded based on visual observation during flowering. Seed weight was measured by weighing 1,000 seeds per line after harvest. Reaction to ascochyta blight of each line was recorded on a plot basis, using a mixed quantitative and qualitative 0–9 score scale as described by Chongo et al. (2004). Seed yield was measured on a plot basis. Seed yield under stressed (late seeding) and non-stressed (normal seeding) conditions were used to predict the stress tolerance index of each line. The different stress indices used to assess stress tolerance included tolerance index (TOL), mean productivity (MP), abiotic tolerance index (ATI), stress susceptibility percentage index (SSPI), and modified stress tolerance index (K₁STI). These have been reported as reproducible indices under severe stress conditions in chickpeas (Farshadfar and Geravandi, 2013). The following formulas were used to calculate the stress indices.

TOL = Yp–Ys (Rosielle and Hamblin, 1981).

MP = (Yp + Ys)/2 (Rosielle and Hamblin, 1981).

ATI = [(Yp–Ys)/(Y̅p/Y̅s)] × 100 (Moosavi et al., 2008).

SSPI = [(Yp–Ys)/(2Y̅p)] × 100 (Moosavi et al., 2008).

K₁STI = (Yp2/Y̅p2) × [(Yp + Ys)/Y̅p2)] (Farshadfar and Sutka, 2003).

In the above formulas, Ys, Yp, Y̅s, and Y̅p represent seed yield in stress and non-stress conditions for each genotype, and mean seed yield in stress and non-stress conditions for all genotypes, respectively.

For field evaluation, both combined analyses across locations and years for each seeding date, and separate analyses on each year (2019 and 2020) for each seeding date were conducted. ANOVA was done using PROC MIXED in which genotypes were considered as a fixed factor and years as a random factor. The LSMEANS statement was used to compute the average phenotypic score for each line. For a separate analysis of each year, ANOVA was done using the PROC MIXED procedure, in which the lines were considered as a fixed factor and replication was considered as a random factor. To estimate the broad sense heritability (H²), variance components were calculated using the SAS PROC VARCOMP procedure (SAS Institute Inc, 1999). The H² of various phenotypic responses at plot level based on individual experiments and over the years were estimated using the following two equations, respectively

$H^2=\frac{{\sigma }^{2}G}{{\sigma }^{2}G+{\sigma }^{2}er} \mathrm{and} H^2=\frac{{\sigma }^{2}G}{{\sigma }^{2}G+{\sigma }^{2}GY+\sigma er}$

where σ²G, σ²Y, σ²GY, and σ²er are estimates of genotype, site-year, genotype by site-year interaction, and error variance, respectively (Singh et al., 1993). Spearmen’s rank correlation coefficients between the seed yields obtained under stress and non-stress conditions and the tolerance indices for each site year were calculated.

Genotyping of mapping population

The interspecific population was genotyped using a 60 K Axiom® SNP array (61,335 SNPs) at Eurofins, WI, United States. Individual plants from each line were grown under controlled conditions in phytotron chambers. Young leaf tissue was harvested from 2 to 3 weeks old plants. These leaf tissues were freeze dried using Labconco FreeZone 6-L Console Freeze Dry System and sent to Eurofins BioDiagnostics, WI, United States, for DNA extraction and genotyping.

Population structure and linkage disequilibrium analysis

The program fastSTRUCTURE (Raj et al., 2014) was used to estimate the most likely number of clusters (K) into which the interspecific population can be grouped, and their degree of admixtures, based on genotypic data of 35,431 SNPs. The value of K was determined based on the lowest prediction error, and the smallest number of iterations for convergence. For each line, the value of Q, which is the probability of belonging to one of the clusters was derived from the matrix of contributions. A shared allele index derived from the dissimilarity matrix estimated from the SNP genotypic data were used to construct an unweighted neighbor-joining tree (Perrier et al., 2003).

Linkage disequilibrium (LD) of each chromosome was calculated as the correlation between marker-pairs among the SNP markers of each chromosome. The calculated Pearson correlation coefficient [r] was used to calculate LD decay by Quantile regression (R package “quantreg”; Koenker, 2017). The LD decay was calculated by plotting r² values as a function of genetic distance.

Association mapping

The software program GAPIT (Genome Association and Prediction Integrated Tool—R package: Lipka et al., 2012) was used to determine the association between SNP genotypes and the tolerance indices of the population. The phenotypes measured in individual environments were used for association analysis. Simultaneously, the best linear unbiased predictors (BLUP) of each trait of four environments tested were calculated using the “Ime4” package of the R3.6.1 software,² considering environmental effects as fixed and genotype as random: y ∼ (1 |rep% in % env) + (1|env) + (1 |lines) + (1 |env: lines), where rep% in % env represents replications were nested within the environments. The formula for the best linear unbiased estimate (BLUP) of phenotype is: y = Xb + Zu + e, where y, b, u, and e represent the observed phenotype, fixed effect vector, random effect vector, and residual, respectively, and X and Z represent incidence matrices. For the association analysis, the kinship coefficient matrix (K) values calculated by GAPIT by identity-by-state (IBS), were used. Multiple models of the GAPIT program including mixed linear model (MLM), multi-linear mixed model (MLMM), compressed mixed linear model (CMLM), general linear model (GLM), settlement of mixed linear models under progressively exclusive relationship (SUPER), fixed and random model circulating probability unification (FarmCPU) and factored spectrally transformed linear mixed models (FaST-LMM’; Tang et al., 2016; Kaler et al., 2020) were tested for comparison of the association analysis (Supplementary Figure 1). By comparison of the Q–Q plots of each model, which are drawn by plotting the observed and expected log₁₀ p-values, the MLMM model was chosen to report marker-trait associations.

Results

Phenotypic evaluation

Various phenotypic traits were recorded across the locations and seeding dates in 2019 and 2020 (Tables 1A,B). ANOVA showed significant differences among the interspecific lines for days to flowering, days to maturity, seed yield, flower color, plant height, seed size, and ascochyta blight severity across all locations for normal (SD1) and late seeding (SD2) treatments in 2019 and 2020. There were significant effects of the environment (year) and genotype by environment (year) interaction in all the phenotypic traits recorded for both seeding dates. The interaction of genotype and site year was not significant for plant height in SD2 (Table 1B). Some of the phenotypic traits were not recorded in all four field trials because of a lack of resources and the restrictions due to COVID-19, especially during the 2019–2020 trials at Yuma, AZ. Plot yield was among the phenotypic traits measured in all years, locations, and seeding dates. Moreover, in previous studies seed yield was found to be an important trait to measure tolerance under stress and non-stress environments. The seed yield data were used for the calculation of stress tolerance indices and association analysis. The mean value and range of population for each trait are presented in Table 2.

TABLE 1

Table 1. ANOVA for days to flowering (DTF), days to maturity (DTM), seed yield, flower color, plant height, seed size (1,000 seed weight), and ascochyta blight rating under (A) normal seeding (SD1) and (B) late seeding (SD2) for the interspecific population was evaluated under field conditions at Lucky Lake and Moose Jaw, Saskatchewan in 2019, Yuma, Arizona in 2019–2020, and Lucky Lake, Saskatchewan in 2020.

TABLE 2

Table 2. Mean and range values of the interspecific population for characters assessed under field conditions.

Correlation of tolerance indices under field conditions

Late seeding (stress condition) lowered the yield in comparison to normal seeding (non-stress condition) in all site years, except Lucky Lake in 2020 where the difference in seed yield was marginal but significant (p < 0.02; Table 2). Spearman’s rank correlation coefficients between the seed yields obtained under normal and late seeding and the stress tolerance indices for each year were calculated (Table 3). The correlation between the seed yield ranks under normal (non-stress) and late seeding (stress) conditions were positive and significant (p ≤ 0.001) for all site years, indicating a change in the ranking of the genotypes for seed yield production due to stress.

TABLE 3

Table 3. Spearmen’s rank correlation between different stress tolerance indices for each site years.

Positive significant relationships were observed between the seed yield under non-stress conditions (Yp) and TOL, MP, ATI, SSPI, and K₁STI in all site years. On the other hand, the seed yield under stress conditions (Ys) was positively correlated with MP and K₁STI, but negatively correlated with TOL, ATI, and SSPI in all site years.

Genotyping

A total of 35,429 SNPs were used to determine the population structure. Of these SNPs, 32,228 were located on the eight chromosomes of chickpea and 3,101 SNPs were on scaffolds and are not located on any of the chromosomes. After filtering for a minimum allele frequency (MAF) of 0.05, 20,679 SNPs were used for association analysis.

Linkage disequilibrium analysis

LD decay was calculated as the Pearson correlation coefficient (r²) between marker pairs of each chromosome. The LD decay differed among the eight chromosomes. The r² _max,90, which is calculated as the maximum r² in the 90th percentile of each chromosome, for chromosomes 1 to 8 is between 0.24 and 0.26. The physical distance in Mb at which LD of each chromosome has decayed to half of r² _max,90 is 0.51, 0.15, 0.25, 0.21, 0.38, 0.23, 0.24, and 0.43 Mb, for chromosomes 1 to 8, respectively (Figure 1).

FIGURE 1

Figure 1. Linkage disequilibrium (LD) decay plots of eight chromosomes of chickpea. (A) Chromosome 1, (B) Chromosome 2, (C) Chromosome 3, (D) Chromosome 4, (E) Chromosome 5, (F) Chromosome 6, (G) Chromosome 7, and (H) Chromosome 8.

Population genetic structure

To determine the population structure, the most likely number of clusters (k) was tested from 2 to 20. A k-value of 15 is best suited to describe the genetic structure of 195 chickpea interspecific lines along with the two parents CDC Consul and CDC Leader, and six CDC cultivars used as checks, CDC Orion, CDC Cory, ILC 533, ICC 4958, and ILC 3279. The probability of membership of individual lines in each cluster was estimated by admixture analysis (Figure 2; Supplementary Table 2). The grouping of the interspecific population into phylogenetic clusters based on the neighbor-joining (NJ) tree differed from the population structure analysis that the population is divided into eight clades (Figure 3).

FIGURE 2

Figure 2. The population structure of 203 chickpea interspecific lines and accessions based on k = 15.

FIGURE 3

Figure 3. Genetic relatedness among the 203 chickpea accessions, estimated by neighbor-joining method.

Association mapping

The seed yield of interspecific lines measured both in normal and heat stress treatments, and the heat stress indices calculated based on the seed yield in both the treatments, were used to identify the associated SNP markers. The BLUP values calculated based on the measurements of each trait/index in four station-years were used for association analysis. SNP markers were identified for association with seed yield both in normal and heat stress environments, and the indices calculated based on the seed yield in these two treatments are listed in Table 4. SNP marker Ca2_34600347 was identified for association with grain yield both in normal and stress (late seeding) treatments. The markers for ATI and TOL are located on chromosome 1, while the markers for other indices K₁STI and SSPI are located on chromosome 4. The Manhattan plots representing the genome-wide marker-trait associations, concerning the seven traits/indices measured, and their corresponding Q–Q plots are presented in Figure 4.

TABLE 4

Table 4. Selected marker-trait associations, identified based on BLUP values of phenotypes measured in four station-years during 2019 and 2020.

FIGURE 4

Figure 4. Manhattan plots of −log10 p-values and the corresponding quantile–quantile (Q–Q) plots of the association analysis for (A) seed yield under non-stress and (B) seed yield under stress conditions and yield indices (C) ATI, (D) K₁STI, (E) MP, (F) SSPI, and (G) TOL using a mixed linear model for four locations: Lucky Lake and Moose Jaw, Saskatchewan, 2019, Yuma, Arizona, 2019–2020, and Lucky Lake, Saskatchewan, 2020. For Manhattan plots: y-axis, −log10 p-values; x-axis, chromosome numbers. For Q–Q plots: y-axis, observed −log10 p-values; x-axis, expected −log10 p-values.

Discussion

The wild relatives of chickpea have been an invaluable source for improving elite chickpea germplasm through resistance to biotic stress (Sharma et al., 2006; Singh et al., 2008). Moreover, crop wild relatives were reported to be an important resource for genetic improvement for various abiotic stresses (Singh et al., 2008; Coyne et al., 2020; Bakir et al., 2021). Collection of various wild chickpeas, mainly C. reticulatum was reported to have biotic and abiotic stress tolerance (Singh et al., 2008; von Wettberg et al., 2018; Bakir et al., 2021). C. reticulatum is in the primary gene pool of chickpea (Bakir et al., 2021). It is known for its crossing-compatibility with domesticated chickpea (C. arietinum) and generally produces fertile progeny because of good chromosome pairing (Bakir et al., 2021). Therefore, in the current study, we have used C. reticulatum to introgress stress tolerance in cultivated kabuli and desi chickpea (C. arietinum) germplasm. In this study, we developed progeny from interspecific crosses between adapted elite cultivars, CDC Leader and CDC Consul, and 20 C. reticulatum accessions. Further lines with improved response to ascochyta blight and acceptable seed quality were selected and used in this study. Interestingly, all these selected lines were derived from CDC Leader (kabuli chickpea) and 17 wild accessions. The design of our population is the same as the nested association mapping (NAM) design. Generally, the bi-parental population has a lack of mapping precision and low genetic diversity. Therefore, using NAM helps to capture additional genetic diversity and increase genetic recombination.

Chickpea is a dry season food legume that is mostly grown on residual moisture after the rainy season. Toward the end of the growing season chickpeas often experience terminal drought stress. Moreover, if the sowing is delayed, the crop may deal with heat stress during the reproductive phase (Gaur et al., 2019; Maphosa et al., 2020; Rani et al., 2020). Low moisture and high temperatures during the flowering and early pod filling stage can substantially influence chickpea yield by forcing early maturity resulting in low biomass and a low number of pods and seeds per plant (Maphosa et al., 2020; Rani et al., 2020). Therefore, to understand the performance of interspecific populations under a suboptimal environment, the population in this study was planted at two different seeding dates (normal and late seeding). The main objective of late seeding was to expose the plants to stress conditions such as higher temperatures and low moisture during flowering. Previous studies have also used the same approach to understand the effect of stress on crop development and yield in canary seeds (Miller, 1999). Furthermore, it was observed that late seeding has a negative effect on chickpea yield (Machado et al., 2006; McKenzie et al., 2006).

The current study evaluated various phenotypic traits such as days to flowering, days to maturity, seed yield, flower color, plant height, seed size, and ascochyta blight disease resistance on both the seeding dates. There were significant effects of environment and genotype by environment interaction in all the phenotypic traits documented for both seeding dates which clearly shows that a suboptimal environment negatively impacts crop development in chickpea. Similarly, other studies have shown the negative effect of stress on various phenotypic traits grown under stress and non-stress environments (Arif et al., 2021; Jeffrey et al., 2021; Jha et al., 2021). Yield is a crucial trait and an important indicator to define tolerance between stress and non-stress conditions and has been used to describe the performance of any genotype while screening in various environments (Kaloki et al., 2019; Pour-Aboughadareh et al., 2019). Moreover, in the current study seed yield was among the phenotypic traits measured in all years, locations, and seeding dates. Therefore, we have used seed yield to calculate various stress tolerance indices. Furthermore, crop yield under stress and non-stress conditions and stress tolerance indices were used for association analysis. Seed yield was significantly reduced in the late seeded population in all the locations as compared to the normal seeded population except Lucky Lake, 2020 which could be because of less stressed conditions while the late seeded population in Lucky Lake. Similar to our findings, previous studies have also shown reduced yield under late seeded populations (Machado et al., 2006; McKenzie et al., 2006). The calculation of spearman’s rank correlation coefficients between seed yields acquired under normal and late seeding and the stress tolerance indices for each year was done. We found that the correlation between seed yield under stress and non-stress conditions was positive and significant which shows that the yield production was different in all the genotypes under non-stress conditions. Similar results were also found in previous studies (Moosavi et al., 2008; Farshadfar and Geravandi, 2013). Moreover, seed yield under non-stress conditions was positively correlated with TOL, MP, ATI, SSPI, and K₁STI which shows that selection of the genotypes based on these indices will improve the yield under non-stress conditions. Some studies have also found a positive correlation between non-stress conditions and MP, K₁STI (Farshadfar and Geravandi, 2013) and TOL, MP, SSPI, ATI (Moosavi et al., 2008). On the other hand, the seed yield under stress conditions was positively correlated with MP and K₁STI but negatively correlated with TOL, ATI, and SSPI indicating that the selection based on a higher value of MP and K₁STI will improve seed yield while selection based on TOL, ATI, and SSPI will reduce seed yield under stress conditions. Therefore, correlation analysis between seed yield and stress tolerance indices can be a good criterion for screening the best genotypes and indices used.

Linkage drag is one of the major bottlenecks to using wild species in the genetic improvement of crop plants including chickpea. Often, a large portion of entire chromosomes is affected by the linkage drag, which is a hindrance to the introgression of desirable alleles. In the current study, the LD decay observed in each chromosome is from 0.15 to 0.51 Mb, indicating that the interspecific population used in the current study is suitable for the objectives of this study. Genome-wide association study (GWAS) was used to understand the genetic basis of complex traits in chickpea (Kohli et al., 2020) and heat stress in other legumes (Tafesse et al., 2020, 2021). Since we intend to identify the heat stress loci in wild chickpea, and as well to introgess the same in an elite chickpea genetic background, we developed an inbred line population derived from crosses of wild accessions and elite chickpea cultivars. The 195 inbred line population used for association mapping of heat stress tolerance was an interspecific population derived from crosses of CDC Leader with 17 C. reticulatum accessions.

In the current study, we used the measurement of seed yield in heat stress conditions as an important criterion to identify heat stress tolerance in chickpea interspecific lines. Understanding the genetic basis of yield in heat stress conditions is important to develop heat stress tolerant high yielding chickpea varieties. Previous studies have identified quantitative trait loci (QTLs) associated with yield (Barmukh et al., 2021), outlining component traits and heat tolerance (Thudi et al., 2014) in chickpea. Barmukh et al. (2021) identified yield-associated QTLs on LG1 and LG4 using a bi-parental mapping population of ICC 4958 × DCP 92–3. In the current study, we identified markers associated with seed yield and heat tolerance indices mainly on LG1 and LG4, while markers associated with seed yield in heat stress were distributed on multiple linkage groups including LG4. Yield is a complex trait and is determined by several components such as pod weight, pod number, and the number of reproductive nodes, etc. Tafesse et al. (2020, 2021) identified the correlation between multiple components of heat stress and SNP markers associated with these components in a GWAS study on peas. We have used the measurement of seed yield in heat stress conditions as a key indicator of heat tolerance, and the indices calculated based on seed yield in normal and heat stress conditions were among the other indicators used. The SNP marker Ca2_34600347, which was identified for its association with seed yield both in normal and heat stress environments, could be a valuable marker for the marker-assisted selection of heat tolerance in chickpea.

Conclusion

Wild relatives are a source of novel alleles to improve chickpea adaptation to suboptimal environments. In the current study, we used an interspecific population derived from C. reticulatum accessions. The late seeding and exposure of this population to suboptimal growth conditions has significantly reduced the seed yield. The study characterized the variation of the interspecific population for its yield performance and yield-related indices under suboptimal conditions. Individual lines that were identified for superior performance in suboptimal conditions can be used as trait donors in breeding for stress tolerance. The trait-associated SNP markers identified can be used for marker-assisted selection in the breeding pipeline. Future detailed analysis of these SNPs would allow us to identify the genes underlying tolerance to abiotic stress, providing new alleles and molecular markers to use in chickpea crop improvement.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

SK conducted all the experiments, analyzed phenotypic data, and wrote the manuscript. KG analyzed genotypic data, conducted association analysis, and wrote the manuscript. BT conceived, designed, and supervised the research and reviewed the manuscript. All authors contributed to the article and approved the submitted version.

Funding

We acknowledge the financial support from the Saskatchewan Ministry of Agriculture through the Agricultural Development Fund.

Acknowledgments

We acknowledge Donna Lindsay for technical assistance in developing the interspecific population, Carmen Breitkreutz for technical assistance with field work, and Brent Barlow and Alison Sackville for technical assistance in the preparation of field experiments.

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.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.933277/full#supplementary-material

SUPPLEMENTARY FIGURE 1 | Q–Q plots of the association analysis for seed yield under non-stress (Yp), seed yield under stress conditions (Ys) and yield indices ATI, K1STI, MP, SSPI, and TOL from various GWAS models (CMLM, FarmCPU, FaST-LMM, GLM, MLM, MLMM, SUPER).

SUPPLEMENTARY TABLE 1 | Daily temperature and precipitation recorded at Lucky Lake and Moose Jaw, Saskatchewan, Canada in 2019, Yuma, Arizona, USA in 2019–2020 and Lucky Lake, Saskatchewan, Canada in 2020 during growing season.

Footnotes

1. ^https://www.canada.ca/en/environment-climate-change/news/2021/06/environment-and-climate-change-canada-warns-western-canadians-about-dangerous-record-high-temperatures.html

2. ^www.r-project.org

References

Arif, A., Parveen, N., Waheed, M. Q., Atif, R. M., Waqar, I., and Shah, T. M. (2021). A comparative study for assessing the drought-tolerance of chickpea under varying natural growth environments. Front. Plant Sci. 11:607869. doi: 10.3389/fpls.2020.607869

PubMed Abstract | CrossRef Full Text | Google Scholar

Bakir, M., Sari, D., Sari, H., Waqas, M., and Atif, R. M. (2021). “Chickpea wild relatives: potential hidden source for the development of climate resilient chickpea varieties,” in Wild Germplasm for Genetic Improvement in Crop Plants. eds. M. T. Azhar and S. H. Wani (Cambridge, MA, USA: Academic Press), 269–297.

Google Scholar

Barmukh, R., Soren, K. R., Madugula, P., Gangwar, P., Shanmugavadivel, P. S., Bharadwaj, C., et al. (2021). Construction of a high-density genetic map and QTL analysis for yield, yield components and agronomic traits in chickpea (Cicer arietinum L.). PLoS One 16:e0251669. doi: 10.1371/journal.pone.0251669

PubMed Abstract | CrossRef Full Text | Google Scholar

Bohra, A., Kilian, B., Sivasankar, S., Caccamo, M., Mba, C., McCouch, S. R., et al. (2021). Reap the crop wild relatives for breeding future crops. Trends Biotechnol. 40, 412–431. doi: 10.1016/j.tibtech.2021.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Buckler, E. S., Holland, J. B., Bradbury, P. J., Acharya, C. B., Brown, P. J., Browne, C., et al. (2009). The genetic architecture of maize flowering time. Science 325, 714–718. doi: 10.1126/science.1174276

CrossRef Full Text | Google Scholar

Chongo, G., Gossen, B., Buchwaldt, L., Adhikari, T., and Rimmer, S. (2004). Genetic diversity of Ascochyta rabiei in Canada. Plant Dis. 88, 4–10. doi: 10.1094/PDIS.2004.88.1.4

CrossRef Full Text | Google Scholar

Coyne, C. J., Agrawal, S. K., von Wettberg, E. J. B., Marques, E., Berger, J., Redden, R. J., et al. (2020). Potential and limits of exploitation of crop wild relatives for pea, lentil, and chickpea improvement. Legum. Sci. 2:e36. doi: 10.1002/leg3.36

CrossRef Full Text | Google Scholar

Devasirvatham, V., Gaur, P. M., Raju, T. N., Trethowan, R. M., and Tan, D. K. Y. (2015). Field response of chickpea (Cicer arietinum L.) to high temperature. Field Crop Res. 172, 59–71. doi: 10.1016/j.fcr.2014.11.017

CrossRef Full Text | Google Scholar

Devasirvatham, V., and Tan, D. K. Y. (2018). Impact of high temperature and drought stresses on chickpea production. Agronomy 8:145. doi: 10.3390/agronomy8080145

CrossRef Full Text | Google Scholar

Dutra, W. F., Guerra, Y. L., Ramos, J. P. C., Fernandes, P. D., Silva, C. R. C., Bertioli, D. J., et al. (2018). Introgression of wild alleles into the tetraploid peanut crop to improve water use efficiency, earliness and yield. PLoS One 13:e0198776. doi: 10.1371/journal.pone.0198776

PubMed Abstract | CrossRef Full Text | Google Scholar

FAO (2019) Chickpea production data. Available at: http://www.fao.org/faostat/en/#data/QCL/visualize

Google Scholar

Farshadfar, E., and Geravandi, M. (2013). Repeatability of drought tolerance indices in chickpea genotypes over stress and non-stress environments. Acta Agron. Hung. 61, 123–137. doi: 10.1556/AAgr.61.2013.2.4

CrossRef Full Text | Google Scholar

Farshadfar, E., and Sutka, J. (2003). Multivariate analysis of drought tolerance in wheat substitution lines. Cereal Res. Commun. 31, 33–40. doi: 10.1007/BF03543247

CrossRef Full Text | Google Scholar

Gangurde, S. S., Wang, H., Yaduru, S., Pandey, M. K., Fountain, J. C., Chu, Y., et al. (2020). Nested-association mapping (NAM)-based genetic dissection uncovers candidate genes for seed and pod weights in peanut (Arachis hypogaea). Plant Biotechnol. J. 18, 1457–1471. doi: 10.1111/pbi.13311

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaur, P. M., Krishnamurthy, L., and Kashiwagi, J. (2008). Improving drought-avoidance root traits in chickpea (Cicer arietinum L.), current status of research at ICRISAT. Plant Prod. Sci. 11, 3–11. doi: 10.1626/pps.11.3

CrossRef Full Text | Google Scholar

Gaur, P. M., Samineni, S., Krishnamurthy, L., Kumar, S., Ghanem, M. E., Beebe, S., et al. (2015). High temperature tolerance in grain legumes. Legume Perspect. 7, 23–24.

Google Scholar

Gaur, P. M., Samineni, S., Thudi, M., Tripathi, S., Sajja, S. B., Jayalakshmi, V., et al. (2019). Integrated breeding approaches for improving drought and heat adaptation in chickpea (Cicer arietinum L.). Plant Breed. 138, 389–400. doi: 10.1111/pbr.12641

CrossRef Full Text | Google Scholar

Gorim, L. Y., and Vandenberg, A. (2017). Root traits, nodulation and root distribution in soil for five wild lentil species and Lens culinaris (medik.) grown under well-watered conditions. Front. Plant Sci. 8:1632. doi: 10.3389/fpls.2017.01632

PubMed Abstract | CrossRef Full Text | Google Scholar

Hajjar, R., and Hodgkin, T. (2007). The use of wild relatives in crop improvement: a survey of developments over the last 20 years. Euphytica 156, 1–13. doi: 10.1007/s10681-007-9363-0

CrossRef Full Text | Google Scholar

Hussain, M. M., Rauf, S., and Warburton, M. L. (2019). Development of drought-tolerant breeding lines derived from Helianthus annuus × H. argophyllus interspecific crosses. Plant Breed. 138, 862–870. doi: 10.1111/pbr.12731

CrossRef Full Text | Google Scholar

Jeffrey, C., Trethowan, R., and Kaiser, B. (2021). Chickpea tolerance to temperature stress: status and opportunity for improvement. J. Plant Physiol. 267:153555. doi: 10.1016/j.jplph.2021.153555

PubMed Abstract | CrossRef Full Text | Google Scholar

Jha, U. C., Chaturvedi, S. K., Bohra, A., Basu, P. S., Khan, M. S., and Barh, D. (2014). Abiotic stresses, constraints and improvement strategies in chickpea. Plant Breed. 133, 163–178. doi: 10.1111/pbr.12150

CrossRef Full Text | Google Scholar

Jha, U. C., Nayyar, H., Palakurthi, R., Jha, R., Valluri, V., Bajaj, P., et al. (2021). Major QTLs and potential candidate genes for heat stress tolerance identified in chickpea (Cicer arietinum L.). Front. Plant Sci. 12:655103. doi: 10.3389/fpls.2021.655103

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaler, A. S., Gillman, J. D., Beissinger, T., and Purcell, L. C. (2020). Comparing different statistical models and multiple testing corrections for association mapping in soybean and maize. Front. Plant Sci. 10:1794. doi: 10.3389/fpls.2019.01794

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaloki, P., Devasirvatham, V., and Tan, D. K. (2019). “Chickpea abiotic stresses: combating drought, heat and cold,” in: Abiotic and Biotic Stress in Plants ed. A. B. Oliveirade (London, UK: IntechOpen).

Google Scholar

Koenker, R. (2017). Quantile regression: 40 years on. Annu. Rev. Econ. 9, 155–176. doi: 10.1146/annurev-economics-063016-103651

CrossRef Full Text | Google Scholar

Kohli, P. S., Verma, P. K., Verma, R., Parida, S. K., Thakur, J. K., and Giri, J. (2020). Genome-wide association study for phosphate deficiency responsive root hair elongation in chickpea. Funct. Integr. Genomics 20, 775–786. doi: 10.1007/s10142-020-00749-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Korte, A., and Farlow, A. (2013). The advantages and limitations of trait analysis with GWAS: a review. Plant Methods 9:29. doi: 10.1186/1746-4811-9-29

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakmes, A., Jhar, A., Penmetsa, R. V., Wei, W., Brennan, A. C., and Kahriman, A. (2022). The quantitative genetics of flowering traits in wide crosses of chickpea. Agriculture 12:486. doi: 10.3390/agriculture12040486

CrossRef Full Text | Google Scholar

Lipka, A. E., Tian, F., Wang, Q., Peiffer, J., Li, M., Bradbury, P. J., et al. (2012). GAPIT: genome association and prediction integrated tool. Bioinformatics 28, 2397–2399. doi: 10.1093/bioinformatics/bts444

PubMed Abstract | CrossRef Full Text | Google Scholar

Machado, S., Humphreys, C., Tuck, B., and Corp, M. (2006). Seeding date, plant density, and cultivar effects on chickpea yield and seed size in eastern Oregon. Crop Manag. 5, 1–7. doi: 10.1094/CM-2006-0621-01-RS

CrossRef Full Text | Google Scholar

Maphosa, L., Richards, M. F., Norton, S. L., and Nguyen, G. N. (2020). Breeding for abiotic stress adaptation in chickpea (Cicer arietinum L.): A comprehensive review. Crop Breed. Genet. Genom. 2:e200015. doi: 10.20900/cbgg20200015

CrossRef Full Text | Google Scholar

McKenzie, R. H., Middleton, A. B., and Bremer, E. (2006). Effect of seeding date and rate on desi chickpea in southern Alberta. Can. J. Plant Sci. 86, 717–721. doi: 10.4141/P05-123

CrossRef Full Text | Google Scholar

Merchuk-Ovnat, L., Fahima, T., Ephrath, J. E., Krugman, T., and Saranga, Y. (2017). Ancestral QTL alleles from wild emmer wheat enhance root development under drought in modern wheat. Front. Plant Sci. 8:703. doi: 10.3389/fpls.2017.00703

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, P. R. (1999). Effect of varying seeding date on crop development, yield and yield components in canaryseed. Can. J. Plant Sci. 80, 83–86. doi: 10.4141/P99-047

CrossRef Full Text | Google Scholar

Moosavi, S. S., Yazdi-Samadi, B., Naghavi, M. R., Zali, A. A., Dashti, H., and Pourshahbazi, A. (2008). Introduction of new indices to identify relative drought tolerance and resistance in wheat genotypes. Desert 12, 165–178. doi: 10.22059/JDESERT.2008.27115

CrossRef Full Text | Google Scholar

Perrier, X., Flori, A., and Bonnot, F. (2003). “Data analysis methods,” in Genetic Diversity of Cultivated Tropical Plants. eds. P. Hamon, M. Seguin, X. Perrier, and J. C. Glaszmann (Enfield, USA: Science publishers), 43–76.

Google Scholar

Placido, D. F., Campbell, M. T., Folsom, J. J., Cui, X., Kruger, G. R., Baenziger, P. S., et al. (2013). Introgression of novel traits from a wild wheat relative improves drought adaptation in wheat. Plant Physiol. 161, 1806–1819. doi: 10.1104/pp.113.214262

PubMed Abstract | CrossRef Full Text | Google Scholar

Pour-Aboughadareh, A., Yousefian, M., Moradkhani, H., Vahed, M. M., Poczai, P., and Siddique, K. H. M. (2019). iPASTIC: an online toolkit to estimate plant abiotic stress indices. Appl. Plant Sci. 7:e11278. doi: 10.1002/aps3.11278

PubMed Abstract | CrossRef Full Text | Google Scholar

Raj, A., Stephens, M., and Pritchard, J. K. (2014). fastSTRUCTURE: variational inference of population structure in large SNP data sets. Genetics 197, 573–589. doi: 10.1534/genetics.114.164350

PubMed Abstract | CrossRef Full Text | Google Scholar

Rani, A., Devi, P., Jha, U. C., Sharma, K. D., Siddique, K. H., and Nayyar, H. (2020). Developing climate-resilient chickpea involving physiological and molecular approaches with a focus on temperature and drought stresses. Front. Plant Sci. 10:1759. doi: 10.3389/fpls.2019.01759

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosielle, A. A., and Hamblin, J. (1981). Theoretical aspects of selection for yield in stress and nonstress environments. Crop. Sci. 21, 943–946. doi: 10.2135/cropsci1981.0011183X002100060033x

CrossRef Full Text | Google Scholar

SAS Institute Inc. (1999). SAS/STAT® User’s Guide, Version 8, SAS Institute Inc., Cary, NC.

Google Scholar

Sharma, H. C., Bhagwat, M. P., Pampapathy, G., Sharma, J. P., and Ridsdill-Smith, T. J. (2006). Perennial wild relatives of chickpea as potential sources of resistance to Helicoverpa armigera. Genet. Resour. Crop. Evol. 53, 131–138. doi: 10.1007/s10722-004-1951-4

CrossRef Full Text | Google Scholar

Singh, M., Ceccarelli, S., and Hamblin, J. (1993). Estimation of heritability from varietal trials data. Theor. Appl. Genet. 86, 437–441. doi: 10.1007/BF00838558

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, R., Sharma, P., Varshney, R. K., Sharma, S. K., and Singh, N. K. (2008). Chickpea improvement: role of wild species and genetic markers. Biotechnol. Genet. Eng. Rev. 25, 267–314. doi: 10.5661/bger-25-267

PubMed Abstract | CrossRef Full Text | Google Scholar

Tafesse, E. G., Gali, K. K., Lachagari, V. B. R., Bueckert, R., and Warkentin, T. D. (2020). Genome-wide association mapping for heat stress responsive traits in field pea. Int. J. Mol. Sci. 21:2043. doi: 10.3390/ijms21062043

PubMed Abstract | CrossRef Full Text | Google Scholar

Tafesse, E. G., Gali, K. K., Lachagari, V. B. R., Bueckert, R., and Warkentin, T. D. (2021). Genome-wide association mapping for heat and drought adaptive traits in pea. Genes 12:1897. doi: 10.3390/genes12121897

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, Y., Liu, X., Wang, J., Li, M., Wang, Q., Tian, F., et al. (2016). GAPIT version 2: an enhanced integrated tool for genomic association and prediction. Plant Genome 9. doi: 10.3835/plantgenome2015.11.0120

PubMed Abstract | CrossRef Full Text | Google Scholar

Thudi, M., Upadhyaya, H. D., Rathore, A., Gaur, P. M., Krishnamurthy, L., Roorkiwal, M., et al. (2014). Genetic dissection of drought and heat tolerance in chickpea through genome-wide and candidate gene-based association mapping approaches. PLoS One 9:e96758. doi: 10.1371/journal.pone.0096758

PubMed Abstract | CrossRef Full Text | Google Scholar

Varshney, R. K., Song, C., Saxena, R. K., Azam, S., Yu, S., Sharpe, A. G., et al. (2013). Draft genome sequence of chickpea (Cicer arietinum) provides a resource for trait improvement. Nat. Biotechnol. 31, 240–246. doi: 10.1038/nbt.2491

PubMed Abstract | CrossRef Full Text | Google Scholar

Varshney, R. K., Thudi, M., Roorkiwal, M., He, W., Upadhyaya, H. D., Yang, W., et al. (2019). Resequencing of 429 chickpea accessions from 45 countries provides insights into genome diversity, domestication and agronomic traits. Nat. Genet. 51, 857–864. doi: 10.1038/s41588-019-0401-3

PubMed Abstract | CrossRef Full Text | Google Scholar

von Wettberg, E. J. B., Chang, P. L., Başdemir, F., Carrasquila-Garcia, N., Korbu, L. B., Moenga, S. M., et al. (2018). Ecology and genomics of an important crop wild relative as a prelude to agricultural innovation. Nat. Commun. 9:649. doi: 10.1038/s41467-018-02867-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Holland, J. B., McMullen, M. D., and Buckler, E. S. (2008). Genetic design and statistical power of nested association mapping in maize. Genetics 178, 539–551. doi: 10.1534/genetics.107.074245

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Yasmin, F., and Song, B.-H. (2019). Neglected treasures in the wild — legume wild relatives in food security and human health. Curr. Opin. Plant Biol. 49, 17–26. doi: 10.1016/j.pbi.2019.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar