Wheat (Triticum aestivum L.) is amongst the most prominent cereal crops in the world. Wheat is cultivated on more than 224 million hectares (Mha) and provides 20% (Bhatta et al., 2017) of dietary nutrients for 40% world population (Ray et al., 2013; Arzani and Ashraf, 2017; Food and Agriculture Organization of the United Nations, 2019). Moreover, 68% of global wheat is directly utilized for human food while the remaining 32% for livestock consumption and other miscellaneous uses such as biofuel (Mohanty and Swain, 2019). The world population is increasing at the rate of 1.5% per annum and is expected to reach 9.7 billion by 2050 (Lobell et al., 2011). However, a 0.9% annual increase in bread wheat yield currently is not sufficient to ensure global food security (Ray et al., 2013). Pakistan is the fifth densely populated country in the world. Wheat was cultivated on an area of 8,825 thousand ha of land, which produced 24,946 MT of grain in 2020 (Pakistan Economic Survey 2019-20).

Wheat genetic gains can be improved by optimizing crop husbandry (Sener et al., 2009) and harnessing genetic diversity in native germplasm resources (Nielsen et al., 2014; Govindaraj et al., 2015). Dissection of genetic diversity is a prerequisite for plant breeding experiments, such as domestication, inheritance, conservation, and evaluation of wheat germplasm (Peterson et al., 2014). Narrow genetic bases, limited genetic diversity, and continuous reduction in arable farmland, as well as various climate-associated anomalies in the form of biotic and abiotic stresses, pose a continuous threat to world food security in developing countries (Fischer et al., 2014). Moreover, wheat has a long evolutionary history from the days of its early domestication from Einkorn (one of the primitive wheat ancestors) to modern bread wheat and assimilates plenty of genetic variation during this long period of evolution. Farmers’ selection, uniform varietal seed production, continuous selfing, use of modern breeding techniques, domestication, and stringent selection pressure lead to prompt genetic erosion and squat gene pool, which causes continuous losses of favorable alleles in currently used wheat germplasm (Haudry et al., 2007; Sofalian et al., 2008).

Strategic trait-based wheat breeding is the most viable and sustainable solution for crop improvement. The success of a wheat breeding program can be directly associated with the availability of valuable genetic diversity in the program (Rufo et al., 2019; Voss-Fels et al., 2019). Genetic resources such as wild relatives, gene bank accessions, landraces, advanced breeding lines, and induced and natural mutants are considered indispensable genetic resources for maintaining genetic diversity and crop improvement (Ogbonnaya et al., 2013; Das et al., 2016). Natural variation selection, wide crossing, new gene introduction, genetic hybridization, induced mutagenesis, and horizontal and vertical gene transfer can contribute to enriching the genetic diversity of modern wheat varieties to meet the challenges of climate change and global food security (Arya et al., 2013).

DNA markers are considered indispensable tools for the genetic characterization of plants. Various types of DNA-based markers, including randomly amplified polymorphic DNA (RAPD), amplified fragment length polymorphisms (AFLP), restriction fragment length polymorphisms (RFLP), and simple sequence repeats (SSRs) (Röder et al., 1995), were successfully used for genetic characterization in plants. Nowadays, single nucleotide polymorphism (SNP) markers generated through next-generation sequencing (NGS,) such as targeted amplicon sequencing (TAS), Illumina bead chip array, DArT, Genotyping-by-Sequencing (GBS), and kompetitive allele-specific PCR (KASP), are gradually replacing other old marker systems due to high throughput and low cost per data point. GBS is a reduced representation sequencing method that identifies SNP for genotyping and discovers new SNP compared to other array-based genotyping technologies, namely, DArT/Illumina bead chip. An array-based genotyping chip targets and identifies a pre-labeled specific number of SNPs markers. GBS-based sequencing could be effectively used in the breeding of complex genome crops without any prior sequencing information or even in the absence of reference genome in many orphan crops (Getachew et al., 2019). GBS is a technology that can reduce genome complexity by using two restriction enzymes (PstI/MspI) and simultaneously discover and genotype genome-wide variations in complex genome crops, namely, bread wheat (Elshire et al., 2011; Poland et al., 2012a; Poland et al., 2012b). Currently, GBS has been successfully used in many crops to unblock the GD including wheat, barley, rice, maize, cassava, potato, and soybean (Elshire et al., 2011).

SNPs are the most abundant polymorphic markers in both plant and animal genomes (Rimbert et al., 2018). The bi-allelic nature, high level of polymorphisms, ubiquitous presence, uniform distribution across genomes, automated data acquisition, and analysis make SNPs the most suitable marker for genome-wide marker analysis (Verma et al., 2015). Quick advancement in NGS techniques significantly improved sequencing throughput and reduced sequencing cost, making the genome-wide SNP analysis a time- and cost-effective tool for genomic studies (He et al., 2014). Therefore, the NGS-based markers have been widely used for retrieving genetic diversity (GD), harnessing population structure (PS), studying linkage disequilibrium (LD), mapping quantitative trait loci (QTL), conducting genome-wide association studies (GWAS), and genomic selection (GS) in various crops (You et al., 2018).

Although the cost of whole-genome sequencing significantly reduced during the last decade, still it is not feasible to completely sequence all genotypes for routine screening of breeding materials especially for those species with a huge genome such as wheat in most breeding programs. Mourad et al. (2020) genotyped a panel of 103 spring wheat accessions that were collected from five continents and deposited in the USDA gene bank to evaluate the GD, PS, and LD patterns. The panel was genotyped with 36,720 high-quality SNPs, and the whole population was divided into three subpopulations on the basis of analysis of molecular variance (AMOVA), structure analysis, kinship, and principal component analysis. No high LD was observed on a whole bread wheat genome, but at the sub-genomic level, the D genome showed the highest LD decay value compared to the A and B sub-genomes.

Alipour et al. (2017) utilized 16,506 polymorphic GBS-SNPs to dissect the GD of an Iranian wheat diversity panel of 369 wheat genotypes. The B subgenome has the highest number of mapped SNPs compared to the other A and D subgenomes, respectively. The whole population was divided into three subgroups: one for cultivars and two for landraces. Heslot et al. (2013) analyzed 38,412 GBS-SNPs in 365 soft winter wheat to harness GD in advanced breeding lines. Yang et al. (2020) characterized a population of 180 bread wheat genotypes from Asia and Europe to determine the indigenous PS and GD using 24,767 high-quality polymorphic SNPs using the GBS approach and to determine GD of the subjected population. The polymorphic information content (PIC) value of markers ranges from 0.1 to 0.4. Based on cluster and structure analysis, the whole diversity panel was divided into two groups: group 1 comprises European and partial Asian and group 2 consists of the Middle East and partial Asian accessions.

Most Pakistani wheat germplasm has not been characterized using DNA (SNP) markers. Local germplasm is always a key source of resistance against biotic and abiotic anomalies. The primary objective of this study is to investigate the extent and the pattern of GD, PS, LD, PIC, phylogeny, and kinship of the subjected population that will be useful for the selection of breeding parents for various stress breeding strategies through GWAS and association mapping (AM).

Genetic diversity among the wheat germplasm is exceptionally important for the genetic improvement of wheat cultivars. Wheat widely grows in South Asia (Pakistan, India, and Bangladesh), producing more than 15% of the world’s wheat production. Pakistan ranks seventh in world wheat production. However, the genetic diversity of most Pakistani wheat germplasm has not been well characterized using genome-wide SNPs markers. Information on the genetic diversity of the germplasm would be crucial for breeding programs that can play a critical role for wheat breeders to make efficient use of available germplasm resources in crop improvement (Hawkes, 1981). In this study, 184 Pakistani spring wheat accessions were genotyped using genome-wide SNPs generated by GBS technology. Most SNPs were mapped on the B genome (42%), followed by the A (36%) and D (22%) sub-genomes (Figure 1A), with the maximum number of SNPs on chromosome 2B (9,126) and the minimum number (2,660) on chromosome 4D (Figure 1B). These findings are consistent with previous studies (Alipour et al., 2017; Bhatta et al., 2017; Eltaher et al., 2018; Rufo et al., 2019; Kumar et al., 2020; Yang et al., 2020).

In the current study, the number of SNPs from the D sub-genome (22%) was about half of these present on the B (42%) sub-genome, which agrees with some of the previous reports that the B genome harbors twice as many of the mapped markers as these present on the D sub-genome (Wang et al., 2013; Iehisa et al., 2014; Edae et al., 2015). However, a much lower ratio of the polymorphic markers was found on the D sub-genome, as low as five times compared to the A or B sub-genomes (Allen et al., 2013; Cavanagh et al., 2013). Low polymorphism in the D genome in most of the germplasm was observed not only for SNPs markers but also for RFLP, AFLP, SSR, and DArT markers (Liu and Tsunewaki, 1991; Röder et al., 1998; Peng et al., 2000; Chao et al., 2007; Akhunov et al., 2010; Würschum et al., 2011; Poland et al., 2012a; Berkman et al., 2013; Marcussen et al., 2014; Nielsen et al., 2014; Shavrukov et al., 2014; Edae et al., 2015; Alipour et al., 2017; Eltaher et al., 2018; Rufo et al., 2019). In wheat, A and B are predecessor sub-genomes that accumulated more genetic recombination, duplication, mutation, and gene flow events than those in the relatively newer D subgenome (Berkman et al., 2013). These results suggest that Pakistani wheat genotypes, like the Iranian germplasm (Alipour et al., 2017), hold more genetic diversity on the D genome than other sources of germplasm. Therefore it is a valuable source for crop improvement against various climate anomalies (Jia et al., 2013). Pakistani wheat accessions as T. aestivum L. ssp. Sphaerococcum indigenous from northern Pakistan and northwestern India an early flowering, yellow rustresistant, semi-dwarf plant with a semispherical grain shape (Percival, 1922). This notion points out the ancestral relationship between global populations of Puccinia striiformis f. sp. tritici (Pst) and the putative origin of Pst in the Himalayan, as well as their neighboring plains and foothill regions. The presence of putative alleles, high level of diversity, ability to have the sex-related structure in (Pst) races, and clear signature of recombination further support this hypothesis (Ali et al., 2014). Archeological remains and evolutionary studies indicated that hexaploid spring wheat was already grown in the region of India and Pakistan somewhat between 4000 and 2000 BC, which indicates that the region (India and Pakistan) has the longest cultivation and interactions history with Pst. This statement also indicates that Pakistani wheat accessions were evolved with comparatively more recombination and ultimately enriched source of genetic diversity (Tengberg, 1999; Habib et al., 2020).

PIC is another important parameter for the selection of markers for breeding programs. Botstein et al. (1980) coined a scale range from 0 to 1.0 to categorize multilocus markers such as SSR according to their PIC value. In the current study, PIC values ranged from 0.004 to 0.37, which is low to moderately high with the mean of 0.158 for the whole genome that was smaller than 0.25 reported by previous literature (Chao et al., 2009; Novoselović et al., 2016; Alipour et al., 2017; El-Esawi et al., 2018; Eltaher et al., 2018; Alemu et al., 2020; Mourad et al., 2020). Among the three sub-genomes, the mean PIC for the D sub-genome (0.14) was slightly lower than that of the A or B sub-genomes (0.16) (Table 1) due to limited hybridization event and gene flow in Aegilops tauschii during evolution (Chao et al., 2009; Lopes et al., 2015; Liu Y. et al., 2017; Eltaher et al., 2018; Tyrka et al., 2021). However, much higher PIC (0.26–0.33) were reported in other studies (Allen et al., 2013; Baloch et al., 2017; Chen et al., 2018; Kumar et al., 2020; Yang et al., 2020).

Gene diversity (expected heterozygosity, He) and PIC are considered as primary measures for the dissection of genetic diversity and shedding light on the mutation rate, as well as the evolutionary pressure on a specific allele in a population over the period of time (Botstein et al., 1980; Shete et al., 2000). Overall GD of a population is mainly explained by the distribution of informative markers on a genome (Nielsen et al., 2014). Genetic diversity provides gene diversity of haploid bi-allelic SNPs and ranges of genetic distance, as well as average heterozygosity among individuals of a population (Nei, 1978). As expected, in the current study, the overall mean GD value was greater than PIC (Table 1). In the absence of more polymorphic alleles and ubiquitously even distribution of allele frequency of markers, PIC is always lower than its GD. The PIC values of SNP markers affect the classification of informative markers: highly informative, moderately informative, and slightly informative. Highly informative markers could be used in studying population genetics and GD in various plant species (Salem et al., 2021).

Understanding GD and population structure is a prerequisite to comprehending the genetic variability of germplasm before it can be used in a biotic or abiotic stress breeding program. In the current study, genome-wide mean GD was 0.18, ranging from 0.005 to 0.97 in the panel evaluated, which is slightly lower than 0.26 by Yang et al. (2020) and 0.29 by Mourad et al. (2020) but slightly higher than 0.14 observed by Alipour et al. (2017). The genome-wide heterozygosity (HZ) was 0.0736, ranging from 0 to 0.89. Mean HZ varied greatly among sub-genomes, with the lowest value (0.064) for B and the highest (0.085) for the D sub-genome. A similar trend of HZ values was also reported in other studies (Eltaher et al., 2018; Rimbert et al., 2018; Liu et al., 2019; Kumar et al., 2020).

AMOVA data suggested significant GD within sub-groups (80%), which might result from continuous selection for specific traits under certain environments by wheat breeders. A low level of genetic variability between populations (20%) (Table 2) may be due to schematic gene flow in the form of germplasm exchange among different countries and regions (Arora et al., 2014). The high genetic diversity (80%) within sub-groups means that diversity was attributed to variation within population sub-groups. The results of clustering and AMOVA suggested that inter-cluster cultivars crossing may be useful for developing promising agronomic and disease resistance associated traits in local germplasm. Higher diversity within sub-population than among sub-populations was in line with the findings from other studies (Alipour et al., 2017; Bhatta et al., 2018; Eltaher et al., 2018; Kumar et al., 2020). The Nm value is another important parameter to estimate the gene flow among subpopulations. Nm haploid value of 1 or >1 indicates the low rate of gene flow between subpopulations. In our experiment, the Nm haploid value of 0.497 indicates minor gene flow among three sub-populations because the current germplasm collection does not represent a diverse geographical area. The current population was mainly collected from the Punjab province of Pakistan.

LD is a non-random co-segregation of single/multiple loci among the same or different chromosomes of a genome. The magnitude estimation of LD decay helps determine the resolution of association mapping and the number of SNPs needed for effective association studies and marker-assisted selection (Davey et al., 2011). LD extent varies with sub-genomes or even with chromosomes of the same sub-genome and depends on various factors such as recombination rate, allele frequency, genetic drift, population structure, natural selection, and chromosomal rearrangements (Sukumaran et al., 2015; Habib et al., 2020). Higher values of LD among pairs of markers indicate slower decay rates on a genome or a chromosome, suggesting that fewer markers are needed for GWAS. The current study was conducted on a diverse set of 184 Pakistani spring wheat accessions, and the LD decay at 0.34 Mbp for the whole genome, and 0.29 Mbp, 0.2 Mbp, and 0.07 Mbp for all subgenomes A, B, and D, respectively, observed (Figure 5). In the current study, the A sub-genome showed the highest value of LD (0.29) which shows relatively fewer markers needed for genome coverage and comparatively slower LD decay due to its evolutionary history of origination. However, at the chromosomal level, chromosome 3A showed the highest value (0.12) (Table 4) of LD and chromosome 1A showed the minimum value of LD decay (0.09).

In the case of the B sub-genome chromosomes, 3B and 2B showed maximum and minimum values of LD decay of 0.13 and 0.01, respectively. In the case of D, the sub-genome chromosomes 2D and 4D showed LD values of 0.16 and 0.01, respectively (Table 4). Chromosomal level LD decay helps determine the level of decay, LD hotspots, and number of markers needed for genome coverage on various regions of chromosomes. Results of this study were consistent with the findings of Chao et al. (2010) and Ladejobi et al. (2018) reported the LD decay of 6.4, 4.5, 4, and 4.98 Mbp for A, B, D, and whole-genome, respectively. Aleksandrov et al. (2021) performed a study on Bulgarian bread wheat germplasm to dissect the population structure and linkage disequilibrium on two populations, including one population of old germplasm and one set of modern semi-dwarf cultivars. Old germplasm showed an LD decay of 3.6, 3.3, 3.1, and 3.3 Mb in subgenomes A, B, and D and the whole-genome, respectively. These results support the results of the current study. Li et al. (2019) reported the lowest LD decay value at the D genome, such as 0.35, 0.75, and 0.25 Mb at the A, B, and D subgenomes, respectively. However, the highest LD value was observed in the B sub-genome of the population. Zhang et al. (2013) reported that the highest LD decay rate for the A sub-genome ranges from 25 to 30 cM and on the B and D sub-genomes, LD decay ranges from 15 to 20 cM). Sehgal et al. (2017) reported the LD decay of values of 10, 4, 8, and 5 cM for A, B, D, and whole-genome of the wheat population, respectively. These results also confirm our results of the highest LD decay at the A sub-genome.

Usually, the highest LD value was detected in the D genome in some previous reports (Edae et al., 2015; Sukumaran et al., 2015; Liu W. et al., 2017; Ayana et al., 2018; Bhatta et al., 2018; Jamil et al., 2019; Sandeep et al., 2020). The higher LD value for the D sub-genome might be due to limited and shallow infusion of Aegilops tauschii with tetraploid wheat in evolutionary history. The slow LD decay rate in the D sub-genome suggests that more markers may be needed for GWAS, association mapping, and marker-assisted selection (MAS) than those that are needed for the A and B sub-genomes. High LD value in the A sub-genome indicates a few markers needed to cover the genome.

Analysis of population structure, PIC, LD, PCoA, and NJ phylogenetic tree was used to dissect the genetic diversity of 184 Pakistani genotypes. The whole wheat population was divided into three distinct subpopulations, with group 1 and group 3 dominated by the local lines and group 2 somewhat intermediate type mixture of local and CIMMYT lines. Most of the local genotypes were distributed in all three sub-groups with CIMMYT lines because local accessions were developed by crossing one or both CIMMYT lines as breeding parents. This shows that a huge amount of exotic blood mixing among local germplasm leads toward the ample amount of genetic diversity in local germplasm. Pakistan continuously faces problems such as sudden heatwaves, unusual rainfall, water shortage, and unexpected weather changes due to prompt climate change. Most Pakistani spring wheat genotypes are prone to yellow stripe rust because pathogens get resistant to hot, humid conditions due to unexpected rainfall at the vegetative stage. The results of population structure, phylogenetic tree based genetic distribution, moderately informative nature, and the number of SNP markers, PIC value, and LD regions of this population could be used for the detection of genetic diversity in the local germplasm, and this high genetic diversity would be used for association mapping studies and the selection of crossing parents in national wheat breeding programs. The results of the current genetic diversity study provided the information necessary to broaden the genetic bases and conservation studies. This information will urge national wheat breeders to speed up wheat breeding by using the genetic diversity of local and exotic germplasm. The results of the current study will help dissect dissecting genetic diversity in local germplasm, as well as for GWAS-related studies, to improve the local wheat germplasm of the country.