Screening of heat-tolerant Ethiopian chickpea accessions: Assessment of phenological and agromorphological traits and genomic relationships

Abstract

A major agronomic challenge for chickpea (Cicer arietinum L.) production is temperatures above 35 °C, which causes reduced fertility and seed development. This study was aimed at assessing the phenotypic variation of chickpea genotypes under variable heat stress conditions. Chickpea genotypes were grown in heat-stressed locations in both Ethiopia and India to assess phenotypic variation for heat tolerance. In addition, genomic relationships among the genotypes were assessed using genome-wide single nucleotide polymorphism (SNP) markers. A total of 121 genotypes were assessed at three field sites in Ethiopia, under heat stress and nonstress conditions, and 57 genotypes were assessed under high-heat-stress conditions at the International Crops Research Institute for the Semiarid Tropics (ICRISAT) in Hyderabad, India. Data for five phenological and seven agromorphological traits were recorded. The results showed that the chickpea genotypes were severely affected by excessive heat at Hyderabad as compared with those planted under non-heat-stress conditions in Debre Zeit, Ethiopia. At extremely high temperatures, chickpea plants exhibited reduced floral initiation, arrested seed and pod development, shortened life cycles, and reduced plant height, seed size, grain yield, and yield-related traits. Across stressed and nonstressed environments, there were highly significant differences among the genotypes for most of the traits (ANOVA, P ≤ .001). Under heat-stressed environments, DZ-Cr-0034 was found to be a highly tolerant, whereas DZ-Cr-0026 was found to be a highly sensitive genotype. Genetic relationships among the genotypes were determined using 5,722 SNPs, revealing a single group of Ethiopian genotypes with small number of cultivars showing introgression from Middle Eastern germplasm. This study clearly demonstrated that there is genetic variability in chickpea for heat tolerance that can be harnessed to meet expected shifts towards warmer climatic conditions.

Abbreviations
• D50F
• days to 50% flowering
• DFF
• days to first flowering
• DFP
• days to first podding
• DPM
• days to physiological maturity
• GY
• grain yield
• ICRISAT
• International Crops Research Institute for the Semiarid Tropics
• PB
• plant biomass
• PCA
• principal component analysis
• PCoA
• principal coordinates analysis
• PHI
• percentage of harvest index
• PHT
• plant height
• SNP
• single nucleotide polymorphism
• TPP
• total number of pods per plant
• TSP
• total number of seeds per plant.

INTRODUCTION

In sub-Saharan Africa, chickpea (Cicer arietinum L.) is predominately cultivated in Ethiopia, Sudan, Eritrea, Kenya, Tanzania and Malawi (Bejiga & Van der Maesen, 2006). Chickpea was likely introduced into sub-Saharan Africa through Ethiopia, and indeed Vavilov designated Ethiopia as a secondary center of diversity for chickpea (Sokolkova et al., 2020). Ethiopian production accounts for 90% (Verkaar et al., 2017) of the total production in sub-Saharan Africa and 3% of the world's production (FAO, 2018) which makes Ethiopia sixth among the world's leading chickpea producers. Chickpea is a source of dietary protein, a legume that improves soil fertility, and an important cash crop for Ethiopia's small farm holders, while also generating foreign currency through export markets (Bejiga & Eshete, 1996).

High temperatures (30–35 °C) limit growth and yield, with primary impact on flower development, fertility, and pod development, and overall impact on plant biomass (PB) and N₂ fixation of chickpea (Devasirvatham, 2012; Krishnamurthy et al., 2011; Summerfield et al., 1984; Summerfield et al., 1990). Devasirvatham (2012) reported that high temperatures affect several chickpea physiological and agronomical traits. Several reports describe reduced yields when altered planting times expose the crop to high temperatures, although this experimental approach confounds several different factors such as day length and temperature. Kalra et al. (2008) reported a yield reduction from 53 to 330 kg ha⁻¹ among early-maturing genotypes, and Upadhyaya et al. (2011) calculated a 10–15% grain yield (GY) loss for each 1 °C rise above optimum temperatures in India. Only two heat-tolerant genotypes were identified in the international mini-core collection of 300 genotypes that represents the global diversity of chickpea (Upadhyaya et al., 2011). Similarly, Devasirvatham (2012) reported that among 167 chickpea genotypes, only two were heat tolerant at ≥38 °C. In Ethiopia, although 22 improved chickpea varieties were released for drought tolerance (Fikre, 2014), no variety has been released yet for heat stress tolerance. Chickpeas are increasingly grown in Ethiopia's lowland regions, including Melka Werer, Humera, Kobo, and South Omo, which exposes the crop to higher temperatures (Bitew & Asargew, 2015).

Understanding the genetic basis of heat tolerance traits is needed in chickpea, both for genomic-guided crop improvement and to identify underlying genes (Devasirvatham, 2012; Paul et al., 2018). As a prelude to genetic analyses, one must first identify genotypes that differ in heat tolerance traits. Genome-wide association studies have the potential to identify causal genetic intervals but require large and properly structured study populations, detailed and accurate phenotyping, and significant resources for genomic analysis. Alternatively, one can understand the degree of genetic variation in tested materials so that suitable parents can be selected to develop biparental populations for trait–locus mapping and molecular breeding (Basu et al., 2019; Kujur et al., 2015).

There are few reports on heat tolerance in chickpea, and only a single report considers heat tolerance in Ethiopia using 18 early-maturing Ethiopian and Indian genotypes (Molla et al., 2018). Although Molla et al. (2018) demonstrated the feasibility of discovering heat-tolerant genotypes, there remains an urgent need to identify additional sources for heat tolerance in breeding material. The objectives of the present study were to assess the phenotypic variation in a multilocation field trial for heat tolerance and to determine genomic relationships among the chickpea genotypes.

MATERIALS AND METHODS Plant materials

A total of 121 and 57 chickpea genotypes were used at three experimental sites in Ethiopia (Kobo, Werer, and Debre Zeit) and at one field site in India (International Crops Research Institute for the Semiarid Tropics [ICRISAT], Hyderabad), respectively (Supplemental Table S1). The chickpea genotypes used in Ethiopia were 81 landraces including a local check, 34 improved lines, and six commercial cultivars. Among 81 landraces, 76 originated from low and mid-altitudes (<1,900 m asl) and five genotypes were from high altitudes (>2,400 m asl). The 57 genotypes used in India (Hyderabad) include 10 landraces, 40 improved lines and commercial cultivars, and seven standard heat-sensitive and -tolerant checks (Supplemental Table S1). Plants were gown at three experimental locations in Ethiopia (Kobo [2017], Werer [2018], and Debre Zeit [2016]) and at one location in India (Hyderabad, ICRISAT [2017]).

Debre Zeit is located at 1,900 m asl (08°44′ N, 38°58′ E). The soil type is black (Vertisols) having pH value ranging from 6.0 to 7.5. Kobo is located at 1,490 m asl (12°09′67′′ N, 39°37′48′′ E). It is dominated by Eutric Vertisol and Eutric Fuvisol soils with pH values range from 7.4 to 8.5. Werer is located at 750 m asl (9°16′ N, 40°9′ E). The soil types are Salic Fluvisols, Eutric Fluvisols, and Eutric Vertisols. The pH of the soils ranges from 7.4 to 8.4. Hyderabad, India, is located at 545 m asl (17.53° N, 78.27° E). The soil type is black (Vertisol). The pH of the soil is alkaline and ranges from 8.35 to 8.85.

Experimental design and field management

An α-lattice design with two replications was used at Hyderabad, whereas three replications were used at three locations in Ethiopia. Plants were grown from March to June 2016 at Debre Zeit, March to June and February to May 2017 at Hyderabad and Kobo, respectively, and February to May 2018 at Werer. Individual plots with two rows were 1 m by 0.6 m and contained 20 seeds of a single genotype planted at 0.1-, 0.3-, and 1-m spacing between plants, rows, and replications, respectively. Due to a shortage of space, only a single row with 10 seeds was planted at Hyderabad, and spacing between plants and replications were the same as above.

In all locations, plants were grown using supplemental irrigation to avoid drought stress. Plots were managed by frequent manual weeding. Chemical sprayings were applied for pod borer (Helicoverpa armigera), Ascochyta blight (Ascochyta rabiei), and Fusarium wilt (Fusarium oxysporum). Root rot (Rhizoctonia solani) was the most frequent disease in Hyderabad, and we controlled it by soil solarization and regular removal of infected individuals.

Core Ideas

• Heat-tolerant chickpea genotypes were screened under three heat-stressed environments.
• Chickpea genotypes were evaluated under non-heat-stress environment for comparison.
• Genetic variability of genotypes for phenological and agromorphological traits were determined.
• The correlations among traits and genotypes were demonstrated.
• The genomic relations among the genotypes were assessed using SNPs.

Temperatures of the experimental field sites

Debre Zeit was used as a control (non-heat-stressed environment), whereas the remaining field sites were heat-stressed environments. The monthly average minimum and maximum temperatures at Debre Zeit varied from 13.3 to 13.8 °C and 26.8 to 31.4 °C, respectively. Hyderabad had the most extreme monthly average minimum and maximum temperatures, which ranged from 19 to 26.5 °C and 35.7 to 39.7 °C, respectively. The monthly average minimum and maximum temperatures at Kobo were 12.4–17 °C and 28.9–33 °C, respectively (Supplemental Figure S1a). The monthly average minimum and maximum temperatures at Werer were 18.1–23.3 °C and 33–37.1 °C, respectively (Supplemental Figure S1b). Maximum daily temperatures at Kobo, Werer, and India (Hyderabad) were above 36, 40, and 45 °C, respectively. At Hyderabad, the length of the day (14 h) was longer than the night (10 h) during the field experiment. All Ethiopian field sites (Kobo, Debre Zeit, and Werer) had similar day and night lengths of 13 and 11 h, respectively.

Phenotyping

Twelve phenological and agronomic traits were recorded, namely, days to first flowering (DFF), days to 50% flowering (D50F), days to first podding (DFP), grain-filling period (GFP), days to physiological maturity (DPM), plant height (PHT), GY, 100 seeds weight (HSW), total number of pods per plant (TPP), total number of seeds per plant (TSP), PB, and percentage of harvest index (PHI). Two additional traits, heat visual scores (1–5, tolerant to sensitive) and percentage of pod set, were also scored at the Werer field site only. Five individual plants were randomly selected for most traits, and values for each trait were calculated as an average. Data for GY, PB, and PHI were taken on a per-plot basis.

DNA extraction

Leaf samples from 15-d-old seedlings of 54 genotypes of chickpea were collected from the greenhouse of University of California-Davis. Seventy- to one hundred-milligram leaf samples were placed in a 12- × 8-well strip tube with strip cap (Marsh Biomarket) together with two 4-mm stainless steel grinding balls (Spex CertiPrep) for DNA extraction using a Retsch mixer mill and Qiagen 96-well PlantDNeasy kits. DNA quality and quantity were assessed using a NanoDrop 1000 spectrophotometer.

Genotyping-by-sequencing library development and sequencing

Individual genome samples were prepared for genotyping using a restriction site associated DNA RAD-seq (restriction site associated DNA sequencing) approach, as described in von Wettberg et al. (2018). Briefly, DNA was digested with HindIII and NlaIII adaptors added as barcodes and to facilitate en masse polymerase chain reaction (PCR) amplification. Library quality was assessed using an Agilent 2100 bio-analyzer. Pooled samples were sequenced as 100-bp reads on an Illumina HiSeq4000 at the University of California-Davis Genomics Facility. Calling of single nucleotide polymorphisms (SNPs) was performed with the Genome Analysis Tool Kit pipeline, following Genome Analysis Tool Kit best practices (McKenna et al., 2010), with filtering thresholds described in von Wettberg et al. (2018) and Sani et al. (2018). Raw sequencing data can be found in National Center for Biotechnology Information (NCBI) Umbrella BioProject PRJNA353637, which includes Illumina data from RAD sequencing of over 500 Ethiopian landrace accessions (PRJNA507628).

Data analysis

All phenological and agronomic traits are considered as dependent variables; however, genotype, block, replication, and location (field site) were taken as independent variables (Gomez & Gomez, 1984). For all traits, raw data means were calculated. For the combined ANOVA, locations and genotypes were considered as random and fixed variables, respectively. Thus, ANOVA was computed based on a mixed linear additive model, using the general linear model procedure of the SAS 9.0 software (SAS Institute, 2002) for single and combined location analyses of phenotypic data (Gomez & Gomez, 1984). The total variation (${\mathrm{\sigma }}_{\mathrm{p}}^2$) for each trait was partitioned into variance components—genotypic (${\mathrm{\sigma }}_{\mathrm{g}}^2$), environmental (${\mathrm{\sigma }}_{\mathrm{e}}^2$), and genotype × environmental interaction (${\mathrm{\sigma }}_{\mathrm{ge}}^2$) variance—using the VARCOMP procedure of SAS. These variance components, coefficients of variation, and genetic advance (GA) per location and combined over locations were estimated based on the methods of Singh and Chaudhary (1985). Broad-sense heritability (H²) was estimated as suggested by Eckebil et al. (1977). In addition, cluster and principal component analyses were done using JMP software, version 15 (SAS Institute, 2019).

We assessed genetic distance and principal coordinates analysis (PCoA) to visulaize the relationship among genotypes. We examined percentage of polymorphism and heterozygosity (genetic diversity) to assess the extent of molecular variability in the landraces, lines, and improved varieties. All the molecular data were analyzed using GenALEX 6.5 (Peakall & Smouse, 2012).

RESULTS Genotypic response to heat stress

The mean, minimum, and maximum DPM of genotypes under heat-stressed environments was 7, 8, and 10 d shorter than in non-heat-stressed environments, respectively (Table 1). At Hyderabad, the daily average temperatures were ≥40 °C at the reproductive stage. As a consequence, two genotypes (DZ-2012-ck-0026 and ICC 4567) did not initiate flowers and pods at all. For instance, accessions 235721, 41268, 41289, Habru, and DZ-2012-ck-235 attained DFF but did not manifest DFP and D50F. Delays of DFP responses were also observed in ICCV 7102, DZ-2012-ck-231, DZ-2012-ck-208, DZ-2012-ck-229, and Minjar, possibly due to the severity of heat stress. The minimum and maximum PHT of genotypes at Kobo and Hyderabad were shorter by 12.4 and 16.2 cm and by 16.9 and 19.2 cm than genotypes at Debre Zeit, respectively. The mean, minimum, and maximum GY of genotypes at Kobo was 2.5-, 10-, and 2-fold lower than at Debre Zeit. At Hyderabad, the GY values were 39-, 15-, and 16-fold lower than at Debre Zeit (Table 1).

TABLE 1 Mean, SE, minimum, maximum, and range values of traits for two heat-stressed environments (Kobo, Hyderabad) and a non-heat-stressed (Debre Zeit) environment

	Kobo				Debre Zeit				Hyderabad
Trait	Mean ± SE	Min.	Max.	Range	Mean ± SE	Min.	Max.	Range	Mean ± SE	Min.	Max.	Range
DFF, d	42.1 ± 0.27	37.3	48	10.7	47.0 ± 0.32	40.3	53.7	13.4	41.6 ± 0.6	35	55	20
D50F, d	45.0 ± 0.28	39.7	50.7	11	51.7 ± 0.29	46.7	64.7	18	44.5 ± 0.9	39	66	27
DFP, d	47.8 ± 0.3	43	53.7	10.7	57.6 ± 0.32	50.3	68.7	18.4	37.1 ± 2.9	46.5	61.5	15
GFP, d	40.9 ± 0.34	37	46.3	9.3	41.7 ± 0.46	34.3	56.3	22	34.6 ± 1.1	24.5	49.5	25
DPM, d	85.9 ± 0.47	79.3	97	17.7	93.4 ± 4.7	88.3	109.7	21.4	85.5 ± 0.74	80.5	101.5	21.5
PHT, cm	24.3 ± 0.4	18.6	34.8	16.2	40.2 ± 0.43	31	51.7	20.7	21.9 ± 0.44	14.8	32.5	17.7
TPP	15.2 ± 0.99	5.3	31.7	26.4	32.0 ± 1.8	10	115	105	0.5 ± 0.09	0	2.1	2.1
TSP	19.07 ± 1.4	5.7	45.3	39.6	38.2 ± 2.1	13.3	129.7	116.4	0.6 ± 0.13	0	2.8	2.8
HSW, g	11.3 ± 0.35	6	26.6	20.6	13.3 ± 0.35	8.2	23.3	15.1	7.2 ± 1.0	0	21.5	21.5
PB, g	146.3 ± 9.7	27.7	446	418.3	228.3 ± 7.3	98	659	561	30.7 ± 3.0	0	87.0	87
GY, g	28.8 ± 2.6	1.5	152.3	150.8	71.0 ± 2.9	14.6	255	240.4	1.8 ± 0.46	0	15.7	15.7
PHI	19.3 ± 1.1	2.2	52.6	50.4	31.8 ± 1.0	7.8	61.9	54.1	6.6 ± 1.3	0	38.8	38.8

Note. DFF, days to first flowering; D50F, days to 50% flowering; DFP, days to first podding; GFP, grain-filling period; DPM, days to physiological maturity; PHT, plant height; TPP, total number of pods per plant; TSP, total number of seeds per plant; HSW, hundred-seed weight; PB, plant biomass; GY, grain yield; PHI, percentage of harvest index.

Similar to other heat-stressed environments, the maximum temperatures, which reached 39°C at Werer (Supplemental Figure S1b), had no effect on the germination of the genotypes. Some heat-tolerant genotypes and checks survived and had reasonable pod sets (40–80%; Supplemental Figure S2). Both HVS and percentage of pod set data revealed that only 5% of the accessions used at Werer were highly heat tolerant (data not shown).

ANOVA

An ANOVA revealed that at Hyderabad, there were highly significant (P ≤ .001) differences among genotypes for nine agronomic and phenological traits. Similar significant differences were found among genotypes for all traits at Kobo. An ANOVA revealed significant difference among sites (P ≤ .01) except for DPM and TPP. Significant genotypes × environment interactions (P ≤ .05) were found for all traits except TPP, TSP, and PB. Moreover, at Debre Zeit, ANOVA revealed that genotypes differed significantly (P ≤ .001) for all traits except PHT (Table 2). The environmental, genotypic, and phenotypic variances of genotypes at Kobo were very high for PB, GY, and other yield-related traits. The environmental variance at Hyderabad for phenological traits was relatively high. Similarly, PB had the highest genotypic and phenotypic variances. In the pooled data of stressed environments, PB had the highest environmental, genotypic, genotype × environmental interaction and phenotypic variances. All traits had higher coefficient of phenotypic variances than the corresponding coefficient of genotypic variances. At Hyderabad, coefficients of genotypic and phenotypic variances for GY and yield-related traits were extremely higher than the phenological traits (Table 3).

TABLE 2 ANOVA for 12 traits under heat-stressed, pool of stressed, and nonstressed environments

	Variable
Source	DFF	D50F	DFP	GFP	DPM	PHT	TPP	TSP	HSW	PB	GY	PHI
Hyderabad
Rep(1)	87.3***	286.4***	419.2	327.6**	99.2	65.7*	3.1***	3.7***	135.2**	3,896.8***	41.9**	135.5*
Block(Rep = 2)	20.9*	46.9	288.7	40.1	76.3	75.3	0.001*	0.01	33.4	1,713.1**	0.5	6.6
Gen(56)	33.3***	61.0***	608.0***	73.2*	114.2	31.4*	0.8***	1.6***	97.2***	788.5***	18.1***	168.5***
Kobo
Rep(2)	10.5*	121.3***	190.3***	185.8***	223.7***	576.6***	4,478.2***	7,964***	43.6***	149,115.9***	13,175.5***	872.3***
Block(Rep)(30)	2.7	5.6**	11.7***	19.7***	26.7***	27.0***	223.5***	500.0***	8.3***	12,866.9***	1,011.8***	67.9*
Gen(120)	19.1***	17.0***	17.8***	16.3***	45.4***	18.9***	121.1***	236.4***	15.5***	13,678.2***	1,098.7***	255.3***
Debre Zeit
Rep(2)	60.1**	19.2*	8.7	60.1**	11.5	427.0***	68.05***	164.5***	1.06*	561.8*	204.9**	16.3
Block(Rep)(30)	26.9***	11.2*	16.7***	26.9***	9.7*	49.2**	22.9**	26.1**	0.37	97.0	50.0	8.3
Gen(120)	45.5***	16.9***	17.1***	45.5***	37.8***	31.2*	869.5***	1,170.2***	19.4***	16,532.2***	2,390.8***	293.2***
Pool of two heat-stressed environments (Kobo and Hyderabad)
Loc(1)	161.7**	1,294.6***	3,154.3***	930.2***	26.6	521.6***	8,880.7	11,152.8***	1,689.8***	788,641.2***	44,808.0***	7,259.1***
Rep(Loc = 3)	13.8**	72.4***	51.8	172.8***	51.2**	226.4***	393.5***	451.0***	54.8*	12,486.6	1,752.8***	182.0***
Gen(47)	34.7***	58.0***	383.9***	51.5***	46.5*	28.2*	38.5	75.5	85.1***	7,264.6	1,075.6***	293.7***
G×E(45)	20,720.7***	36.3***	455.7***	62.8***	33.0***	29.5*	34.2	67.8	74.9***	8,452	891.1***	105.9***

Note. DFF, days to first flowering; D50F, days to 50% flowering; DFP, days to first podding; GFP, grain-filling period; DPM, days to physiological maturity; PHT, plant height; TPP, total number of pods per plant; TSP, total number of seeds per plant; HSW, hundred-seed weight; PB, plant biomass; GY, grain yield; PHI, percentage of harvest index; Rep, replication; Loc, location; numbers in parentheses = degrees of freedom; Gen, genotype; G×E, genotype × environment interaction.

*

Significant at the .05 probability level.

**

Significant at the .01 probability level.

***

Significant at the .001 probability level.

TABLE 3 Variance components, heritability, genetic advance, and genetic advance as a percentage of the mean of traits under heat-stressed, pool of stressed, and nonstressed environments

[Table omitted. See PDF]

Note. DFF, days to first flowering; D50F, days to 50% flowering; DFP, days to first podding; GFP, grain-filling period; DPM, days to physiological maturity; PHT, plant height; TPP, total number of pods per plant; TSP, total number of seeds per plant; HSW, hundred-seed weight; PB, plant biomass; GY, grain yield; PHI, percentages of harvest index; EV, environmental variance; GV, genotypic variance; PV, phenotypic variance; GEV, genotype × environment interaction variance; GCV%, percentage of coefficient of genotypic variance; PCV%, percentage of coefficient of phenotypic variance; GECV%, percentage of coefficient of genotype × environment interaction variance; H^2, percentage of broad-sense heritability; GA, genetic advance; GAM , genetic advance as percentage of the mean.

Principal component and cluster analyses

At Hyderabad, a principal component analysis (PCA) biplot revealed that the associations of GY, TPP, TSP, and PHI traits were strongly positive and contributed the most to the phenotypic variability among the genotypes. These traits were negatively associated with DPM and D50F (Figure 1a). At Kobo, strong and positive associations were found among the four phenological traits (DFF, D50F, DFP, DPM) as revealed by PCA biplot. The GY, TPP, TSP, and PHI were positively associated traits. The GY and PHI were negatively associated with D50F (Figure 1b). At Debre Zeit, positive and strong associations were found among GY, PHI, TSP, and TPP traits as revealed by PCA biplot. In addition, DFP, DFF, and D50F were positively and strongly associated traits (Figure 1c).

FIGURE 1. Principal component analysis (PCA) biplots of 12 agronomic and phenological traits for two heat-stressed environments ([a] Hyderabad and [b] Kobo) and a nonstressed environment ([c] Debre Zeit). DFF, days to first flowering; D50F, days to 50% flowering; DFP, days to first podding; GFP, grain-filling period; DPM, days to physiological maturity; PHT, plant height; TPP, total number of pods per plant; TSP, total number of seeds per plant; HSW, hundred-seed weight; PB, plant biomass; GY, grain yield; PHI, percentage of harvest index. The full list of genotypes for each location can be found in Supplemental Table S1

Under heat-stressed environments, GY, PHI, TSP, and TPP vectors had high weighting towards heat-tolerant genotypes (i.e., they were important agronomic traits for the selection of genotypes for heat tolerance). Based on the two principal components, we observed more variability in the heat-stressed environments of Kobo (64%) and Hyderabad (59.1%), and less in the non-heat-stressed environment of Debre Zeit (52.4%). The PCA biplots revealed the presence of genetic variability among the genotypes that scattered between highly heat-tolerant and -sensitive categories under heat-stressed environments (Figure 1a, b).

Cluster analysis of 12 phenotypic traits at Hyderabad showed that 57 chickpea genotypes were categorized into five clusters. Cluster 5 contained only a highly heat-tolerant line DZ-2012-ck-0034 (Figure 2, Supplemental Figure S2g, Supplemental Table S1), whereas Cluster 1 contained heat tolerant genotypes including the standard heat-tolerant check, ICCV 92944, and ICCV 7102 (Figure 2). Cluster 4 contained highly heat-sensitive genotypes DZ-2012-ck-0026 and ICC 4567, whereas genotypes categorized in Cluster 3 were heat sensitive and had high PB but did not set pod at all (Figure 2).

FIGURE 2. Cluster analysis of 57 chickpea genotypes at International Crops Research Institute for the Semiarid Tropics (ICRISAT), Hyderabad, India, based on their phenotypic data. Clusters 5, 4, 3, 2, and 1 were identified by purple (V), light green (IV), orange (III), blue (II), green and red (I) colors in their respective order. The full list of germplasm accessions in each cluster can be found in Supplemental Table S1

Genomic relationships among genotypes

By using data from 5,722 SNPs, PCoA revealed that the Ethiopian chickpea genotypes grouped into two populations: one big population with five clusters and a smaller population which contained commercial cultivars. The percentage of genetic variation explained by Coordinates 1 (16.03%) and 2 (4.95%) of the PCoA was 20.98% (Figure 3). The genomic relationship results from PCoA analysis were consistent with the results of population genetic structure analysis (data not shown). The presence of a high percentage of polymorphism (82%) within landraces might infer the existence of high expected heterozygosity (H_e = 0.1, data not shown) among the landraces for heat tolerance. The second group, containing commercial cultivars, is interpreted as those with introgression from breeding lines outside Ethiopia.

FIGURE 3. Principal coordinates analysis (PCoA) of 54 chickpea genotypes based on pair-wise genetic distance. The red color depicts Population 1 (Pop1, landraces), which was further divided into five clusters based on their genomic relationships, and the green color depicts Population 2 (Pop2, commercial cultivars). The two arrows confirmed genomic distinctions of the two landraces (212589 and 41110) as they were distinct in geographical origin and altitude. Population 1: 41037, 41052, 41055, 41108, 41109, 41110, 41111, 41113, 41114, 41117, 41119, 41120, 41122, 41144, 41147, 41149, 41150, 41155, 41157, 41189, 41267, 41268, 41271, 41279, 41286, 41288, 41289, 41295, 207651, 207654, 207655, 207657, 207658, 207661, 207673, 208829, 208900, 208991, 208993, 208996, 208999, 209000, 209007, 209017, 209093, 209096, 209109, 209110, 212478, 212589, 225876, 235394. Population 2: Habru, Minjar

Furthermore, to confirm the genomic relationships among the genotypes, two landraces (accessions 41110 and 212589) were chosen as an example. Accession 41110 was collected from an extremely lowland area, Amhara region, North Shewa zone, and grown at an altitude of 1,220 m asl. Accession 212589 was collected from Amhara region, South Wello zone, and grown at an altitude of 1,600 m asl. Cluster analysis based on phenotypic data (Figure 2) revealed that accessions 41110 and 212589 were categorized under heat-tolerant and a moderately heat-tolerant genotypes, respectively (Figure 2). Similarly, based on the analysis of genomic relationships (PCoA), the two accessions were also grouped within separate clusters. Accessions 212589 and 41110 were grouped separately from most other accessions, as they are indicated by two arrows (Figure 3), indicating they are not closely related to each other.

DISCUSSION Phenotypic responses of genotypes under heat-stress environments

Assessment of whole plant growth, flowering, pod sets, and grain-filling periods under warmer environments in the field are useful preliminary heat resilient screening parameters (Devasirvatham et al., 2013; Gaur et al., 2014; Jumrani & Bhatia, 2014; Krishnamurthy et al., 2011; Summerfield et al., 1984). In addition, planting time is a useful strategy to detect heat and drought tolerances in chickpea (Devasirvatham et al., 2013; Gaur et al., 2014; Molla et al., 2018).

At Hyderabad and Werer, our two high-temperature environments, we have found more than 20 heat-sensitive genotypes that cannot flower at all and others that can flower but could not set pods. This is likely due to stress caused by high daily average temperature (≥40 °C) that could hinder flower initiations and induce sterility (Krishnamurthy et al., 2011; Summerfield et al., 1984; Upadhyaya et al., 2011). Moreover, heat-sensitive genotypes have shown delayed physiological maturity (Supplemental Figure S2f) compared with heat-tolerant genotypes under heat-stressed environments. This indicates that the transition from vegetative to reproductive phase of the plants might have been affected by heat stress, resulting in poor allocation of assimilates from the leaves into pods and seeds (Krishnamurthy et al., 2011; Upadhyaya et al., 2011).

Grain yield and yield-related traits were highly affected under all heat-stressed environments compared with non-heat-stressed environments. Particularly, at Hyderabad, there were high monthly average maximum/minimum temperatures all the way from germination (35.7/24.1 °C) to maturity (39.7/26.5 °C) (Supplemental Figure S1a), and as a consequence the majority of the genotypes yielded very little and performed poorly in yield-related traits (Table 1). Consistently, at Werer, most genotypes had poor pod sets, with reduced TPP and TSP compared with nonstressed settings (data not shown). In agreement with the present study, Summerfield et al. (1984) found that the greater the proportion of the reproductive period occurs on hot days (≥35 °C), the lower the yields. Plants transferred to temperatures above 35 °C at 50% flowering did not produce pods. The negative effects of heat stress on these traits might also be associated with the reduction in reproductive efficiency, with low levels of carbohydrates and growth regulators released in plant sink tissues and induced early maturity (Jumrani & Bhatia, 2014).

At Hyderabad, both the vegetative and reproductive stages were exposed to heat stress (monthly average temperatures ≥ 35.7 °C; Supplemental Figure S1a). The effects of heat stress on the vegetative stages of chickpea were drying of leaves, shoots, branches, stems, and whole plants, exhibited by stunted growth, delayed reproductive stages, yellowing of leaves, branches, flowers, and pods, and decrease in the number of plants per plot (Supplemental Figure S2c, d). This is because high temperatures might affect some of the physiological processes and might reallocate some of the resources of chickpea genotypes to withstand the heat stress. A single highly heat-tolerant genotype (DZ-2012-ck-0034) that possessed high pod set, TPP, TSP, and GY was found under all heat-stressed environments both in Ethiopia and India compared with the standard and local checks (Supplemental Figure S2i, Figures, 1a, b, and 2). In agreement with these findings, Devasirvatham (2012), Molla et al. (2018), and Upadhyaya et al. (2011) have found very few heat-tolerant chickpea genotypes under greatly heat-stressed conditions.

In the present study, ANOVA revealed highly significant differences (P ≤ .001) among the genotypes for most traits under stressed and nonstressed environments. In line with this, Devasirvatham (2012), Kanouni et al. (2009), Krishnamurthy et al. (2011), Saeed and Darvishzadeh (2017), and Upadhyaya et al. (2011) reported significant differences (P ≤ .01) among chickpea genotypes for most of the traits they studied under stressed and nonstressed conditions. Our results suggest that there are heat-tolerant accessions available in Ethiopian germplasm that can be harnessed in breeding efforts, addressing this increasing abiotic challenge.

Performance of variance components under heat-stressed environments

In the current study, environmental, genotypic, and phenotypic variances became higher for GY and yield-related traits among genotypes at Kobo than at Hyderabad. Geographic location, altitude, the extent of heat stress, and the composition and number of genotypes might account for this difference. In agreement with this study, wide phenotypic variations were detected for GY and yield-related traits among 291 chickpea accessions (Basu et al., 2019). Various authors have observed phenotypic, genotypic, and genotype × environment interactions for traits studied in chickpea genotypes (Krishnamurthy et al., 2011; Paul et al., 2018; Saeed & Darvishzadeh, 2017; Upadhyaya et al., 2011). This study also found large coefficients of phenotypic and genotypic variances for GY and yield-related traits in all heat-stressed environments than in nonstressed environments; this is because the trait value of an individual genotype had large deviations from the mean value of the trait as depicted by PCA biplot analysis (Figure 1a, b) compared with the nonstressed environment (Figure 1c). These results were in agreement with the reports of Paul et al. (2018), where they found a high genotypic coefficient of variation and phenotypic coefficient of variation for traits assessed under heat-stressed environments. Like other stresses, high temperature conditions may cause greater phenotypic variation than less stressful environments, which has been noted in natural populations (Stanton et al., 2000). Further investigation of the impacts of multiple stresses on phenotypic variation in chickpea germplasm is needed.

Genomic relationship among genotypes and their patterns of distribution

The present study used 5,722 genome wide SNPs data to determine genome-wide relationships and patterns of distribution among 54 chickpea genotypes. In the landrace group (population 1), out of 5,722 SNP loci; 572 SNPs (0.1%) of the loci had three alleles per locus.

The landraces of chickpea in Ethiopia had higher genetic variation (H_e = 0.1) compared with the genetic diversity (H_e = 0.02–0.04) of landraces of chickpea in India, Iran, and Pakistan, as reported by Sani et al. (2018). This might be important for agroclimatic adaptation in their respective agroecological zones. In line with this, Sani et al. (2018) reported that with this low amount of genetic diversity (0.02) in Pakistan, landraces of chickpea have longstanding cultivation histories under heat and drought environmental conditions. Penmetsa et al. (2016) reported that cultivated chickpea accessions had low genetic diversity (H_e = 0.065) compared with the wild species, Cicer reticulatum Ladiz. (H_e= 0.332) and Cicer echinospermum P.H. Davis (H_e = 0.301), using a set of SNPs from a GoldenGate assay. Similarly, in the present study, the genetic variation among 52 chickpea genotypes (H_e = 0.1) is approximately three times less than the genetic diversity of these wild species (H_e = 0.30–0.33).

Sani et al. (2018) revealed low genome wide differentiation among Pakistan chickpea landraces even though significant relationships were found between genetic, geographic locations, and climatic factors. Similarly, in this study, our PCoA analysis categorized 54 chickpea genotypes into two groups: desi landraces (52) and commercial cultivars (2, Habru and Minjar) (Figure 3). In the same scenario, 66 (Teshome, 2012) and 93 (Kujur et al., 2015) chickpea accessions were categorized into three distinct populations: wild, kabuli, and desi; 291 chickpea accessions were categorized into two major populations, desi and kabuli (Basu et al., 2019). Likewise, Penmetsa et al. (2016) reported that 247 chickpea accessions were clustered into three populations (K = 3): desi, kabuli, and wild.

Importantly, we found that the few heat-tolerant genotypes in our phenotypic study did not cluster closely together, suggesting that they are not closely related. This means that this critical trait likely segregates across the Ethiopian gene pool. It also indicates that different tolerant accessions may do so by different genetic mechanisms. If this proves to be the case, crossing these accessions could create new varieties with superior heat tolerance. Stacking heat tolerance genes will help increase the climatic resilience of this important crop.

Conclusions and implications to improve chickpea for heat tolerance environments

Choosing appropriate experimental field sites and planting dates under heat-stressed environments and characterizations of phenological and agronomic traits are crucial to select and identify heat tolerant genotypes. Understanding genome-wide relationships among genotypes via SNP markers helps to assess the clustering patterns of genotypes and their relationships. Heat stress is a function of plant genotypes, high temperatures, plant phenology, and soil types. The experimental materials we have examined demonstrate that the production range of chickpea can be stretched in environmental heat gradients and using existing trait variability. Our field observations revealed that heat stress (≥35 °C) mainly affects the reproductive stages. If plant germination occurs at high temperatures (≥35 °C), heat stress has negative effects on vegetative stages of chickpea genotypes. Grain yield, TPP, TSP, and pod set are the best traits for screening of heat-tolerant and -sensitive genotypes under heat-stressed environments. A total of 5,722 genome-wide SNPs were used to determine genome-wide relationships among 54 chickpea genotypes. Principal coordinate analysis has found five different clusters within the landrace group. Each genomic cluster might possess specific adaptive genes to their respective geographical locations. The genotypes that showed better performance in this experiment can be used in crossing blocks to improve the productivity of chickpea in high temperature agroecological regimes.

ACKNOWLEDGMENTS

The authors are very grateful to the United States Agency for International Development (USAID); and Feed the Future Chickpea Innovation Laboratory at the University of California Davis for the full financial support of phenotyping as well as genotyping works that were done in Ethiopia and India, and at the University of California-Davis (USA). Moreover, the authors are thankful for the provision of germplasm accessions by the Ethiopian Biodiversity Institute (EBI), Melkasa and Debre Zeit Agricultural Research Centers, and ICRISAT, Hyderabad, India. They also acknowledge Debre Zeit, Werer, Kobo, and Hyderabad (ICRISAT) research centers for the provision of experimental sites.

AUTHOR CONTRIBUTIONS

Tsegaye Getahun: Data curation; Formal analysis; Investigation; Methodology; Software; Visualization; Writing-original draft. Kassaye Negash: Data curation; Resources. Peter L. Chang: Data curation; Methodology; Software; Validation. Eric von Wettberg: Conceptualization; Formal analysis; Investigation; Methodology; Software; Supervision; Writing-original draft; Writing-review & editing. Noelia Carrasquilla-Garcia: Data curation; Methodology; Validation; Visualization. Pooran M. Gaur: Conceptualization; Investigation; Resources; Supervision; Visualization; Writing-review & editing. Asnake Fikre: Conceptualization; Investigation; Methodology; Supervision; Writing-review & editing. Teklehaimanot Haileslassie: Resources; Supervision; Writing-review & editing. Douglas Cook: Conceptualization; Funding acquisition; Investigation; Project administration; Resources; Supervision; Writing-original draft; Writing-review & editing. Kassahun Tesfaye: Conceptualization; Investigation; Project administration; Resources; Supervision; Writing-review & editing.

CONFLICT OF INTEREST

All authors declared no conflict of interest.