Association analysis of biotic and abiotic stresses resistance in chickpea (Cicer spp.) using AFLP markers

ABSTRACT Chickpea suffers from different biotic and abiotic stresses, which has led to considerable decreases in seed yield. Since the economic value of a cultivar is determined by its phenotypic characteristics, good knowledge of the genetic structure and identification of genomic loci associated with desired traits will facilitate crop breeding. Amplified fragment length polymorphism (AFLP) markers associated with different agro-morphological characteristics, cold and drought tolerance as well as Ascochyta blight resistance were identified in chickpea by using 44 genotypes comprising cultigens, landraces, internationally developed improved lines and wild relatives, evaluated in a randomized complete block design with three replications at the Urmia rainfed research station. Using six AFLP primer combinations, 64 clear and reproducible markers were developed. The polymorphism information content values calculated for each primer pair ranged from 0.155 (EcoR1-ACC/Mse1-CAG) to 0.270 (EcoR1-ACC/ Mse1-CTG) with an average of 0.237. Analysis of molecular variance (AMOVA) revealed that 90% of the total variance was due to differences within populations and 10% due to differences among populations. Using a general linear model and applying multiple testing corrections, 24 AFLP markers associated with genes controlling the studied traits were identified. Some identified markers were associated with more than one trait. The identified markers could be of interest in marker-assisted selection in chickpea breeding programmes.


Introduction
Chickpea (Cicer arietinum L.) as a third major cool-season grain legume in the world is grown in 13.9 million hectares with 13.7 million ton production as of 2014 (http:// faostat.fao.org). Stresses caused by high or low temperatures lead to crop yield losses due to damage to reproductive organs [1], floral bud development [2], pollen viability and fertilization [1], pod production, seeds per plant, seed yield [3], seed filling [4] and seed composition [5]. Chickpea may be exposed to low temperatures and excess water early in the growing season (especially if it is planted in the fall) or high temperatures and inadequate water at the end of the growing season [6]. Tolerance to low temperatures in chickpea primarily depends on the weather conditions and the sowing date. High levels of cold tolerance are known in local cultivars and related species such as Cicer bijugum, Cicer reticulatum, Cicer echinospermum, Cicer pinnatifidum and Cicer judaicum [7]. Wery [8], by evaluating 29 chickpea genotypes for frost resistance in Montpellier, France, identified three frost-resistant genotypes, FLIP 81-293C, FLIP 82-127C and FLIP 82-128C. A key recurring problem in the screening for abiotic stress tolerance is environmental heterogeneity [6].
Drought, as another major abiotic stress in many parts of the world [9], causes heavy production losses, 30%-100% [10] in cool-season food legumes [6,11]. Stoddard et al. [6] reviewed the screening techniques in cool-season food legumes and listed 25 drought-resistant chickpea genotypes. Nevertheless, these genotypes are far from possessing the desirable characteristics to grow in farmer fields due to being small seeded and susceptibile to Ascochyta blight [12]. The quantitative inheritance of drought tolerance and its interaction with the environment are the major challenges that exist in identifying drought tolerance genotypes in chickpea [13].
Ascochyta blight (causal agent Ascochyta rabiei) is considered to be a major constraint for chickpea production in cool and wet regions of the world [14]. About 50% yield reduction has been reported to have occurred due to Ascochyta blight epidemics in 1979-1982 [15]. Most of the cultivated chickpea varieties are only partially resistant to this pathogen [16]. Pyramiding quantitative trait loci (QTLs) using marker-assisted selection (MAS) is generally considered an appropriate strategy for improving the tolerance levels of chickpea in breeding programmes. Comparing to other crop species, such as cereals, the advances in the development of genomic resources in the leading legume crops of the semi-arid tropics, especially chickpea (Cicer arietinum), have been very slow [17]. However, these programmes have significantly gained by integrating MAS within empirical selection [18].
The power of association analysis to identify and characterize loci/genes associated with different complex traits is highly affected by admixtures of populations [39]. It is generally considered that good understanding of the population structure in an association panel is necessary to avoid identifying false positive correlations between markers and traits [40]. The objectives of the present study were to characterize the population structure within chickpea genotypes and to identify AFLP markers associated with different agro-morphological characteristics, cold and drought tolerance as well as Ascochyta blight resistance in chickpea using association analysis.

Plant materials
A population of 44 chickpea genotypes comprising cultigens, landraces, internationally developed improved lines and wild relatives was used for association analysis ( Table 1). The genotypes of the association panel were prepared from six geographically diverse regions including International Center for Agricultural Research in the Dry Areas (ICARDA) (28 genotypes), Turkey (five genotypes), Iran (eight genotypes), USSR (one genotype), Egypt (one genotype) and Spain (one genotype) with Latitude of 5 30 W, 40 02 N to 45 09 E, 37 12 N ( Figure 1).

Genomic DNA isolation, PCR conditions and electrophoresis
Young leaf tissues of two-week-old seedlings were cut and were subjected to vacuum freeze-drying for dehydration. Genomic DNA was extracted from leaf tissues using the modified method described by Dellaporta et al. [44]. An association panel was fingerprinted with six AFLP primer combinations (Table 2). To identify the most informative AFLP primer combination, the means of the polymorphism information content value, the effective multiplex ratio and the resolving power of markers (Rp) were considered [45]. The DNA concentration in each sample was adjusted to 50 ng mL ¡1 . To 6.0 mL of diluted DNA in suitable 1.5-mL tubes, 0.25 mL EcoRI, 0.25 mL MseI, 4.0 mL 10X buffer and 9.5 mL ddH 2 O were added. The tubes were incubated at 37 C overnight to digest. The DNA digestion step was confirmed by running the products in a 1% agarose gel. In the next step, to 20 mL of digested DNA, 0.5 mL EcoRI adapter, 0.50 mL MseI adapter, 0.50 mL adenosine-5'-triphosphate (ATP) (10 mmol L ¡1 ), 1.0 mL T4 DNA ligase, 2.0 mL T4 buffer and 0.5 mL distilled water were added. All tubes were incubated at 37  Four microliters of reaction products were mixed with 1 mL of formamide dye (98% formamide, 10 mmol L ¡1 ethylenediaminetetraacetic acid (EDTA), 0.25% v/v bromophenol blue and 0.25% v/v xylene cyanol FF) and denatured at 95 C for 4 min and quenched on ice. The products were resolved in a 1.2% (w/v) metaphor agarose gel (Sigma-Aldrich, St. Louis, MO) in 0.5X Tris-Borate-EDTA (TBE) buffer, stained with ethidium bromide (1.0 mg mL ¡1 ) and photographed under ultraviolet (UV) light using a Gel-Doc image analysis system (Gel Logic 212 PRO, Woodbridg, NJ). The standard 50-900 bp AFLP DNA ladder (Fermentas, Carlsbad, CA) was used in the flanking lanes. Compared to polyacrylamide-gel electrophoresis, or automated analysis, agarose-gel electrophoresis is the most-appropriate and safe technology for routine analyses of these types of markers.

Population structure and association analysis
The AFLP data were scored as dominant markers for each locus. Population structure was analysed using a model-based Bayesian approach in the software package Structure 2.3.4 [46]. Five independent runs were performed, setting the number of sub-populations (K) from 1 to 10, the burn in time and MCMC (Markov Chain Monte Carlo) replication number both to 100,000, and a model for admixture and correlated allele frequencies.
The K value was determined by the log likelihood for each K; Ln P (D) = L (K), and an ad hoc statistic DK based on the second-order rate of change in the likelihood (DK) [47]. The appropriate value of K was selected when the estimate of Ln P (D) reached a minimum stable value. Inferred ancestry estimates of individuals (Q-matrix) were derived for the selected sub-population [46]. Association analysis was performed to analyse marker-trait association by the structured association approach using ancestry coefficient (Q-values) estimates as covariate in a general linear model (GLM) function using TASSEL 2.1. Multiple testing corrections were performed by adjusting maker probability values for multiple test runs by a permutation test in the TASSEL 2.1 software.

Morphological diversity
Analysis of variance revealed significant differences among the studied genotypes for most of the studied traits ( Regarding the tolerance to Ascochyta blight disease, 50% and 83.3% of desi and kabuli types were tolerant, respectively. In terms of drought stress tolerance, 53.8% and 33.3% of the two types were tolerant, respectively (Online Supplemental Data).
The potential seed yield of chickpea, about 6 t ha ¡1 , is far below the world average yield (0.98 t ha ¡1 ) [48]. This is due to several complex field environmental factors such as the global climate change [49], rainfed sowing, and poor farm management, lack of well adapted cultivars with low levels of tolerance to different biotic and abiotic stresses. High genetic variability was observed among the chickpea genotypes for the studied agromorphological traits as well as resistance to biotic and abiotic stresses. Our findings showed that considerable part of genotypes involved in the current minicore germplasm were drought tolerant. Toker et al. [50] reported that, under drought stress conditions, all accessions of perennial wild Cicer species were significantly superior to the annual wild and cultivated chickpeas, including the best tolerant chickpea, ICC4958. The results showed that the chickpea sowing time in dry lands is a major parameter in drought tolerance by accelerating the plant life cycle before dramatic water deprivation. Drought escape through short and/or better timed life cycles, matched with large total biomass, is the primary form of crop adaptation to low rainfall Mediterranean-type conditions [51]. Therefore, achieving successful production under terminal drought stress depends on the development of short-season (earliness) varieties that enable the crop to escape severe soil-water deficits [52]. Drought avoidance as an alternative mechanism corresponds to the maintenance of high tissue water potential and consists in mechanisms that both reduce water loss from plants, via stomatal control of transpiration, and also maintain water uptake through an extensive and prolific root system [53]. In the present study, we took advantage of both mechanisms, so the suitable drought tolerance varieties showed acceptable yield at stressed conditions. The results showed that Cicer genotypes are well adapted to germination and seedling establishment in moderate conditions, followed by vegetative growth in cooler conditions which are typical following autumn or winter sowing in temperate Mediterranean and cool West Asia regions. The cold tolerance scores of Cicer species indicated that the level of cold tolerance in all accessions of C. reticulatum and C. echinospermum were superior to those of well-known cold tolerant varieties, ILC 8262 and ILC 8617. These results were in agreement with the data reported by Robertson et al. [54] and Singh et al. [55]. They found that most of the accessions of C. bijugum, C. echinospermum and C. reticulatum are tolerant to cold and have significantly higher levels of cold tolerance than cultivars. The results indicated that some desi types, such as Sel 95TH 1716 and Sel 95TH 1744, as well as kabuli types, such as FLIP 93-261C and x03TH21, showed higher cold tolerance during the early seedling stage and withstood temperatures of ¡15.6 C without snow covering, under dry cold conditions. The results also indicated that, for colder areas of highlands, genotypes adapted for winter or autumn sowing could be obtained by the use of genes for cold tolerance from the wild Cicer species. In a preliminary screening for cold tolerance, Toker [7] found that 11 accessions of C. reticulatum and 4 accessions of C. echinospermum were highly cold tolerant and superior than that the best cold tolerant cultigens, ILC 8262 and FLIP 93-53C. This was supported by our findings, and the findings of Singh et al. [56].
The results showed that all wild Cicer genotypes and landraces as well as some kabuli types were resistant to Ascochyta blight. Most of desi types showed susceptibility to disease. Reddy and Sing [15] evaluated 9574 desi and 3836 kabuli germplasm accessions from the world collection of chickpea for resistance to Ascochyta blight both in vegetative and podding stages and reported that 11 kabuli and 6 desi accessions were resistant. A major part of reports focus on strategies to minimize the damage caused by Ascochyta blight [16]. Factors such as using Ascochyta blight-free seeds, one to two years of non-host crops for warm and wet areas, 3-4 years crop rotation for cold and dry areas, optimum sowing date, deep sowing, optimizing plant density, balanced nutrition and alternative sowing patterns should be considered as useful means of reducing this disease [16].

Molecular diversity and population structure
Using 6 AFLP primer combinations, 64 clear bands were generated. The number of polymorphic markers varied from 6 (EcoRI-ACC/ MseI-CAG) to 46 (EcoRI-AAC/ MseI-CAG) for primer combinations. Detailed information on the levels of genetic diversity of genotypes is available elsewhere [45].
Analysis of population structure was carried out by using a model-based clustering approach that assigns individuals into one or more sub-populations probabilistically based on the allele frequencies detected at different loci [46]. The structure clustering results from K = 2-10 are shown in Figure 2. As indicated in Figure 2, the K values (number of sub-populations) steadily kept on increasing until K = 4, indicating that the 44 genotypes in this study are clustered into four sub-populations. Genetic variation among the four identified sub-populations was tested using F-statistics, estimated from   pairwise comparisons as a measure of genetic distance between sub-populations. The F-statistics values between sub-populations were significant (P = 0.01) and ranged from 0.02 to 0.96, supporting the existence of genetic structure. Of all genotypes, 65.9% were assigned into the corresponding subgroups, and the remaining ones were categorized into the 'mixed' genotypes based on their Q-values (Table 1 and Figure 3). In our association panel based on AFLP markers, the r 2 values ranged from 0.082 to 0.333 with an average of 0.027 and D' ranged from 0.001 to 0.999 with an average of 0.356 ( Figure 4). Of 576 AFLP marker pairs, 7.76% showed a significant level of LD (P < 0.05).

Association analysis
By applying the GLM procedure and multiple testing corrections in the TASSEL software, one AFLP marker was identified to be associated with genes controlling the stem number; two markers, with the days from the first effective raining after sowing to 50% flowering and two ones, to 90% maturity; two markers, with partial tolerance to Ascochyta blight; one marker, with the 100-seed weight and two ones, with the yield ( Table 4). None of identified markers for plant height, cold and drought tolerance were statistically significant. Among the identified AFLP markers, six markers were common, whereas 12 ones were specific.
The significant level of LD found in the present study supports the idea of using natural populations for genomic studies. Being a predominantly self-pollinating species, chickpea is expected to show higher LD extent, 150-200 Kb [36] and 500-600 Kb [57], than outcrossing species like maize with (2 Kb) LD extent [58]. Several factors which cause LD content increasing include low recombination rate, small population size, population admixture, inbreeding, genetic isolation between lineages, population subdivision, genetic drift and epistasis [22,59]. The admixture of populations is the most important factor influencing marker-traits associations [39]. Therefore, good understanding of the population structure is a prerequisite in association mapping and necessary to avoid identifying false-positive correlations between markers and traits [40,60]. Using polymorphic bands produced by six AFLP primer combinations, the association panel subdivided into seven sub-populations. The GLM was used to identify genes controlling some agro-morphological traits in the studied chickpea germplasm mini-core collection. Heretofore a considerable number of QTLs have been identified for desirable agronomic traits in chickpea on different linkage groups (LG 2, 3, 4, 6 and 8) [61]. For instance, Iruela et al. [62] identified three QTLs for resistance to Ascochyta blight on linkage group 4 (LG4) (QTL AR1 and QTL AR2) and 2 (QTL AR3). Cho et al. [63] reported two major QTLs on LG 2 close to the GA16 and TA37 loci controlling the resistance to Ascochyta rabiei pathotype I, as well as two QTLs for resistance to pathotype II on LG 4 such that one is linked to CaETR or GAA47 and the other is linked to TA72/ScY17. Bian et al. [64], by comparing three chickpea LGs comprising QTLs for tolerance to A. rabiei, found that LG3 correspond to the sub-centromere region of chromosome C in chickpea.
LG8 is located on the long arm of chromosome H (W-Ca-LG8) and LG4 is located in the sub-centromere region of chromosome B (W-Ca-LG4). Based on bibliography, the majority of AB-resistance QTLs were reported mainly on two LGs, CaLG02 and CaLG04 [61,65]. Two important QTLs (Q3-1 and Q1-1) were reported for drought tolerance from ILC 588 £ Table 4. QTLs detected for the studied agro-morphological traits in chickpea using the GLM procedure. ILC 3279 population which was located on LG3 and LG1 [66,67]. In order to understand the complex nature of drought tolerance in chickpea, Thudi et al. [68] carried out an association mapping study and identified 15 markers significantly associated with five root traits (root dry weight, root length density, root surface area, root volume and rooting depth). Among them, seven markers showed significant relation with a single trait (rooting depth) and two markers (NCPGR7 and DR ç 237) showed associations with more than one trait. They concluded that these two markers were associated with co-localized/pleiotropic QTLs. In the case of cold tolerance, we were not able to identify any reported QTLs. The association studies in chickpea have revealed that the linkagebased QTL analyses can be complemented by LD-based association studies [68]. The markers identified in this study could be transferred to sequence characterized amplified region (SCAR) markers to increase the efficiency of MAS in chickpea breeding programmes for biotic and abiotic stresses.

Conclusions
High genetic diversity was observed for the studied agro-morphological traits in the association panel. Using a GLM, several AFLP markers were detected for the studied traits. The identified markers can be transferred to sequence characterized amplified region markers for facilitating MAS. Detection of molecular markers associated with genes controlling different traits could increase the efficiency of MAS in chickpea breeding programmes for biotic and abiotic stresses.