Identification of SNP Markers Associated with Iron and Zinc Concentrations in Cicer Seeds.

BACKGROUND: Cicer reticulatum L. is the wild progenitor of chickpea Cicer arietinum L., the fourth most important pulse crop in the world. Iron (Fe) and zinc (Zn) are vital micronutrients that play crucial roles in sustaining life by acting as co-factors for various proteins. AIMS AND
OBJECTIVES: In order to improve micronutrient-dense chickpea lines, this study aimed to investigate variability and detect DNA markers associated with Fe and Zn concentrations in the seeds of 73 cultivated (C. arietinum L.) and 107 C. reticulatum genotypes.
METHODS: A set of 180 accessions was genotyped using 20,868 single nucleotide polymorphism (SNP) markers obtained from genotyping by sequencing analysis.
RESULTS: The results revealed substantial variation in the seed Fe and Zn concentration of the surveyed population. Using STRUCTURE software, the population structure was divided into two groups according to the principal component analysis and neighbor-joining tree analysis. A total of 23 and 16 associated SNP markers related to Fe and Zn concentrations, respectively were identified in TASSEL software by the mixed linear model method. Significant SNP markers found in more than two environments were accepted as more reliable than those that only existed in a single environment.
CONCLUSION: The identified markers can be used in marker-assisted selection in chickpea breeding programs for the improvement of seed Fe and Zn concentrations in the chickpea.

INTRODUCTION

The chickpea, a member of family Fabaceae, is one of the most important pulse crops in the diet of millions of people [1]. It is cultivated in over 14.8 Mha area in Asian and sub-Saharan regions in the world with an annual production of 14.24 million tons. Global climate change and the growing human population have resulted in increased hunger and malnutrition [2], and food and nutritional security has become one of the major global issues [3]. Increasing vitamin and mineral content of crops, such as chickpea, is essential for humans to sustain their health because over three billion people in developing countries are suffering from micronutrient deficiencies [3], which has led chickpea breeding programs to focus on biofortification [4].

The chickpea is a member of West Asian Neolithic crop assemblage, associated with the origin of agriculture in the

Fertile Crescent around 10,000 years ago with the oldest archaeological evidence from 7500 B.C. Cicer reticulatum Ladiz. is a wild progenitor of the chickpea, and Cicer arietinum L. is the only cultivated species in genus Cicer. Since both species are interfertile, C. reticulatum is used as a primary gene pool [5]. Wild Cicer species contain genes for disease resistance, drought tolerance, and desirable yield components, which increase their significance [6]. Supporting this knowledge is a study reporting that transferring some genes from C. reticulatum into cultivated chickpea by the pedigree method increased the seed yield [6]. Two different forms of cultivated chickpea are Desi and Kabuli. The former has small, dark-colored seeds with an angular shape while the latter has large, beige seeds with an owl-head shape [7].

Iron (Fe) and zinc (Zn) are vital micronutrients and play crucial roles in sustaining life by acting as co-factors for various proteins. Fe and Zn deficiencies cause critical physiological disorders, such as low immunity, tissue hypoxia, anemia, stunted growth, and dwarfism [3, 8]. According to the World Health Organization (2002), Fe and Zn are the most widespread micronutrient deficiencies in developing countries, and 66-80% [9] and 50% [10] of the global population are affected by Fe and Zn deficiencies, respectively. Biofortification is the process of increasing the concentration of essential minerals and vitamins in staple crops through conventional plant breeding and genetic modifications in a cost-effective and sustainable way [11]. Through biofortification, essential nutrients, such as Fe and Zn which are not available to all populations can be supplied, thus reducing micronutrient malnutrition [12]. According to economic analyses, Zn biofortification of wheat in Turkey and Fe biofortification of rice in South Asia and Bangladesh increase the cost-to-benefit ratio [13]. In addition, Fe and Zn biofortification help to significantly reduce the prevalence of diseases associated with Fe and Zn deficiency in India [13]. In order to advance breeding strategies, identifying the genetic basis of different traits in plants is an important step, and genetic sources of phenotypic variation are a primary focus. To date, linkage and quantitative trait loci (QTL) analyses have been the main methods for obtaining information on the genetic bases of several important traits in the chickpea, but they are limited in terms of allelic diversity and genomic resolution. Genome-wide association study (GWAS), which was developed to accelerate molecular breeding for the advancement of agronomic traits, eliminates the limitations of QTL mapping. Genome-wide association study is based on linkage disequilibrium and detects QTL in diverse collections of crop germplasm. However, since a set of unrelated genotypes shows restricted linkage disequilibrium, a great number of markers should be used to sufficiently cover the whole genome [1, 14]. High throughput genotyping techniques, such as genotyping by sequencing (GBS) provide the required marker coverage by generating a large number of single nucleotide polymorphisms (SNPs) [2, 14, 15]. Genotyping-by-sequencing is a cost-effective approach that allows simultaneous discovery and genotyping of a large number of SNPs with high quality and low cost, and also simpler and faster than previous methods. In addition to the requirement of a lower amount of input DNA, GBS is available for species with complex genomes with or without reference sequences [16].

The first association mapping study related to major genomic loci regulating Fe and Zn concentration in the chickpea was conducted [1] with 1,536 SNP markers by Illumina GoldenGate assays. The authors reported eight SNPs associated with Fe and Zn concentration using 94 chickpea genotypes, of which 23 were Desi and 71 were Kabuli. In a later study, 24,620 SNPs were obtained by the GBS method, 16,591 SNPs used in GWAS analysis and 16 genomic loci/genes were defined related to seed Fe and Zn concentrations using GWAS [3]. The researchers also confirmed 11 of these genes by QTL mapping. However, to date, no association mapping study has been conducted to investigate Fe and Zn concentrations in the chickpea seeds using both C. arietinum and wild progenitor C. reticulatum genotypes. In contrast to previous studies, the current research used a higher number of genotypes to obtain increased genetic variation and a higher number of SNP markers to deepen the analysis. The objectives of this study were to identify the population structure, genetic diversity and the markers associated with Fe and Zn concentrations in the chickpea. The knowledge gained from this research will allow scientists to increase the Fe and Zn concentrations of the chickpea with various molecular breeding techniques.

MATERIALS AND METHODS

Plant Materials and DNA Extraction

A total of 180 genotypes, 107 C. reticulatum and 73 C. arietinum, were used for genetic evaluation (Supplementary Tables and ). C. reticulatum genotypes were supplied by the USAID project [17] and C. arietinum genotypes were supplied by the GeneBank of the Aegean Agricultural Research Institute in Izmir, Turkey. The seeds of each genotype were sown in 5x1 m2 plots in the field. The trial design consisted of three blocks, each of 216x5 m2. Experimental soil tillage was done by plow 30 cm in depth followed by disk harrow, and seeding was applied using a rototiller. The field experiments were conducted in two locations in Turkey (the experimental field of the Department of Field Crops at Ege University in Bornova, Izmir and Sanliurfa), and for two years (2015 and 2016 growing seasons, sowing dates of Bornova were 15 November for the first year and 20 November for the second year, and sowing dates of Sanliurfa were 22 November for the first year and 28 November for the second year). Fertilizer was applied at a rate of 40 kg/ha for Nitrogen and 70 kg/ha for Phosphorus. The experiment was set up in a complete randomized block design (CRBD) with three replications in each location and year. The genotypes were accepted as a fixed variable, and years and locations constituted the random variables. The weather data and soil profile of Izmir and Sanliurfa are given in Supplementary Table . In terms of weather, Sanliurfa is warmer, but the total monthly rainfall and average rainy days are higher in Izmir. In terms of the soil profile, pH, total salt (%), organic matter and clay are very similar in Izmir and Sanliurfa, but there are considerable differences in receivable elements.

For DNA extraction, young leaves from four- to six-week-old seedlings were harvested and placed in Eppendorf tubes in liquid nitrogen, then stored at -80°C until processed further. The leaf tissues from each individual were ground to a fine powder in liquid nitrogen with a tissue lyser (Technogene, Turkey). Then, genomic DNA was extracted from 100 mg of leaf tissue using a DNA isolation kit (Plant Genomic DNA Purification Kit Cat No: K0792 Fermentas, Germany). The isolated DNA was visualized by electrophoresis on 1% agarose gel, and the purified DNA was quantified by a Nanodrop ND-1000 spectrophotometer (Thermo Sci. Co.).

Micronutrient Analysis

The seeds of each genotype were washed first with tap water, then with ddH2O water. Following this step, the seeds were dried in an oven at 65°C for 24 hours. The dried seed samples were ground using an analytic mill (IKA, A11, Staufen, Germany) in accordance with the protocol described by Kacar and Inal [18]. All procedures were repeated three times (technical replications). Two g of each ground sample was placed into a 150 mL Erlenmeyer flask and 24 mL of HNO3-HClO4 (4:1) was added; then, the Erlenmeyer was placed on a hot plate (250°C) for three hours to decompose the samples. An atomic absorption spectrometer (Varian, SpectrAA 220/FS, California, USA) was used to determine the Fe and Zn values in all seed samples. The Fe and Zn values were converted into mg/kg concentrations in seeds. The genotype-environment interactions for the Fe and Zn concentrations were identified by analysis of variance (ANOVA) using R statistical software. These concentrations were evaluated as dependent variables. The genotypes, locations, and years were considered as the fixed effects and replications as random effects [1].

Molecular Marker Analysis

SNP analysis was implemented by the GBS approach in accordance with the procedure described by von Wettberg et al. [17]. Firstly, genomic DNA was extracted and diluted to 10 ng/µl. Genomic DNA was digested using HindIII and NlaIII restriction enzymes to create a 4-bp overhang to promote efficient adapter ligation. After the adapters were ligated to the restriction sites, purification and fragment size selection processes were applied using AMPure Beads, and amplification was performed by PCR. The sequencing of the amplicons was implemented using the Illumina HiSeq4000 sequencing system using a 100 bases reads protocol at the University of California Davis Genomics Facility. The raw sequencing data was analyzed by the GATK pipeline as described by Mc-Kenna et al. [19] by mapping them according to the C. arietinum CDC Frontier reference by Varshney et al. [20] in Haplotype Caller program. Then, the data was filtered according to the recommendations of the GATK Best Practices [21]. All sequencing data can be found in the NCBIBioProjects PRJNA353637 and PRJNA416006.

Diversity and Population Structure Analysis

The population structure analysis was performed as described earlier [22] using STRUCTURE software [23] by adopting the Markov Chain Monte Carlo approach. For this purpose, the number of subpopulations (K) ranged from 1 to 8 based on admixed and correlated allele frequencies, and each was repeated 10 times. Each number of subpopulations was characterized by a set of allele frequencies at each locus. For each run, burn-in period iterations and the number of iterations were set to 100,000. The run with maximum likelihood was utilized to appoint individual genotypes into sub-populations. Population structure analyses were carried out with all 180 genotypes (Analysis A) and also separately for each species (C. reticulatum (Analysis B) and C. arietinum (Analysis C) genotypes). Principal component analysis (PCA) was conducted in R statistical software [24] to reveal the spatial representation of the relative genetic distances between the genotypes. A neighbor-joining phylogenetic tree (NJTree) based on standard genetic distance was also drawn in the same software.

Association Analysis

The marker genotype data from the SNP markers obtained from GBS were used in the analysis. The SNPs with minor allele frequencies (<5%) and those more than 10% missing data were excluded. Association analysis was performed on the Fe and Zn concentrations according to a mixed linear model (MLM) (Q+K) method using TASSEL software [25]. For this purpose, phenotypic data, marker genotype data, Q matrix and Kinship matrix (obtained from TASSEL software, given in Supplementary Tables and ) were utilized in the same software. The population structure Q matrix (obtained from STRUCTURE) was used to define the sub-population parentage for each line in the analysis. False discovery rate (FDR) correction [26] and Bonferroni’s correction [27] were applied for the Fe and Zn concentrations. A significant association between a marker and a phenotypic trait was explained by the p-value of ≤ 0.01 demonstrated by the Manhattan plots and quantile-quantile (Q-Q) plots. The magnitude of the QTL effects was evaluated by the R2 ≥0.1 marker. The associated markers were examined using three datasets: analysis A including both C. arietinum and C. reticulatum, analysis B consisting of only C. reticulatum genotypes, and analysis C comprising only C. arietinum genotypes. Association mapping analysis was conducted on four environments and the general mean of the Fe and Zn concentrations separately. The markers which exceeded the threshold of (FDR) for more than two environments were accepted as more reliable than those that exceeded this threshold only in a single environment.

Candidate Gene Analysis

The basic local alignment search tool (BLAST) was used on the National Center for Biotechnology (NCBI) database to identify candidate genes with putative functions based on the sequences of significant SNP markers associated with Fe and Zn concentration in Cicer seeds. In addition, for all SNPs, the genes within 100K to the left and right position on the chromosome were included in the analysis (Supplementary Table ).

RESULTS

Fe and Zn Concentrations in Seeds

The frequency distributions for the Fe and Zn concentrations of the 180 genotypes are presented in Supplementary Fig. (). The means and ranges of Fe and Zn concentrations in 180 genotypes of the Bornova (2015 and 2016) and Sanliurfa (2015 and 2016) are presented in Table . The Fe concentrations in both locations and years showed a normal distribution pattern. In Bornova 2015 and 2016, the maximum number of genotypes formed a group in the range of 60-70 mg/kg Fe concentration (Supplementary Fig. ). In Sanliurfa 2015 and 2016, the maximum number of genotypes formed a group in the range of 50-60 mg/kg Fe concentration (Supplementary Fig. ). The Zn concentration in most experimental cases showed a normal distribution pattern (the distribution of the data from Sanliurfa 2015 appears skewed distribution), and the maximum number of genotypes given on the bar ranged from 40 to 45 mg/kg (Supplementary Fig. ).

The genotype x environment interaction had significant effects on the Fe and Zn concentrations in the seeds (Supplementary Table ). There was a significant interaction between genotype and location for Fe concentration and a significant interaction between genotype and year for Zn concentration. In addition, there were significant genotype and location effects of Fe concentration and significant genotype and year effects of Zn concentration. The overall mean seed Fe and Zn concentrations were 62.43 and 42.94 mg/kg,

respectively (Table ). The mean seed Fe concentration was higher in Bornova compared to Sanliurfa in both years, but did not change within either location between the two years. There was an approximately threefold difference in the ranges of seed Fe concentrations for all experimental cases. The mean seed Zn concentration was very close in both locations for two years, but was slightly higher in 2016 than in 2015. Approximately twofold differences were observed in the range of Zn values across the four cases (Table ). In addition, the correlation between the Fe and Zn concentrations was slightly significantly negative (correlation coefficient, -0.15).

The mean seed Fe concentrations of the C. reticulatum and C. arietinum genotypes were 68.37 and 53.94 mg/kg, respectively (Table ). In all environments, the variability in the seed Fe concentration in the C. reticulatum genotypes was higher than that in the C. arietinum genotypes. The variability in the seed Fe concentration was threefold in C. reticulatum and twofold in C. arietinum (Table ). The mean seed Zn concentrations in the C. reticulatum and C. arietinum genotypes were 42.16 and 44.21 mg/kg, respectively. The mean variability of the seed Zn concentrations was very similar in the C. reticulatum and C. arietinum genotypes in all environments.

Genetic Diversity

Of the 431,274 SNPs generated, 410,406 were omitted from analysis after filtering due to having more than 10% missing data and less than 5% minor allele frequencies. The remaining 20,868 SNPs were aligned to the chickpea genome, of which 17,102 were aligned within the existing chromosomes and 3,766 were in the scaffolds. All 20,868 SNPs used in the analysis were distributed relatively evenly across the chromosomes and suitable for GWAS. The number of markers for each chromosome varied from 962 SNPs (Chr 8) to 2,823 (Chr 6) with an average of 2,318 SNPs per chromosome.

Population structure analyses were carried out to gain an insight into genetic diversity across the population. As a result of the population structure analysis, delta K showed a peak at two means, indicating that 180 genotypes clustered into two groups; red and green representing C. reticulatum and C. arietinum genotypes, respectively (Supplementary Fig. ). The spatial representation of the relative genetic distances between the chickpea genotypes in PCA also revealed two distinct groups (Supplementary Fig. and ), supporting the results of the STRUCTURE analysis. To provide an understanding of the genetic diversity in the chickpea population, an NJTree was constructed based on Nei’s (1972) genetic distance. The results of the analysis of this tree were also in agreement with those of the STRUCTURE and PCA analyses. Supplementary Fig. () shows that C. reticulatum and C. arietinum genotypes strictly diverged from each other. In addition to the population structure analysis with both C. reticulatum and C. arietinum (180 genotypes), the STRUCTURE and NJTree analyses were applied to the C. reticulatum and C. arietinum populations separately. When only C. reticulatum was included in the population structure analysis (Analysis B), two subpopulations were found (Supplementary Fig. ). When the same analysis was applied only to C. arietinum (Analysis C), three subpopulations were detected (Supplementary Fig. ). The population structure of the two Cicer populations varied according to the NJtree analysis given in Supplementary Fig. ().

Association Mapping for Fe Concentration

A total of 23 SNPs for all environments were found to be significantly associated with the seed Fe concentration after the FDR correction. In analysis C, no SNP was found to be associated with the seed Fe concentration. Based on analyses A and B, four SNPs were consistently observed (Table ), of which one (SNP204) was located on Chr 1, two (SNP8254 and SNP8255) on Chr 4, and one (SNP 9478) on Chr 5. Two SNPs were common in Bornova 2015 and Bornova 2016 and two in Sanliurfa 2015 and Sanliurfa 2016. In the analysis based on the general mean of the Fe concentrations, 4 SNPs were found in analyses A and B, four of them SNP4503, SNP4504, SNP4508, and SNP4512 were located on Chr 3. The manhattan plots of Fe concentration are given in Fig. ( and Supplementary Fig. ( and ). The Q-Q plots of Fe concentration are given in Supplementary Fig. ().

Fig. (1)

Manhattan plots of Fe concentration for Bornova 2015, Sanliurfa 2015, Bornova 2016, Sanliurfa 2016, and the general mean of concentrations based on −log 10 (P) values according to analysis A. The lower horizontal line represents FDR correction and the upper horizontal line represents Bonferroni’s correction. (A higher resolution / colour version of this figure is available in the electronic copy of the article).

Association Mapping for Zn Concentration

A total of 16 SNPs for all environments were found to be significantly associated with the seed Zn concentration after the FDR correction. In analysis C, no SNP marker was found to be associated with the seed Zn concentration. Four SNPs were common in analyses A and B, of which one SNP (SNP8284) was located on Chr 4 and three (SNP9528, SNP9529, and SNP10249) were located on Chr 5 (Table ). All of these SNPs were common in Bornova 2016 and Sanliurfa 2016 environments. In the analysis based on the general mean of the Zn concentrations, two SNPs were found to be common in analyses A and B (Table ). Of these SNPs, one (SNP2460) was located on Chr 2 and the other (SNP9097) on Chr 4. The Manhattan plots of Zn concentration are given in Fig. ( and Supplementary Fig. (). The Q-Q plots of Zn concentration are given in Supplementary Fig. ().

Fig. (2)

Manhattan plots of Zn concentration for Bornova 2015, Sanliurfa 2015, Bornova 2016, Sanliurfa 2016, and the general mean of concentrations based on −log 10 (P) values according to analysis A. Horizontal line represents FDR correction. (A higher resolution / colour version of this figure is available in the electronic copy of the article).

In the current analysis, in total four candidate genes, of which one was related to Fe concentration and three were related to Zn concentrations were identified (Table ).

DISCUSSION

Breeding for improved Fe and Zn concentrations in the chickpea can be an efficient and sustainable way of eliminating micronutrient deficiency from human diets. Association mapping studies have been performed on several staple crops to identify micronutrient accumulation in seeds [1, 28]. The current study identified associations between SNP markers and seed Fe and Zn concentrations in chickpea to contribute to the existing literature.

The results of this research revealed that genotypes and environmental factors had a significant effect on both Fe and Zn concentrations in chickpea seeds (Supplementary Table ). Similarly, Diapari et al. reported that there were important interactions between genotypes and the environment [1]. Other researchers also reported that the genotype x environment interaction had an impact on the grain Fe concentration, whereas the grain Zn concentration was associated with soil properties [29]. Supporting our results, Ray et al. showed that the mineral content of the chickpea, lentil, common bean and field pea was affected by the environment [30]. Recently, Khazaei et al. revealed the significant effect of the genotype x environment interaction on the seed Fe and Zn concentrations in the lentil [31]. The differences in the Fe and Zn concentrations we found in different environments (Table ) may be due to soil fertility, husbandry practices, and edaphic factors, such as moisture, pH, and other soil

chemical properties [32]. Furthermore, condition-specific (genetic) factors that control mineral concentration and bioavailability are other important factors that contribute to the differences in micronutrient concentrations between different environments [33]. A non-significant negative correlation was detected between the chickpea seed Fe and Zn concentrations unlike the positive correlation observed in the lentil [28a], barley [32], and pea [28b, 28c]. The most probable explanation for the improved Fe and Zn concentrations in the seeds is independent genetic selection.

The results of this study indicated substantial variation in the Fe and Zn concentrations between the genotypes (Table ). In breeding programs, this finding can be used to increase diversity for the improvement of Fe and Zn biofortification [28a]. The variation in the seed Fe concentration (38.85-100.41 mg/kg, seed Fe concentration) was higher between the genotypes used in this study than reported in other studies as [1] (37.6-68.7 mg/kg seed Fe concentration), and [3] (40.3-67.0 mg/kg Fe concentration). This may be because the genotypes used in the current research (both Cicer reticulatum and Cicer arietinum) had greater genetic variation compared to the other two studies (using only Cicer arietinum genotypes) because wild Cicer increases genetic diversity [34]. The variation in the seed Zn concentration (33.25-51.74 mg/kg Zn concentration) was very similar among the genotypes used in this study, as also found in other studies [1] (22.9-52.7 mg/kg Zn concentration) and [3] (27.9-48.5 mg/kg Zn concentration). Seed Zn concentration in the Cicer reticulatum genotypes varied between 35.85 to 49.7 mg/kg in this research, and 65.24 to 123.51 mg/kg in von Wettberg et al. [17]. In the current study, the level of variation in the Fe concentration was higher in the C. reticulatum genotypes than the C. arietinum genotypes (Table ). C. reticulatum is a wild progenitor of the C. arietinum, and collected from a small geographical region in Southeast Turkey that covered the full distribution of the species [35]. C. reticulatum has a higher level of genetic variation than C. arietinum, because, during the domestication of C. arietium, the species went through several genetic bottlenecks [34]. Supporting these findings, White and Broadley [13] reported that wild germplasm and landraces generally included high variability of macro and micronutrient content. It was also noted that wild species had higher levels of genetic diversity than their cultivated descendants but due to the process of domestication, the diversity at the genome and local levels decreased [35]. The variation level of the Zn concentration was similar between C. reticulatum and C. arietinum genotypes in the current study.

In this study, the Fe concentration was found to be higher in the C. reticulatum genotypes compared to the C. arietinum genotypes. Similarly, the Fe and Zn concentrations in the cereals of landraces and aromatic rice cultivars are reported to be four-fold higher than those in commercial rice cultivars [36]. It has also been shown that the seeds of wheat wild germplasm and landraces have three or four-fold higher concentrations of Fe and Zn compared to bread and durum wheat [37]. This might be because wild genotypes tend to have a lower yield and higher concentration of micronutrients, while cultivated genotypes have a higher yield but lower micronutrient concentration [28c]. Although we found a higher Fe concentration in C. reticulatum, the Zn concentrations were similar in the C. reticulatum and C. arietinum genotypes.

Several genotypes containing a high Fe concentration were detected, all belonging to C. reticulatum. The Bari3_r_079 genotype had the highest Fe concentration with an average of 100.41 mg/kg. Three other accessions, Bari1_r_066, Bari3_r_079, and Bari3_r_110, had higher Fe concentrations with an average of 95.5, 91.24, and 91.05 mg/kg, respectively. Several genotypes had a high Zn concentration; e.g., Gunas_e_063, Sarik_r_064, and Sirna_r_062 accessions belonging to C. reticulatum contained 48.6, 50.03, and 51.73 mg/kg Zn, respectively. Two C. arietinum accessions, namely TR51411 and TR42362 had a Zn concentration of 48.66 and 49.7 mg/kg, respectively, and both were of the Kabuli type.

In this study, the results of the population structure analysis indicated a separation between C. reticulatum and C. arietinum. Furthermore, the population structure analysis performed separately for C. reticulatum and C. arietinum provided an in-depth understanding of the two Cicer species. While two subpopulations were found within C. reticulatum in the current study, von Wettberg et al. [17] detected eight subpopulations.

In this study, a higher number of Fe and Zn associations were found compared to previous studies on the chickpea [1, 3], which might be due to the greater number of molecular markers used and the higher levels of genetic variability for these traits in the screened population, as noted in [28a]. Using both C. reticulatum and C. arietinum genotypes increases the genetic variability in the population. The strength of association mapping studies is dependent on the genetic variability of the screened population and the number of molecular markers [38]. In the current study, SNPs consistently associated with the seed Fe concentration were located on Chr 1, Chr 3, Chr 4, and Chr 5. Similarly, Diapari et al. [1] reported SNPs located on Chr 1, Chr 4 and Chr 5, and Upadhyaya et al. [3] found SNPs on Chr 1, Chr 3, Chr 4, Chr 5, and Chr 7. On the other hand, in this study, SNPs that were consistently associated with the seed Zn concentration were located on Chr 2, Chr 4, and Chr 5. In previous studies, SNPs were identified on Chr 1, Chr 4, and Chr 7 [1] and on Chr 3, Chr 4, Chr 5, and Chr 7 [3]. While several markers were found in analysis A and analysis B, no marker was detected in analysis C. The lack of significant marker-trait associations might be due to the lower levels of genetic variability for these traits in the Cicer arietinum genotypes. In addition, the small sample size of the experimental population adversely affects the results of association mapping studies [39].

A significant association between a trait and an SNP marker can arise from the SNP being within the gene or closely associated with the gene related to the seed Fe and Zn concentrations [1]. In the current study, four putative candidate genes were predicted from the sequences representing homology to the associated SNPs for the seed Fe and Zn concentrations (Table ). The marker SNP204 is associated with protein trichome birefringence-like 21, which may be involved in the specific O-acetylation of cell wall polymers [40]. The function of trichome birefringence is in primary cell wall synthesis in etiolated seedlings [41]. The marker SNP2460 is associated with deoxynucleoside triphosphate triphosphohydrolase SAMHD1 functioning in the defense response and regulation of DNA end resection at stalled replication forks [42]. The marker SNP9097 is associated with pentatricopeptide repeat-containing protein, which is involved in RNA editing event in chloroplast and required for the editing of a single site in rps12 transcript [43]. The marker SNP10249 is associated with serine/threonine-protein kinase ATM isoform. Serine/threonine protein kinase is necessary for activating checkpoint signaling upon genotoxic stresses, such as ionizing radiation or DNA replication stalling. It plays an important role in the perception and response to both developmentally programmed damage during meiosis and stress-induced damage in somatic cells [44]. It is critical for regulating the DNA damage response mechanism at dysfunctional telomeres and preventing the activation of functional telomeres [45].

When the gene annotation analysis was applied to the genes within 100K to their left and right on the SNP region, several candidate genes, two of which are very important, were identified (Supplementary Table ). The marker SNP2460 is associated with zinc transporter 3-like, which confers Zn uptake activities. The zinc transporter-3 gene is expressed in roots in response to Zn deficiency and transports Zn from the soil into the plant [46]. In addition, the marker SNP10249 is associated with DNA-binding One Zinc Finger protein DOF5.6. The DOF family transcription factors play a role in several fundamental processes in higher plants, such as seed maturation, germination, and responses to phytohormones [47]. DOF5.6 has specific functions in vascular tissue development and interfascicular cambium formation [48]. These two genes, which are related to the Zn concentration, should be closely linked in the markers given.

Considerable variation in the micronutrient concentration among staple crops allows improving micronutrient concentration in seeds in breeding programs. The development of genomic tools to identify the genetic factors associated with micronutrient concentrations assist in the selection process in these programs. The current research increased the genomic diversity by using both Cicer reticulatum and Cicer arietinum, and extended environmental range by conducting experiments in two different locations and years in Turkey. Gathering the wild relatives of crops is highly beneficial for breeding programs since they present the extent of adaptation in natural populations.

CONCLUSION

In conclusion, this study provided an insight into the genetic basis of variability in the Fe and Zn concentrations of seeds of a chickpea population containing a diverse set of genotypes. This research identified the population structure and genetic variation in the selected population and determined the SNP markers significantly associated with the seed Fe and Zn concentrations. The identified markers can be used in marker-assisted selection to improve the Fe and Zn concentrations in chickpea seeds with good qualities. Furthermore, the genotypes containing high Fe and Zn concentrations can be used as parents to obtain Fe and Zn-biofortified chickpea genotypes.

Table 1

Minimum, maximum and mean ± standard deviation (SD) values of Fe and Zn concentrations in Bornova (2015 and 2016) and Sanliurfa (2015 and 2016), and all cases combined.

-	Fe Conc (mg/kg)			Zinc (mg/kg)
Experimental Case	All Genotypes	C. reticulatum	C. arietinum	All Genotypes	C. reticulatum	C. arietinum
Bornova 2015	-	-	-	-	-	-
Mean + SD	64.40 ± 13.27	70.82 ± 12.46	55.33 ± 8.12	41.66 ± 4.46	41.08 ± 4.27	42.51 ± 4.62
Range	40.00-112.70	46.00-112.70	40.0-77.80	29.77-41.66	31.58-49.51	29.77-49.99
Bornova 2016	-	-	-	-	-	-
Mean + SD	64.35 ± 13.25	70.81 ± 12.41	55.2 ± 7.98	43.90 ± 5.06	43.55 ± 5.80	44.41 ± 3.69
Range	40.00-111.72	46.02-111.72	40.0-77.8	30.58-61.90	30.58-61.9	35.25-52.24
Sanliurfa 2015	-	-	-	-	-	-
Mean + SD	60.72 ± 13.35	66.60 ± 14.00	52.66 ± 6.69	42.29 ± 5.13	40.08 ± 4.71	45.54 ± 3.86
Range	37.70-133.90	41.45-133.9	37.7-65.5	21.06-49.96	25.76-49.58	21.06-49.96
Sanliurfa 2016	-	-	-	-	-	-
Mean + SD	60.25 ± 11.90	65.5 ± 11.95	52.69 ± 6.69	43.90 ± 5.07	43.56 ± 5.82	44.41 ± 3.69
Range	37.70-100.80	41.50-100.80	37.7-65.5	30.58-62.00	30.58-62.00	35.25-52.24
General mean of concentrations	-	-	-	-	-	-
Mean + SD	62.43 ± 11.58	68.37 ± 10.44	53.94 ± 7.00	42.94 ± 3.32	42.16 ± 3.35	44.21 ± 2.84
Range	38.85-100.41	50.76-100.41	38.85-69.65	33.25-51.74	33.25-51.74	35.85-49.7

Table 2

The list of SNP markers significantly associated with the Fe concentration according to both analyses A including both C. arietinum and C. reticulatum and analysis B consisting of only C. reticulatum genotypes.

Analysis	Marker	Location/year	Chromosome	Position (bp)	P	-LOG(P)	R2
Analysis A	SNP204	Sanliurfa 2015	1	3212701	1.045E-06	5.9810084	0.14348
		Sanliurfa 2016			3.76E-04	3.4248815	0.07337
Analysis B		Sanliurfa 2015			1.55E-05	4.8107062	0.19583
Analysis A	SNP8254	Bornova 2015	4	41702329	0.0001704	3.768454	0.08225
		Bornova 2016			0.0001666	3.7783771	0.08249
Analysis B		Bornova 2015			0.0002548	3.5937665	0.13664
		Bornova 2016			0.0002548	3.5937665	0.13664
Analysis A	SNP8255	Bornova 2015	4	41702981	0.0001704	3.768454	0.08225
		Bornova 2016			0.0001666	3.7783771	0.08249
Analysis B		Bornova 2015			0.0002548	3.5937665	0.13664
		Bornova 2016			0.0002548	3.5937665	0.13664
Analysis A	SNP9478	Sanliurfa 2015	5	27634952	7.494E-06	5.1253153	0.11942
Analysis B		Sanliurfa 2015			7.14E-06	5.146466	0.21287
		Sanliurfa 2016			3.82E-04	3.4179253	0.12828

Table 3

The list of SNP markers significantly associated with the Zn concentration according to both analyses A including both C. arietinum and C. reticulatum and analysis B consisting of only C. reticulatum genotypes.

Analysis	Marker	Location/year	Chromosome	Position (bp)	P	-LOG(P)	R2
Analysis B	SNP8284	Bornova 2016	4	42385136	0.0002372	3.6248304	0.13781
		Sanliurfa 2016			0.0002352	3.6286365	0.13799
Analysis B	SNP9528	Bornova 2016	5	28469147	0.0003239	3.4895891	0.13148
		Sanliurfa 2016			0.0003257	3.4871956	0.13137
Analysis B	SNP9529	Bornova 2016	5	28469199	0.0003239	3.4895891	0.13148
		Sanliurfa 2016			0.0003257	3.4871956	0.13137
Analysis A	SNP10249	Bornova 2016	5	39588845	0.000194	3.7121087	0.08117
		Sanliurfa 2016			1.96E-04	3.7085423	0.08108

Table 4

The list of SNP markers significantly associated with seed Fe and Zn concentrations according to both analyses A including both C. arietinum and C. reticulatum and analysis B consisting of only C. reticulatum genotypes based on the general mean of concentrations.

Trait	Marker	Analysis	Chromosome	Position (bp)	P	-LOG(P)	R2
Zn	SNP2460	Analysis A	2	2347203	0.000375	3.425818	0.0728
		Analysis B			6.52E-05	4.185626	0.16395
	SNP9097	Analysis A	4	37724523	0.000191	3.719058	0.07549
		Analysis B			1.58E-05	4.80148	0.07173
Fe	SNP4503	Analysis A	3	13716969	0.000249	3.604115	0.07754
		Analysis B			0.000267	3.57305	0.13562
	SNP4504	Analysis A	3	13717369	0.000249	3.604115	0.07754
		Analysis B			0.000267	3.57305	0.13562
	SNP4508	Analysis A	3	13777857	0.000249	3.604115	0.07754
		Analysis B			0.000267	3.57305	0.13562
	SNP4512	Analysis A	3	13814197	0.000249	3.604115	0.07754
		Analysis B			0.000267	3.57305	0.13562

Table 5

Results of the BLAST analysis on the NCBI databases indicating candidate genes with putative functions related to the SNP markers associated with Fe and Zn concentrations.

Trait	SNP ID	Chromosome	Enzymes and Proteins Related to Putative Candidate Genes	Species
Fe concentration	SNP204	1	protein trichome birefringence-like 21	Cicer arietinum
Zn concentration	SNP2460	2	deoxynucleoside triphosphate triphosphohydrolase SAMHD1 homolog	Cicer arietinum
	SNP9097	4	pentatricopeptide repeat-containing protein At1g22960, mitochondrial	Cicer arietinum
	SNP10249	5	serine/threonine-protein kinase ATM	Cicer arietinum