The association panel of C. quinoa accessions was grown at three different locations of University of Agriculture, Faisalabad, Pakistan for 3 consecutive years. As heterogeneity is a strong bottleneck for genetic study, each accession was sown in an isolated plot in the field to ensure the homogeneity of the collected seed germplasm. For the genetic analysis, heterogeneity was scored for each plot within a population of each accession. Following the phenotypic cards developed by Stanschewski et al. (2021), plots were excluded for further analysis where >50% of the plants within a plot were segregating, i.e., when the accession’s main phenotype was not observable. Thirty plants for each accession were sown with 15 cm plant spacing in each plot.

The experiment was designed in a randomized complete block design (RCBD) with split-plot arrangements on a total area of 71 ft². The soil was prepared with two plowings (30 cm depth) followed by planking to conserve moisture content for seed emergence. Ridges of 30 cm height were prepared with 75 cm spacing between each ridge. Seeds were sown in a plot manually by dibbing 2–3 seeds for each hole at a depth of 2–3 cm for each ridge with 15 cm plant spacing. Furthermore, four irrigations totaling 330 mm were applied throughout the crop, including the pre-sowing irrigation suggested as ideal for quinoa growing under Pakistani conditions (Akram et al., 2021).

A selective approach was employed to assess phenotypic diversity and ensure homogeneity, focusing on fewer plants while maximizing the uprooting of individuals. The panicles of these selected plants were bagged within 100 μm mesh pollination bags before the flowering stage [at BBCH 60 as described by Sosa-Zuniga et al. (2017)] to avoid heterogeneous seeds. After the flowering stage [at BBCH 70, as described by Sosa-Zuniga et al. (2017)], bags were removed to allow the panicles to expand and grow. After the removal of the bags, each plant was tagged based on phenotypic diversity to ensure the harvesting of panicles for pure seeds.

Eighteen phenotypic traits, including qualitative and quantitative traits, were recorded to assess the variability among quinoa accessions. Phenotypic descriptors (Stanschewski et al., 2021) were used for data recording of qualitative traits. The phenotypic traits included 8 qualitative traits (panicle color, PC; panicle shape, PS; panicle density, PD; panicle leafiness, Pl; stem color, SC; leaf shape, LS; leaf margins, LM; and leaf granule color, LGC and 10 quantitative traits (number of branches, NOB; number of panicles, NOP; plant height, PH; stem diameter, SD; panicle length, PL; days to maturity, DM; days to flowering, DF; 1,000 seed weight, SW; yield per plant, YPP; and productivity, P). The data for days to flowering were collected when approximately half of the plants were at the flowering stage (Tabatabaei et al., 2022). Moreover, days to maturity were recorded when 90% of the plants in a plot attained maturity (El-Harty et al., 2021).

Phenotypic data of qualitative variables were analyzed using PAST software (Version 3.16) for dendrogram and principal component analysis. For quantitative data, hierarchical clustering heatmaps were generated in TBtool software. Spearman correlation analysis was performed in the R package “corrplot” in the R 4.1.3 program. Principal component analysis, eigenvalues, scree plots, and biplot analysis were also performed in the R 4.1.3 program using “factoextra,” “factorMineR,” and “GGbiplot,” respectively. These PCAs were generated from a correlation matrix created in R 4.1.3.

Furthermore, the “Nbclust” package was employed to calculate the optimum number of clusters for individual and variable biplot analysis. The elbow method from the “Nbclust” package was utilized to identify the optimum number of clusters. The data were analyzed to calculate the genetic distance matrix using the Euclidean distance method for creating a cluster dendrogram using the “ggplot2,” “factoextra,” and “ggsci” packages in R 4.1.3.

Fresh and young quinoa leaves (at the juvenile stage) were used for genomic DNA isolation. Fresh leaves were washed with distilled water and ground using a mortar and pestle. The ground material was subjected to DNA isolation using the GeneJET Genomic DNA Purification Kit following the manufacturer’s protocol. The purified DNA was stored at −20 °C for downstream analysis. The purified DNA was run on a 0.8% gel, and its quality and integrity were assessed through a gel documentation system (GDS). The DNA concentration was measured by UV visible NANODROP (8000 Spectrophotometer, Thermo Scientific), and a working dilution of 50 ng/μl was made.

For genotypic analysis, polymerase chain reaction (PCR) was performed on a 96-well thermal cycler (peqSTAR) using 20 ISSR markers (Supplementary Table 2). The PCR mixture of 30 μl included template DNA (50 ng/μl), dNTPs, primer pair (each of 10 μM), MgCl₂, Taq Polymerase, Taq buffer, and nuclease-free water to make the final volume. The PCR thermal profile consisted of an initial denaturation temperature at 94°C for 5 min, followed by 40 cycles of denaturation at 94°C for 1 min, annealing at 51°C for 1 min, and extension at 72°C for 2 min, and a final extension at 72°C for 15 min. PCR products were resolved on a 2.5% (w/v) high-resolution agarose gel (biotech grade, ACTGene) and visualized on a Gel Documentation system (Bio-Rad, USA).

Unambiguous DNA bands were counted and scored in the form of a binary matrix in Excel (MS toolkit) as “1” (presence) and “0” (absence). The collected data were aligned and analyzed using PAST software (Version 3.16). PowerMarker (Version 3.25) software computed major allele frequency, genetic diversity, heterozygosity, and polymorphic information content (PIC). Principal coordinate analysis (PCoA) was performed using DARwin6 software. The similarity matrix was generated in Pop Gen32 (Version 1.32) based on Nei’s original measure. STRUCTURE software (Version 2.3.4) was used to assess the pattern of genetic structure in a population, and STRUCTURE harvester software was used to assess the exact number of subpopulations (K). The number of subpopulations (K) was set from 2 to 10, with 10 replications for each run. The program was run with parameters of 10,000 burn-in periods followed by 10,000 iterations. Furthermore, STRUCTURE harvester software was used to assess subpopulations’ exact number (K) (Earl and VonHoldt, 2012).

The linkage disequilibrium (LD) between all pairs of ISSR markers was calculated with 10,000 permutations using correlation coefficient (r²) in TASSEL software (Version 5.2.90) (Bradbury et al., 2007). The significance (p-value) level of r² was estimated using Fisher’s exact test. The marker pairs in the same linkage group were considered linked markers; otherwise, they were considered unlinked markers. The LD level for linked and unlinked markers was also calculated using TASSEL software. Allelic pairs with p ≥ 0.0001 were considered to be in significant LD. The LD decay was estimated using the “LOESS” function in the R 4.1.3 program, and a trend line illustrating the LD decay was plotted using the r² values against the genetic distance of loci pairs (bp). The crucial threshold below which the LD might be regarded as being caused by physical linkage was calculated, corresponding to the 95th percentile of the distribution of square root (r²). The LD decay value for locus pairs was calculated where the LD curve and the r² threshold intersected.

The phenotypic and ISSR marker data were combined for association mapping. The association study was conducted using a general linear model (GLM) with a Q-matrix and mixed linear modes (MLM) with a population’s genetic structure (Q) and kinship (K) matrix following the permutation tests of 10,000 (Bradbury et al., 2007). A kinship matrix (K) was generated from 20 ISSR markers and 18 phenotypic variables. The Structure Harvester program calculated the population’s genetic structure (Q) matrix. Subsequently, trait-associated markers were subjected to filtration to evaluate their r² values for GLM with significance levels p ≤ 0.01 and MLM with significance levels p ≤ 0.01 and p ≤ 0.05, eliminating markers with lesser statistical significance.

As heterogeneity is a strong bottleneck for genetic study, each accession was sown in isolated plots in the field to ensure the homogeneity of collected seed samples. Among 43 accessions, populations of 18 accessions (CHEN 425, D-12175, CHEN 391, D-12220, CHEN-297, CHEN-128, CHEN-470, CHEN-33, CHEN-71, Ames-13760, PUC-mix-red, Javi, PI-614927, RU-5, D-12014, PI-614885, PI-614882, and Ames-13731) in 18 isolated plots were grown in the field (Supplementary Table 3). Heterogeneity was scored in each plot within the population of each accession. Plots were excluded from the analysis if >50% of the plants within a plot were segregating, i.e., they were observably different, and a plot was excluded when the main phenotype of the accession was not identifiable within the plot.

Fewer plants were selected to assess phenotypic diversity and ensure homogeneity and the maximum number were uprooted. Among the grown accessions, the populations of CHEN-425, D-12175, CHEN-391, D-12220, CHEN-128, CHEN-33, CHEN-470, CHEN-71, CHEN 297, Ames-13760, PUC-mix-red, Javi, PI-614927, and RU-5 were observed to be homogenous within their respective plots. Based on scored percent homogeneity, plants for each accession were selected, tagged, and bagged before BBCH 60 (Supplementary Table 3); however, the remaining plants were uprooted. Moreover, four plots for accessions D-12014, PI-614885, PI-614882, and Ames-13731 showed >50% heterogeneity and were excluded from the experiment.

For phenotyping, cluster analysis was performed on phenotypic data collected from the field-grown accessions and the literature documented for these accessions. The cluster analysis showed that phenotypic data of the tagged plants from the plots of CHEN-425, D-12175, CHEN-391, D-12220, CHEN-128, CHEN-33, CHEN-470, CHEN 297, and CHEN-71 were congruent with the data documented in the literature. However, the phenotypic data for accessions PI-614927, RU-5, PUC-mix-red, JAVI, and Ames-13760 deciphered incongruence with their documented phenotypic data. In cluster analysis, these accessions showed phenotypic divergence from the plants in their respective plots. It might be due to outcrossing and natural hybridization that led to recombinant line development. Therefore, these plants were labeled distinctly and used individually in the genetic study (Figure 1).

The phenotypic traits displayed significant and wide-ranging variations among the 72 quinoa plants. These variations could be important for developing new cultivars with distinctive morphological and agronomic characteristics. The PC demonstrated high variations among mature quinoa plants. Notably, most plants had green panicles (58%) with a high distribution frequency. However, plants characterized by dark-colored or beige/white panicles had the lowest distribution frequency (1%). For PS, 72% of the plants were intermediate, while only 4% were amarantiform. Pl is categorized as less (75%) with a high distribution frequency; conversely, only 6% of the plants had more Pl with the lowest distribution frequency. For PD, 69% of plants were intermediate, 28% compact, and 4% lax. SC showed green as the most frequent color (67%), while the lowest frequency distribution (14%) was observed for yellow stems. Leaf shape was divided into two major categories: 60% of the plants were rhomboidal, and 40% were triangular, which could be polymorphic for the same plant. The most common leaf edges were dentate (57%), with the highest distribution frequency, and only 11% of the plants had entire leaf edges. For LGC, white was the most frequent color in 81% of the plants, while only 1% had purple leaf granules. The Shannon diversity indices (SDI) for eight qualitative traits varied (H > 0.5) for all study descriptors, ranging from 1.706 for white PC to 3.574 for LS.

Qualitative traits and their respective descriptor states, frequency distribution (%), and Shannon diversity index (SDI) among different quinoa plants are summarized in Table 1.

Principal component analysis (PCA) was performed to comprehend how the qualitative traits contribute to the variation among various genotypes. Principal component analysis explained 25.64% of the total variation by PC1 and 20.61% by PC2, accounting for 46.25% of the variance (Figure 2). These axes, which include PC, PS, Pl, PD, LS, LM, LGC, and SC, display significant differentiating qualitative traits. The qualitative traits that contributed more to the variation in PC1 include PC-G (green), Pl-L (less), SC-G (green), PS-In (intermediate), LS-R (rhomboidal), LM-D (dentate), and PD-C (compact) in PC1, while PC-Pi/Pu/R (pink/purple/red), SC-R (red), LM-S (serrate), LS-T (triangular), and PD-I (intermediate) primarily contributed to the variation in PC2. Furthermore, panicle color was highly correlated with stem color among various plants.

Histograms with continuous distribution for all the quantitative variables were generated, displaying a wide range of phenotypic variation across the quinoa accessions (Figure 3). Descriptive statistics, including minimum and maximum range, mean, mean standard error, standard deviation, variance, coefficient of variation, data distribution pattern (skewness and kurtosis), and diversity indices (SDI and E_H), are given in Table 2. For the coefficient of variation, a high level of phenotypic variability was observed in P (49.74%) and Y (45.71%), followed by moderate levels of variability in NOP (29.07%), NOB (27.77%) and PL (22.23%) and low levels of variability in PH (16.21%), DF (10.39%), SD (7.84%), DM (6.58%), and SW (4.74%). Most of the quantitative variables observed a low level of phenotypic variation. For the distribution pattern, skewness and kurtosis coefficients were formulated to construe the nature of the distribution for quantitative variables. The skewness value ranged from 0.5 to −0.5 for all the quantitative variables except for P (1.20). Likewise, all quantitative traits exhibited kurtosis <0, except for P (2.76). The Shannon diversity indices (SDI) for 10 quantitative traits varied (H > 0.5) for all study descriptors, ranging from 4.16 for production to 4.28 for SW and DM. Moreover, the evenness (E_H) results likewise showed a high variation index value for all quantitative traits and ranged from 0.89 to 1.00. The highest E_H value was recorded at 1.00 for SD, SW, and DM, whereas the lowest E_H was 0.89 for P.

Spearman correlation analysis explains the significant relationship among the quantitative traits. The Spearman correlation coefficient (r) at the significance level (p ≤ 0.001) for quantitative traits is represented in Figure 4, which shows the highest positive significant correlations between NOB and NOP (r = 0.96, p ≤ 0.001), P and YPP (r = 0.95, p ≤ 0.001), followed by DM and DF (r = 0.80, p ≤ 0.001). There were also significant positive correlations between DM and PH (r = 0.46, p ≤ 0.001) and between PH and PL (r = 0.66, p ≤ 0.001). However, a negative significant correlation was observed between PH and YPP (r = −0.50), followed by DM and SW (r = −0.45).

For variable principal component analysis, the results showed that DM, DF, and PH primarily contributed to PC1 (Figure 5A1); however, NOB and NOP showed major contributions to the overall phenotypic diversity in PC2 (Figure 5A2). Furthermore, P and Y were the distinctly contributing variables in PC3 (Figure 5A3). Furthermore, the score and loading plots of the PCA plot based on the quantitative data are presented in Figure 5B. In variable PCA, the first two dimensions (PC1 – 23.5% and PC2 – 22%) accounted for the largest variation in the overall variance. The PCA was represented as a correlation circle where quantitative variables were graphically represented by normalized vectors, depicting the correlation and contribution of variables (Figure 5B). In PC1, DF, DM, PH, PL, NOB, NOP, and SW contributed positively, whereas Y, P, and SD contributed to positive variance in PC2. In the correlation circle, each variable’s significance was strongly correlated with the length and direction of the vectors. Two variables (NOB and NOP) showed higher magnitudes with longer vector lengths (more variance), followed by DF, DM, and PH. Conversely, SW, SD, P, Y, and PL showed lower magnitudes with shorter vector lengths (less variance). Depending on the vector angles, two variables (NOP and NOB) had closer vector angles that showed a strong positive correlation. Furthermore, Y, SD, and SW also showed strong correlations with each other but a negative association with PL.

For individual PCA, S33-01, SD121-08, S71-03, S128-04, S470-05, SD122-05, S71-06, SD121-10, SPrecm-03, S33-04, S470-04, SD122-02, SPrecm-07, S71-01, SPMrecm-03, S470-01, SArecm-01, S291-01, S470-09, S128-01, S425-01, SPrecm-02, S128-08, S128-10, SD122-04, S297-02, S33-09, S470-03, S128-02, SPMrecm-04, and SD121-05 plants primarily contributed to PC1 and PC2 (Figure 6A). The score and loading plots of the PCA plot based on the quantitative data of quinoa plants are presented in Figure 6B. In individual PCA score plot, the first two dimensions (PC1 – 23.5% and PC2 – 22%) accounted for the largest variation in the overall variance. The results showed that all plants were effectively separated in all quadrants based on the first two principal components that showed phenotypic variability among them (Figure 6B).

The comparison of eigenvalues and proportion of variances for 10 principle components among quinoa accessions is given in Supplementary Table 4. A scree plot for quantitative traits was created based on the 10 PCs. The results showed that the first four PCs based on the quantitative data (eigenvalue >1) accounted for 76.1% of the total variance among all the accessions (Figure 7). Individually, PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, PC9, and PC10 represented 23.5%, 22%, 18.2%, 12.4%, 9.9%, 6.4%, 4%, 2%, 1.4%, and 0.3% of the total variance, respectively.

For PCA biplot analysis, the elbow method was used for the cluster optimization to delineate the level at which most phenotypic data is retained with well-explained phenotypic variability among quinoa accessions. The results showed that the optimal number of clusters was k = 3 to perform the cluster analysis, which shows the variations based on quantitative traits within the clusters (Figure 8).

Cluster analysis was performed for quantitative data that grouped all the quinoa plants into three distinct clusters (cluster I, cluster II, and cluster III) (Figure 9). The clustering tree (cluster dendrogram) showed that similar plants tend to group in the same cluster. The highest number of plants was grouped in cluster II, which comprised 30 plants, followed by cluster III with 29 plants, although cluster I showed the least number of plants (Figure 9).

In the PCA biplot, PC1 and PC2 accounted for 23.5% and 22% of the total phenotypic variance, respectively. The results showed that the plants from each accession were disseminated across all quadrants of the PCA ellipse plot and presented diverse clustering in the PCA biplot (Figure 10A). However, plants correlated with a particular varietal cluster primarily employed certain quadrants based on quantitative variables. From 72 quinoa plants, 13 plants (S71-04, SD121-07, SD121-08, SD121-09, SD121-10, SD122-03, S128-04, S128-04, S470-02, S470-08, SJrecm-02, SJrecm-04, SRUrecm-01, and SPMrecm-04) were closely grouped in cluster I and predominantly correlated with Y, SW, P, and SD, which showed that these plants might be crucial for generating high-yielding quinoa varieties (Figure 10B). However, most of the plants (S71-02, S71-04, S71-07, S128-03, S128-05, S128-08, S128-09, S33-03, S33-05, S33-07, S33-10, S33-07, SD121-06, SRUrecm-01, SJrecm-01, SJrecm-03, SJrecm-05, SJrecm-06, SPrecm-03, SPMrecm-02, SArecm-01, S470-01, and S470-05) in cluster II were correlated with DF, DM, PL, and PH, demonstrating that these plants were either early maturing or late maturing along with longer or shorter plant height and panicle length. In cluster III, 17 plants (S297-02, S71-06, SPrecm-01, SPrecm-05, SPrecm-07, SMPrecm-01, S33-02, S33-06, S33-08, S470-09, S470-10, S128-02, S128-06, S128-07, S128-08, SD121-04, and S425-01) were highly correlated with NOP and NOB than Y, SW, P, and SD.

Heatmap cluster analysis was performed for quinoa plants and their quantitative traits, including NOB and NOP, DM, and YPP. The heatmap along the dendrogram showed the grouping of quinoa plants based on particular quantitative traits (Figures 11A–C). The clustering heatmaps (cluster dendrogram) categorized quinoa plants into two groups based on NOP and NOB (Figure 11A). Group 1 is the largest subdivided into two subgroups; subgroup 1 included 10 plants that varied in NOP and NOB from 2 to 6, while subgroup 2 included 36 plants that ranged from 6 to 10. Group 2 had 26 plants with the highest number of branches and panicles, ranging from 12 to 16. In group 2, two plants (S71-06 and SPrecm-07) had maximum number of branches and panicles (15). For DM, the clustering heatmap showed two major groups (Figure 11B). Group 1 was subdivided into two subgroups; subgroup 1 included 17 plants that matured early, ranging from 145 to 160 days, while, subgroup 2 included 23 plants that required 160–170 days to attain maturity. Group 2 was further divided into two subgroups. Subgroup 2 included 21 plants that reached the maturity stage in 170–185 days. Subgroup 1 included 11 late-maturing plants, requiring 185–190 days for maturity. Heatmap showed the most early maturing plants were SPrecm-04 and SD122-02 that attain maturity in only 148 days followed by S33-09 (149 days). A clustering heatmap for Y showed two major groups (Figure 11C). Group 1 was subdivided into two subgroups; subgroup 1 comprised 26 plants with the lowest Y value ranging from 5 to 15, and subgroup 2 included 40 plants ranging from 20 to 30. Group 2 included only six plants (SD121-07, S128-04, S33-09, S470-02, SPMrecm-04, and SPrecm-05) with the highest value for Y, ranging from average value reside within 30 to 45 g/1,000 seeds.

The phenotypic data and analysis showed strong correlation of panicle shape, color, density and leafiness (panicle architecture) with yield and productivity in quinoa. The accession SD-121-07 exhibited highest yield per plant (42.72 g; YPP) possessing green, glomerulate, compact panicle with less leaves followed by accessions S128-04 (36.88 g), S33-09 (36.24 g), S470-02 (34.44) with non-significant difference in yield having green, intermediate, intermediate panicles with less leaves. It showed that change in panicle shape and density from glomerulate shape, compact density panicle to intermediate shape, intermediate density panicles resulting in yield reduction. The accessions SJrecm-03 and SJrecm-05 produced less yield, 5.21 and 5.49 g, respectively. Both exhibited pink, intermediate shape, intermediate density panicles with less leaves. In addition to this, the results also inferred that the accessions with panicle shape amarantiform exhibited less yield (Supplementary Table 5).

Linkage disequilibrium (LD) analysis was performed for 20 ISSR markers using pairwise squared-allele frequency correlations (r²) to assess the LD level among different marker pairs (Figure 15). Based on significant LD (p ≤ 0.05, 0.0001), the marker combinations (%) were calculated from a total number of possible marker combinations. As a result, a total of 7,140 pairwise combinations were identified. Among these, 1,654 (23.16%) locus pairs showed significant LD at p ≤ 0.05, while 745 (10.43%) locus pairs showed strong LD at p ≤ 0.0001. Based on r² at a significance level of p ≤ 0.0001, 4.1% and 2.7% of locus pairs showed significant LD with r² < 0.2 and 0.3 > r² > 0.2, respectively. Furthermore, 1.5% and 0.3% of locus pairs were in significant LD, with 0.5 > r² > 0.3 and r² > 0.5 (Supplementary Table 9). LD decay was assessed by plotting the r² value against the genetic distance (bp) of locus pairs, and subsequently, a trend line was created to show the pattern of LD decay (Figure 16). LD relies on r² values, and its critical threshold value (r² = 0.12) was measured corresponding to the 95th percentile of the coefficient square (indicated by the red line in Figure 16), which showed that LD could be due to physical linkage. The LD decay value at the intersection of the LD curve and r² threshold was attained at 5 bp for the whole genome.