In the face of increasing climate variability and the growing global demand for sustainable protein sources, legume crops are gaining renewed attention for their role in food security and environmental resilience. Among them, chickpea (Cicer arietinum L.) has emerged as a valuable crop, especially in southern Europe, due to its favorable nutritional profile, drought tolerance and adaptability to low-input systems (Ilić et al., 2024). Chickpea is an important pulse crop cultivated and consumed worldwide, ranking third in global legume production with an annual output exceeding 11.5 million tons (Merga and Haji, 2019). It is a rich source of carbohydrates, protein and essential vitamins, with protein content ranging from 18.2% to 26.7% and its protein quality is often considered superior to that of many other pulses (Sodavadiya et al., 2023; Sahu et al., 2022). Chickpea is predominantly grown under rainfed conditions, relying on residual soil moisture with little precipitation during the growing season. As the season progresses, soil moisture declines, leading to terminal drought, which is one of the most critical abiotic stresses affecting about two-thirds of the global chickpea area. This shortened cropping window necessitates the selection of genotypes with appropriate growth duration and improved mechanisms for drought avoidance (dehydration postponement) or drought tolerance (dehydration resistance), to maintain stable yields under stress-prone conditions (Gaur et al., 2019; Karim et al., 2021).
       
Phenological traits, particularly days to flowering, are key indicators of chickpea adaptation and yield stability in such environments (Mallikarjuna et al., 2019). However, chickpea phenology is more complex than in many other crops, influenced by local factors such as temperature, photoperiod and soil moisture. These interact with genotype-specific responses, contributing to variability in development and productivity across growing conditions (Gimenez et al., 2025).
       
For local breeding programs targeting yield improvement and stress resilience, the identification of genetically diverse and agronomically promising chickpea genotypes is essential. Morphological and phenological traits remain important indicators of adaptability and productivity. At the same time, genotype × environment interaction (GEI) represents a significant source of variation in yield across locations, underscoring the need for multi-environment trials to identify genotypes with both high performance and broad adaptability (Farshadfar et al., 2013; Arshad et al., 2003). In addition to field-based evaluation, molecular markers, such as SSRs and SNPs, provide complementary insights into genetic diversity. They enable the detection of duplicates, assist in designing efficient sampling strategies and support more precise selection in breeding, thereby accelerating cultivar development (Mir et al., 2022; Abbo et al., 2005).
       
Our work reports the results of testing a collection of chickpea accessions in two agroecological zones of Serbia, a country where chickpea has not yet become a major crop. This study was conducted in two complementary phases during a single growing season. The first phase focused on the characterization of chickpea genotypes to assess morphological and phenological variability. In the second phase, a subset was evaluated under three distinct agronomic settings: two contrasting input systems at one location and a third environment at a separate location, to investigate yield performance, stability and adaptability. In addition, molecular profiling of the selected genotypes was performed to strengthen the phenotypic assessment and contribute to a more comprehensive understanding of genotype differentiation. To the best of our knowledge, this is the first study in Serbia to combine phenotypic and molecular characterization of chickpea genotypes under contrasting agronomic conditions. Given the increasing interest in chickpea cultivation across Southeast Europe, the results of this study may contribute to broader efforts aimed at improving drought resilience and productivity in similar agroecological settings, while also emphasizing the potential for expanding chickpea production in the Balkan region, where it remains underexploited, despite favorable agroecological conditions.

Field observations
 
The study was conducted in two phases at the Institute of Field and Vegetable Crops, Novi Sad and Institute for Vegetable Crops, Smederevska Palanka (Serbia). In the first phase, 16 chickpea genotypes of diverse origin were evaluated under uniform field conditions at a single location (Rimski Šančevi, hereafter referred to as RŠ, Republic of Serbia: 45°20'N, 19°51'E) to assess their morphological and phenological diversity (Table 1). The number of genotypes was constrained by the limited availability of chickpea germplasm, as this crop is not traditionally grown in Serbia and is only recently being introduced into local research and potential production systems. The studied genotypes represent inherited germplasm accessions maintained in a collection; they are neither released cultivars nor advanced breeding lines and detailed pedigree information is not available due to their historical origin. Sixteen genotypes were chosen as they represent the full set of available chickpea germplasm accessible for field evaluation in Serbia at the time of the study, with selection criteria focused on maximizing phenotypic diversity. The fields at Rimski Šančevi were sown on chernozem soil (WRB, 2024). The observations were carried out in 2024 and involved genotypes belonging to the Kabuli or intermediate seed type, characterized by light-colored seed coats and various sizes and shapes. Randomized complete block design was used, with three replications per genotype. The sowing density was 50 plants per square meter and each plot measured 5 m². Standard phenotypic, assessment including plant height (PH), phenological stages (beginning of flowering (BF), 50% flowering (X50FL), 100% flowering (X100FL), end of flowering (EFL) and maturity (MAT)), pod and seed traits (number of pods per plant (PPP), seeds per pod (SPP) and 100 seed weight (SW)), protein content in seeds (P%) and seed yield (SY). Protein content (P%) in chickpea seeds was determined using Near-Infrared Reflectance Spectroscopy (NIRS) with an Antaris II spectrometer (Thermo Scientific, USA). Seed samples were ground prior to analysis and protein concentration was expressed as a percentage of dry matter based on established calibration models.


       
Based on the results from the trial during the preceding year (2023), five genotypes with distinct morphological profiles and promising agronomic traits were selected for further evaluation in the second phase. This included three field trials carried out during the same growing season but under different agronomic management systems. Main yield components, i.e., number of pods per plant (PPP), number of seeds per plant (SPP) and seed yield per plant (SY), were assessed as principal indicators of agronomic performance. The first two trials were performed at the same location (RŠ), under two contrasting management systems: one with full conventional agrotechnics (RS1) and the other under a low-input approach (RS2). Neither trial was irrigated, relying solely on rainfall and the chickpea’s inherent drought tolerance. The full agrotechnical system included standard fertilization and chemical protection, whereas the low-input system involved minimal external inputs and reduced soil disturbance. The third experiment was conducted at Smederevska Palanka (SP), Republic of Serbia (44°22WRB, 2024N, 20°57WRB, 2024E; 103 m a.s.l.), on a flat site surrounded by the rivers Kubršnica and Jasenica. The soil is classified as alluvial smonitza, formed through river sediment deposition and long-term anthropogenic influence, with relatively high humus content, abundant plant-available potassium, good phosphorus availability and a neutral pH reaction (Dugalić and Gajić, 2012; Ugrinović, 2015).
       
Climatic data for all experimental sites during the 2024 growing season were collected from the local meteorological station near the trial site, the Annual Bulletin for Serbia for 2024 (RHMZ, 2025) and are presented in Fig 1 to support the interpretation of genotype performance under variable environments.


       
Statistical analyses were conducted using R software (R Core Team, 2024) and XLSTAT (Addinsoft, 2024). Descriptive statistics and Pearson correlation coefficients were calculated to evaluate variability and trait associations. Principal component analysis (PCA) was used to determine key traits contributing to genotype differentiation, while heatmap visualization was employed to explore phenotypic diversity patterns. Genotype × environment (G×E) interaction was analyzed to assess the yield stability and adaptability of the selected genotypes across different agronomic conditions. In R, data processing and visualization were performed using the tidyverse (Wickham et al., 2019) and ggplot2 (Wickham, 2016) packages. Correlation analysis and descriptive statistics were performed using psych (Revelle, 2023), PCA was conducted and visualized using FactoMineR and factoextra (Lê et al., 2008; Kassambara and Mundt, 2020) and heatmaps were generated with pheatmap (Kolde, 2019). In the second phase of the experiment, the agricolae package (de Mendiburu, 2021) was used for agronomic comparisons and GGE biplot analysis was carried out with GGEBiplotGUI (Frutos, 2014). Two-way analysis of variance (ANOVA) was performed in the program InfoStat ver. 2020 (Di Rienzo et al., 2020), to examine differences between the means of the chickpea genotypes for the studied traits and to identify main effects (G, E, G×E) contributing to overall variation.
 
DNA polymorphism analysis
 
DNA was extracted from freshly collected leaves using modified CTAB method (Torres et al., 1993). Six plants were randomly sampled from each line. Final air-dried DNA pellets were dissolved in 50 μl deionized water. The concentration of DNA in samples was evaluated using an Evolution 201 spectrophotometer (Thermo Scientific) and the manufacturer’s software. For the subsequent PCR, DNA was dissolved to a concentration of 100 ng/μl.
       
To overcome the potential heterogeneity of the lines, bulk samples were prepared by pooling 1.5 μg aliquots of DNA from each individual plant and these samples were used for subsequent PCR. To reveal polymorphism, we used iPBS retrotransposon-targeted markers (Kalendar et al., 2010). As retrotransposons are abundant and ubiquitous in genomes, iPBS markers allow rapid detection of numerous polymorphic loci with a low taxon specificity, i.e. the same set of markers can be used for numerous plant species. In particular, 12-18-merous primers 2076, 2079, 2222, 2230, 2373, 2378, 2383 and 2415 (Kalendar et al., 2010) were synthesized by Microsynth AG and used in our work.
       
Each PCR microtube contained 1 μl (1 U) of Taq DNA polymerase, 2 μl of Taq 10× buffer with KCl, 1.6 μl of 25 mM MgCl₂ (all reagents produced by Thermo Scientific), 300 ng of DNA, 2 μl of primer (20 ng/μl), 1.25 μl of bovine serum albumin (Sigma-Aldrich) and nuclease-free water to 20 μl.
       
PCR was conducted on a 2720 Thermal Cycler (Applied Biosystems). The following reaction program was used: 3 min at 95°C; 40 cycles of 30 s at 94°C, 30 s at Tm (appropriate for each primer), 2 min at 72°C; a final extension of 5 min at 72°C. The PCR products were separated in 2% agarose (SeaKem LE Agarose, Lonza) gel and visualized and photographed in a TFP-M/WL transilluminator (Vilber Lourmat). The photographs were analyzed using GelAnalyzer 23.1.1 (designed by Istvan Lazar Jr., PhD and Istvan Lazar Sr., PhD, CSc). After the automatic detection of lanes and quantification of fragment sizes, these were additionally edited manually to exclude artifacts. Fragments differing in length by less than 5 bp were considered identical. The presence or absence of each fragment was encoded in a binary matrix, which was then analyzed using PAST 4.03 (Hammer, 2001).

Morphological traits and correlations studies
 
The diversity of 16 chickpea genotypes based on the studied phenotypic traits is presented in Table 2. For all studied accessions, plant height ranged from 38.6 cm (SRBCIC 4) to 52.5 cm (SRBCIC14), with a mean value of 45.9 cm. Most of the accessions began to flower around 76 days from sowing, while they completed flowering after 106 days. On average, full maturity was reached in about three months after the initial sowing date. Low CV values for phenological traits among the genotypes indicate their uniformity in relation to the development of different phenological stages. In general, accessions produced 38.6 pods per plant (range 18 to 76.3) and 1.2 seeds per pod (range 1 to 2) with substantial variation (CV = 36.6% and 33.9%, respectively). Seed yield per plant varied from 5.1 g (SRBCIC4) to 29.2 g (SRBCIC3), offering large variability (CV = 53.3%). Studied accessions displayed an average 100 seed weight of 27.3 g (range from 18.9 g (SRBCIC1) to 45.3 g (SRBCIC9)), while protein content, as determined by NIRS, spanned between 21.2% (SRBCIC21) to 25.5% (SRBCIC3), with low variability (CV = 5.7%).


       
Strong positive correlations were observed between seed yield per plant and pods per plant (r = 0.88, p<0.01) and between end of flowering and maturity (r = 0.86, p<0.01) (Fig 2). Protein content correlated with seeds per pod (r = 0.72, p<0.01) and phenological traits, while seed weight correlated negatively with seeds per pod (r = -0.50, p<0.05). Weak but reliable correlations were found between pods per plant and plant height (r = 0.35, p<0.01) and between pods per plant and seeds per pod (r = 0.39, p<0.01) (Fig 2). In Serbia, the main determinant of yield is the number of pods per plant, which varied considerably across accessions. It agrees with the results of Yücel et al. (2006), whereas in some trials this trait was found correlating much weaker with plant yield (Ali and Ahsan, 2012).


       
The structure of phenotypic variation among the studied chickpea genotypes were assessed using principal component analysis (PCA) and a heatmap. The first two principal components accounted for 34.76% and 21.1% of the total variation, respectively (Fig 3), reflecting a considerable phenotypic diversity consistent with previous reports (Jain et al., 2023; Sellami, 2021). Several genotypes clustered together based on similar phenotypic profiles. Genotypes SRBCIC9, SRBCIC12, SRBCIC14 and SRBCIC15, located in the upper-left quadrant of the biplot, exhibited the longest stems, largest seed size and highest yield per plant, in agreement with studies linking seed size and stem length with yield (Jain et al., 2023). SRBCIC5, SRBCIC6 and SRBCIC17 clustered near the phenological traits vectors, showing no significant differences in vegetation length or stage duration, which aligns with observations in other chickpea collections (Sellami et al., 2021). Genotypes SRBCIC1, SRBCIC2 and SRBCIC3 were distinguished by either the highest number of seeds per pod or the highest protein content, but the lowest 100-seed weight, reflecting variability in agronomic and nutritional traits reported in previous studies (Jain et al., 2023). The remaining six genotypes (SRBCIC4, SRBCIC7, SRBCIC8, SRBCIC10, SRBCIC11, SRBCIC21) did not form distinct clusters and showed intermediate values for most traits, suggesting balanced phenotypic characteristics suitable for further evaluation. A heatmap (Fig 4) based on standardized trait values and seed shape was used to explore the relationships between chickpea genotypes with respect to seed shape, the most efficient qualitative trait to distinguish the studied lines. The dendrogram revealed two main clusters, each comprising genotypes with similar phenotypic profiles. The first cluster grouped accessions predominantly characterized by the ‘pea-shaped’ seed type, showing lower values for traits such as number of seeds per pod, seed weight and seed yield per plant. The second cluster included mostly ‘angular’ and ‘owl head’ seed types, which tended to have higher values for traits related to seed size and phenological characteristics (e.g., MAT, EFL and X100FL). The separation of genotypes according to seed morphology indicates a potential association with specific agronomic traits. Seed shape is important because it differentiates genotypes and correlates with key traits such as seed size and yield. In our study, ‘owl head’ seeds showed higher 100-seed weight and yield than ‘pea-shaped’ and ‘angular’ seeds, suggesting that seed morphology may serve as a practical marker for selecting superior genotypes.




 
DNA polymorphism
 
In total, 370 bands were detected in electrophoretic gels, with only eight of them (2.2%) being monomorphic, i.e., present in all 16 lines. The binary matrix of bands distribution in 16 genotypes is available in Supplement 2. Cluster analysis and principal component (PC) analysis both revealed several genotype groups (Fig 5A, B). However, this grouping found relatively low statistical support. The first two PCs together explained only 22.5% of the total variance. In cluster analysis, only the distal nodes were supported by bootstrap values exceeding 50 (Fig 5B). The following genotypes were grouped in both analyses; SRBCIC14 (angular seeds) + SRBCIC21 (pea-shaped); SRBCIC58 (owl head and pea-shaped); SRBCIC9 + SRBCIC14 + SRBCIC15 + SRBCIC17 (owl head); SRBCIC11 + SRBCIC12 (owl head, although this grouping was not supported in cluster analysis). SRBCIC10 occupied a somewhat isolated position among the studied cultivars in PCA (Fig 5A). With a single exception, it is primarily the seed type that matches the iPBS-based grouping of genotypes. The relatively low support for this grouping may be due to the fact that, during the maintenance of germplasm collection, different accessions are grown in open field plots. Although chickpea is preferentially self-pollinated, its flowers are nectariferous and visited by numerous insect species, some of which may contribute to cross-pollination (Latif et al., 2019). For chickpea seed production, it is recommended to sow various accessions separated by 5-10 m isolation distance (Gaur et al., 2010). Mechanical harvesting is also associated with the risk of cross-cultivar seed contamination. The only phenotypic markers that differentiate the studied germplasm are seed traits, such as colour and shape, that can be used to evaluate phenotypic uniformity. Taken together, these factors may explain the relatively low support (bootstrap often not exceeding 50) of iPBS-based grouping, which, however, agrees with seed phenotypes. This offers prospects for further differentiation of the examined chickpea collection by applying approaches such as single-seed descent, which would reduce genetic heterogeneity and allow for more precise genotyping. The resulting sublines may prove even more beneficial for breeding and more genetically uniform than the original accessions. Interestingly, the iPBS-based grouping of accessions at least partly fits their similarity found during the PCA of agriculturally valuable characteristics (Fig 3). For example, SRBCIC1-3 formed a cluster in both analyses. The most reliable phenotypic characteristics that match this grouping are seed type, such as angular in the case of SRBCIC1-4, owl head type in SRBCIC14 and SRBCIC15 etc. This suggests that seed features are the most efficient phenotypic markers that differentiate the accessions and also partially underlie the breeding value of the studied lines. Taken together, phenotypic, biochemical (protein) and molecular data confirm genotypic diversity and highlight potential breeding value for chickpea improvement under variable environmental conditions in Serbia.


 
Agronomical performance of selected chickpea accessions
 
Five accessions, i.e., SRBCIC1, SRBCIC2, SRBCIC14, SRBCIC17 and SRBCIC21, were selected for more detailed evaluation of their agronomical performance. Three main yield components, number of pods per plant (PPP), number of seeds per plant (SSP) and seed yield per plant (SY), were assessed at three different localities/agronomical conditions (Fig 6). The only statistically reliable difference was found for the number of pods per plant of SRBCIC1 between Rimski Šančevi 1 (RS1) and Smeredevska Palanka (SP) (Dunn test, p<0.05). In general, agronomic conditions at the location SP were the most suitable for the production of the studied chickpea accessions. Climatic data indicated that SP had slightly lower mean temperatures but substantially higher rainfall during March to June (SP: 35.5-71.9 mm vs. RS: 5.3-26.3 mm), which likely contributed to higher yields. The largest seed yield, number of pods and seeds per plant were observed at this location for selected genotypes. Only exemptions were accessions SRBCIC1 and SRBCIC2, which produced slightly more seeds per plant at the location Rimski Šančevi 1 (RS1). On the contrary, all accessions displayed the lowest yield at the low input location, RŠ2, with a two- and three-fold decrease in yield compared to locations RS1 and SP, respectively. Agronomical evaluation of the selected chickpea accessions revealed substantial variation in key yield components across locations. ANOVA indicated that environment (location, E) was the main factor driving trait variation, accounting for 87.2%, 49.1% and 67.3% of the variance for pods per plant, seeds per plant and seed yield, respectively (Table 3). Genotype (G) was significant only for seeds per plant (31.3%), while G×E interactions were not significant, suggesting environmental differences largely influenced performance (Table 3). Within this limited set of genotypes, small-seeded lines (SRBCIC1 and SRBCIC2) showed relatively stable yields across locations, whereas heavier-seeded accessions (SRBCIC14, SRBCIC17) were more responsive to favorable conditions. These findings should be interpreted cautiously due to the small number of genotypes evaluated. Conversely, heavier-seeded genotypes appear more responsive to favorable conditions and achieve higher yields when rainfall and management are adequate. Overall, these findings indicate that local climate, especially precipitation, strongly affects chickpea productivity. For regions with lower and more erratic rainfall, low-input management of stable, small-seeded accessions may be a reliable strategy, whereas heavier-seeded genotypes can be exploited under conditions with higher water availability.

The study demonstrates that environmental conditions, particularly precipitation, are the main drivers of variation in chickpea yield components. Small-seeded accessions showed high stability across variable environments and are therefore suitable for low-input or drought-prone regions. In contrast, heavier-seeded genotypes achieved higher yields under favorable conditions, indicating their responsiveness to improved agronomic practices. Morphological variation, including plant height, seed size and phenological traits, was substantial among the studied accessions and PCA and heatmap analyses revealed clear phenotypic groupings that corresponded to seed shape and other agronomically relevant traits. Genotypic analysis confirmed that seed traits are reliable phenotypic markers, although some heterogeneity exists due to open-field maintenance and potential cross-pollination. Overall, these findings emphasize the importance of integrating phenotypic and environmental data when selecting chickpea accessions for stable and productive cultivation in different climatic conditions.

The present study was supported by the European Cooperative Programme for Plant Genetic Resources, ECPGR: “Exploring grain legumes diversity for sustainable European agri-food systems (ExploDiv)’’ 2023-2025 and the Ministry of Science, Technological Development and Innovation of the Republic of Serbia, grant number: 451-03-136/2025-03/200032. 
 
Disclaimers
 
The views and conclusions expressed in this article are solely those of the authors and do not necessarily represent the views of their affiliated institutions. The authors are responsible for the accuracy and completeness of the information provided, but do not accept any liability for any direct or indirect losses resulting from the use of this content.

The authors declare that there are no conflicts of interest regarding the publication of this article. No funding or sponsorship influenced the design of the study, data collection, analysis, decision to publish, or preparation of the manuscript.