An association panel consisting of a collection of 284 diverse germplasms obtained and extracted from ICRISAT (Upadhyaya, 2003) along with four checks, BG 372, BG 3022, BG 547, and BG 1105, were used for GWAS analysis on nodulation traits. The experimental trials for the association panel were conducted at four environmental locations, namely, IARI, New Delhi (location 1; 28°38'24.0252” N latitude, 77°10'26.328” E longitude and 228.6 m AMSL) having sandy clay loam soils; Sam Higginbottom University of Agriculture, Technology and Sciences (SHUATS), Naini, Prayagraj (location 2; 25°24'41.27” N latitude, 81°51'3.42” E longitude and 98 m AMSL) with clay loam to sandy loam soil; KVK, Vaishali, Dr. Rajendra Prasad Central Agricultural University (RPCAU), Samastipur (location 3; 25°86'29.679”N latitude, “85°78'10.263”E longitude to accurately reflect the geographical location, and 52 m AMSL) with sandy loam soil; and IARI Regional station, Pusa, Bihar (location 4; 25°54'56.16” latitude, 85°40'24.9564” longitude, and 52 m AMSL) with alluvial soils during the year 2020–2021. Experimental trials consisted of 284 germplasm lines. The field in all the locations was divided into four blocks with four checks in each block. Each check was replicated three times in all the blocks. The experimental design was an augmented random block design with a row length of 5 m, a row spacing of 60 cm, and a plant-to-plant distance of 10 cm (Supplementary Table 1). Observations were taken for the following six traits associated with nodulation: number of nodules per plant (NON/plant), nodule fresh weight per plant [NFW (mg)/plant], number of pods per plant (NOP/plant), number of seeds per plant (NOS/plant), seed yield [SY (g)/plant], and 100 seed weight [100 SW (g)/plant]. Because NON is the ability of the nodulated plant to initiate symbiotic associations with rhizobia, a higher nodule count generally indicates enhanced opportunities for nitrogen fixation. Furthermore, efficiency also depends on nodule activity; therefore, observations on NFW were recorded. This trait is directly associated with metabolic activity and nitrogen fixation capacity. Larger and healthier nodules indicate effective rhizobial colonization and improved nitrogen assimilation.

The phenotyping for nodulation traits was done at 60 days after sowing (DAS), where plants were at 50% flowering stage. No external inoculation of Rhizobium was done to capture the natural variation of nodulation and associated traits. Five randomly selected plants from each genotype were uprooted with the adhered soil mass using a hand hoe by digging 20 cm or even deeper into the soil based on the plant growth and root length observed. Particular care was taken not to disturb the root nodule system during the sampling and removal of adhered soil. Root and shoot systems were separated. Roots with intact nodules were washed and counted for NON and stored in butter paper bags to further take NFW. NOP, NOS, SY, and 100 SW were taken by randomly selecting five plants for each genotype. Yield and 100 SW were measured in grams. The “corrplot” (Wei and Simko, 2021). R package was used to estimate Pearson’s correlation among the measured traits, while the “Factoextra” R package was used to undertake the principal component analysis (PCA) for the filtered data. The frequency distribution plots were generated using the “ggplot2” package in the R environment. Best linear unbiased predictors (BLUPs) were generated using Plant Breeding Tools Version 1.4 (PB Tools, 2014). Bartolome et al. (2014).

The genotypic/SNP marker data for the 284 accessions of the association panel were obtained in the HapMap format from the database of the Centre of Excellence in Genomics and Systems Biology (CEGSB), ICRISAT, Patencheru, Telangana (https://cegresources.icrisat.org/cicerseq/; Varshney et al., 2021). Out of a total of 2,470,880 raw SNPs, 355,546 SNPs were obtained for the association panel (284) and further analyzed. The following parameters were used for SNP filtering: missing data with less than 20%, minor allele frequency (MAF) cutoff at 0.05, and an additional filter for the rate of heterozygosity of less than 10% and after calculating the threshold value of the working set of SNPs (355,546) from the Bonferroni correction as 6.85. We utilized the vcf2gwas tool (Vogt et al., 2021) for filtering SNPs.

Identification of marker–trait associations (MTAs) from genome-wide markers was conducted following GWAS analysis. The generated genotypic data were integrated with multi-locational phenotypic data recorded for the traits NON, NFW, NOP, NOS, SY, and 100 SW. During the phenotypic analysis of multi-location augmented trial design, the checks were considered as fixed effects while the remaining factors (location, block, and new entries/genotypes) were treated as random effects. The phenotypic data were first converted into BLUPs using statistical methods [analysis of variance (ANOVA)] performed using SAS v9.4 mixed procedure Cody, R., (2018). The population structure was assessed using neighbor-joining phylogenetic tree constructed through TASSEL software and visualized through ITOL software and PCA. For PCA, we considered the first three principal components as covariates in the Genome Association and Prediction Integrated Tool (GAPIT) in R software. The extent of linkage disequilibrium (LD) between the SNP markers was analyzed by calculating the LD (r²) values in TASSEL. The only r² values with p < 0.05 within each chromosome were considered for LD decay analysis. LD decay plots were drawn by using r² and physical distances (bp) measured by using the script following R version 4.1.1. Chromosome-wise SNP distribution with the “SNP density plot” was generated using the web tool SR-Plots (https://www.bioinformatics.com.cn/en). To identify significant MTAs, GWAS was carried out utilizing the GAPIT3 package Wang and Zhang, 2021) by employing Bayesian-information and Linkage-disequilibrium Iteratively Nested Keyway (BLINK) and fixed and random model circulating probability unification (FarmCPU) models (Liu et al., 2016). The Manhattan and quantile–quantile (Q–Q) plots were generated from qqman version 0.1.8 (Turner, 2018) using the rMVP package (0.99.17; https://github.com/xiaolei-lab/rMVP). The BLASTn search for the chickpea genome of GCA 000331145.1 on the NCBI database (https://blast.ncbi.nlm.nih.gov) was performed to investigate the genomic locations of the significant SNPs in order to find their corresponding genes.

ANOVA, frequency distribution, and correlation for the phenotypic data of the association panel collected under all the environments were statistically analyzed and the results are presented below (Table 1).

Thus, notable statistically significant variations were observed for all the studied traits (p < 0.01). Frequency distribution depicted that all the studied traits were distributed normally in the population (Figure 1).

Combined descriptive statistics like minimum, maximum, mean, standard deviation, and coefficient of variation for all the studied traits for all environments were calculated. Based on the descriptive statistics, NON ranged from 5.53 to 80, NFW ranged from 141.53 to 4,646.0, NOP ranged from 7.12 to 123.68, NOS ranged from 7.60 to 113.40, SY ranged from 16.0 to 73.71, and 100 SW ranged from 7.12 to 123.68. Pearson correlation coefficient among different traits indicated that there was a significant positive correlation between NON and NFW with the yield of the studied genotypes along with positive correlation among all the studied traits (Figure 2).

Based on mean of genotypes across the locations, the highly nodulating genotypes found in the association panel are presented (Figure 3). The premise was to identify a set of accessions (ICC 7390, ICC 15, ICC 8348, and ICC 2474) that can be incorporated in breeding programs for enhancing nodulation traits without adversely affecting other agronomic traits screened in different environments.

Population structure was assessed by admixture (Figure 4A), PCA (Figure 4B), and phylogenetic tree (Figure 4C). The admixture plot provides valuable insights into the genetic diversity and structure of the study population. The ancestral populations (K) on the X axis vs. cross-validation on the Y axis were used to determine the optimal number of genetic clusters or k in a data set when performing admixture analysis, and we found good CV at k = 16. PCA revealed that the population could be divided into three groups with significant genetic variations among the genotypes studied (Figures 4A, B). The presence of the three sub-populations was further confirmed from the neighbor-joining diversity tree (Figure 4C).

The SNP density per chromosome ranged from 1,102 to 8,802 SNPs. The LD across the genome was estimated from the HapMap file containing 355,546 SNPs, and the average LD estimated across the genome was 635.9 kb (Figures 4D, E, respectively). Chromosome 4 contained a greater number of SNPs and chromosome 8 had the lowest number of SNPs. In our study, a huge number of high-quality SNPs were found in GWAS that increased the probability of detecting all the possible causal variants of the traits under study.

A genome-wide association analysis for traits under study was conducted using BLINK and FarmCPU models. The Bonferroni correction threshold value of –log10 > 7.0 (p-value) was used as the cutoff, which is highly significant to identify the significant SNPs associated with the studied traits. The markers considered to be significantly associated with tested traits were represented by illustrating the Manhattan plots (Figures 5, 6).

This study identified a total of 632 significantly associated SNPs by both BLINK and FarmCPU models (Supplementary Tables 2–7) and Manhattan plots (Supplementary Figures 1–3). MTA for the trait NON resulted in the identification of 27 SNPs from BLINK and 38 SNPs from FarmCPU models having 15 markers common between these two models (Table 2).

The MTAs Ca7_32613892, Ca2_825900, and Ca2_825902 explained phenotypic variance (PVE) with 48.80%, >25%, and >25%, respectively (Table 3). The MTAs of NFW demonstrated 96 SNPs from the BLINK model and 3 SNPs from the FarmCPU model (Supplementary Table 3), and the results from the BLINK model are shown in Table 4, as well as Manhattan plots along with Q–Q plots (Figures 5, 6).

The MTAs Ca1_28905467 and Ca1_29852105 explained 3.2% and 9.2% of PVE for respective SNPs of NFW. The MTAs of NOP recognized 43 SNPs from BLINK and 60 SNPs from FarmCPU, respectively, having 17 SNPs common with the BLINK model. MTAs of NOS identified 16 and 59 SNPs from the BLINK and FarmCPU models, respectively, with 12 SNPs being common. Furthermore, for the trait 100 SW, we have found 4 and 286 MTAs from the BLINK and farm CPU models, respectively. Thus, the maximum number of SNPs was identified for 100 SW followed by NFW, NON, and NOP.

The SNPs found for more than two locations were considered as stable SNPs in our study (Table 3). The two SNPs Ca1_36724701 and Ca7_18587603 present on chromosomes 1 and 7, respectively, belonging to the trait NON were found to be stable at locations 1, 3, and 4. Seven SNPs were found to be stable at locations 1 and 3 for NON. The stable SNP 2_825902, which is found on chromosome 2, explained 27.33% of the PVE. Furthermore, 7 stable SNPs for NFW, 5 for NOP, 6 for NOS, and 32 for 100 SW were identified. The PVE explained by the stable MTAs ranged from 1.11% to 27.33% for NON, 5.7% to 34.46% for NFW, 1.12% to 8.95% for NOP, 1.23% to 11.71% for NOS, 1.23% to 31.45% for seed yield, and 0.92% to 28.65% for 100 SW. Variation in PVE % was observed across different locations. This may be attributed to the genotype-by-environment (G×E) interactions, which influence the relative performance of the genotypes under different environmental conditions.

The SNPs found common for two or more traits were considered as pleotropic (Table 5). The three SNPs (Ca1_10074058, Ca2_825900, and Ca2_825902) were found to be common for the traits NON and NFW. The SNP Ca7_18587607 located on chromosome number 7 was found to be common for the traits NOP, NOS, and 100 SW. We also found 11 significant and common SNPs for the traits NOP and 100 SW. Significant SNPs contribute substantially to the variability of the trait in the population being studied. Researchers prioritize these SNPs for further investigation and consider them as potential candidates for explaining the genetic basis of the trait.

The significantly associated SNPs with different traits were further used for the identification of putative candidate genes based on the position of SNPs and flanking regions (Kawahara et al., 2013). The 10-kb upstream and 10-kb downstream sequences from the SNP positions were retrieved from the NCBI database and functionally annotated based on CDC Frontier v1 functional annotations (Varshney et al., 2013). The SnpEff-4.3T open-source program was also used for variant annotation and prediction of significant SNP effects, and we found that majority of the SNPs in our study were present as intron variants. The Ca7_32113892 identified for NON was present within the genomic regions of the protein encoding a calumenin-like isoform X2 and a calumenin B-like isoform X1. These have a function in NOD factor export and, in turn, are shown to be involved in nodulation. The SNP located on Ca7_45149648 localized with the protein transcription factor GTE4-like protein that functions as an activator of gene expression upon infection with Pseudomonas syringae and helps in the upregulation of salicylic acid (SA)-mediated immune defense genes was associated with NON. The GTE4-like protein is closely associated with bacterial infection and may be involved in the regulation of the nodulation process. Some of the important proteins along with their molecular functions are listed in Table 6.