A panel of 422 hexaploid wheat accessions was examined for molecular characterization analyses. Seed samples were obtained from the Plant Genetic Resources Institute, National Agricultural Research Center, Islamabad, Pakistan. Based on the given samples information, the 422 wheat accessions (hereafter referred to as the whole population, WP) were classified into three subpopulations, including 145 synthetic-derived wheats (SYN-DER, developed by crossing primary synthetic hexaploid wheats with advanced lines and elite cultivars of Pakistan and CIMMYT), 128 commercially released Pakistan cultivars (PC, i.e., genotypes that are unique and stable, and have been selected for agronomic traits), and 149 advanced lines (AL, i.e., group of lines developed for better agronomic characteristics). Detailed information on the 422 wheat accessions can be found in Supplementary Table S1.

Five viable seeds of each accession tested in this study were planted in 5 cm diameter plots. Genomic DNA was isolated and purified from fresh leaf samples of 2-week-old seedlings using the cetyltrimethyl ammonium bromide (CTAB) method (Doyle and Doyle, 1987). DNA of all the samples was sent to the Cornell University Biotechnology Resource Center for GBS genotyping. The GBS method was performed according to the protocol proposed by Poland et al. (2012) using a two-enzyme (MspI-PstI) approach. A GBS analysis pipeline in TASSEL version 5.0 (Trait Analysis by aSSociation Evolution and Linkage) was used for SNP calling (Bradbury et al., 2007). A total of 133,738,39 GBS-SNPs were retrieved based on the “Chinese Spring” reference genome v.1.0 from International Wheat Genome Sequencing Consortium (IWGSC). Identified GBS-SNPs were named as “chromosome number_physical position, that is, 1A_555961328. More than 50% of the GBS-SNPs were removed from the dataset due to their missing rate being higher than 50%. Unmapped GBS-SNPs were also excluded from the dataset. The remaining SNPs were imputed using Beagle version 5.1 with default parameters (Browning et al., 2018). Then, 18,589 GBS-SNPs with heterozygosity <0.2 and minor allele frequency (MAF) exceeding 0.05 were retained using TASSEL version 5.0 (Bradbury et al., 2007) for the follow-up analysis.

PowerMarker version 3.25 software (Liu and Muse, 2005) was used to calculate population genomics parameters, including MAF, heterozygosity (He), genetic diversity (GD), and polymorphism information content (PIC) for the WP and each of the three predefined subpopulations (i.e., SYN-DER, PC, and AL). To investigate the patterns of nucleotide variations, transition (Ts) to transversion (Tv) mutation statistics, Tajima's D tests, and nucleotide diversity (π) were evaluated using VCFtools version 0.1.15 (Danecek et al., 2011). The population structure of the WP was assessed using NJ-tree and PCA. NJ analysis was conducted using TASSEL version 5.0, while PCA analysis was performed using an R “SNPrelate” package (Zheng et al., 2012). TASSEL version 5.0 (Bradbury et al., 2007) was used to perform linkage disequilibrium (LD) among pairs of SNPs of each subpopulation by estimating squared allele frequency correlation (r²) of alleles. The LD decays within WP and three subpopulations were evaluated, as was the distance among pairs of SNPs with non-linear regression using a custom R script.

The relative kinship analysis implemented in the GAPIT (Genomic Association and Prediction Integrated Tool) R package (Lipka et al., 2012) was performed to reveal the genetic identity (or genetic relationship) between any two given accessions. Negative kinship coefficient values between two accessions, indicating the existence of a weaker genetic relationship than would be expected between two random accessions, were set to zero. Roger's genetic distance was estimated using BIO-R software version 2.0 (Pacheco et al., 2016). Negative genetic distance values were replaced by zero. The analysis of molecular variance (AMOVA) and pairwise F_ST analyses were performed using Arlequin 3.5 software to estimate genetic differences between predefined subpopulations (Excoffier et al., 2005).

EigenGWAS implemented in the GEAR software (freely available from https://github.com/gc5k/GEAR), was used to identify genes/QTL that underlying population genetic differences and to detect candidate regions of the wheat genome under selection in any genetic population (Chen et al., 2016). The EigenGWAS is a single marker regression method based on the PCA. It is similar to a typical GWAS method; however, the phenotype is replaced with an individual-level eigenvector (EV) derived from the genotypic data. Briefly, EigenGWAS involved three steps: first, 18,589 high-quality GBS-SNPs were used to generate the genomic relationship matrix for WP; second, the first 10 eigenvalues and their corresponding eigenvectors (i.e., EV1-EV10) were calculated; and third, marker effects were estimated by regressing each GBS-SNP for a selected eigenvector from the second step. More detailed instructions can be found on the “GEAR” software website (https://github.com/gc5k/GEAR/wiki/EigenGWAS). To exclude the effect of genetic drift (Devlin and Roeder, 1999), the p was adjusted using a genomic control factor (λ_GC), denoted as P_GC, and was used to identify genomic regions under selection. To determine the threshold of significance of genomic regions under directional selection, the first EV was reshuffled 1,000 times to simulate the null distribution. The 95th quantile of the 1,000 most significant P_GC was calculated using 1,000 permutations to determine the cutoff. After log₁₀ (p) transformation, 5.0 was applied for -log₁₀(P_GC) of EigenGWAS analyses in all 10 EVs to declare as regions under selection.

To exclude the strong effect of LD, significant SNP loci within 5 Mb both up- and downstream based on the LD level of the WP were merged as potential selected regions. Functional annotations of the target GBS-SNPs were performed using SnpEff software (Cingolani et al., 2012). The wheat IWGSC RefSeq Annotation v1.0 as a “ggf3” file format was downloaded from the EnsemblPlants database at https://plants.ensembl.org/. The PANTHER Overrepresentation Test (release 20210224) using Gene Ontology (GO) database (release 2021/05/01) using Triticum aestivum database as a reference list. GO analysis included biological process, molecular function, and cellular component. The raw p < 0.01 was set as the threshold to declare the significant differential expression.

A total of 18,589 high-quality GBS-SNPs were well distributed across the genome in the current diversity panel (Table 1). The GBS-SNPs covered a physical distance of 14,053.03 megabase (Mb), with an average density of 1.26 Mb per SNP. The number of GBS-SNPs identified were 7,423 (40%), 9,035 (49%), and 2,131 (11%) in A, B, and D genomes, respectively (Table 1). Among genomes, the highest and lowest number of GBS-SNPs were recorded on chromosomes 2B (1575 SNPs) and 4D (128 SNPs), respectively. The marker density for the D-genome (0.53 Mb per SNPs) was lower than that for the A- (1.51) and B- (1.74) genomes (Table 1). Chromosome-wise marker density varied from 0.25 (4D) to 2.05 (7A and 7B). Chromosome sizes ranged from 473.05 Mb (6D) to 829.74 Mb (3B).

The genetic diversity parameters including MAF, He, GD, and PIC were calculated using 18,859 GBS-SNPs for the WP and each of the three predefined subpopulations per chromosome in each genome of the panel. The frequency distribution of GBS-SNPs for MAF, He, GD, and PIC is presented in Figure 1. The details of GBS-SNPs per chromosome and across genomes are presented for values of MAF, He, GD, and PIC for WP and subpopulations (Supplementary Table S2). For the WP, as expected the MAF value across genomes ranged from 0.05 to 0.5 with an average of 0.21. As expected, the subpopulations still contained GBS-SNPs with MAF ranging from 0 to 0.05 (Figure 1A). It suggests that some of the common alleles in the WP were rare (MAF < 0.05) in the subpopulations. The PC subpopulation had a higher number of rare GBS-SNPs than the other two subpopulations (i.e., SYN-DER and AL). The numbers of GBS-SNPs with MAF ranged from 0 to 0.05 were 2,015, 2,513, and 1,034 in SYN-DER, PC, and AL subpopulations, respectively (Figure 1A). Viewing the WP, rates of GBS-SNP heterozygosity varied from 0 to 1.99, with an average of 0.019 (Figure 1B). The averaged heterozygosity rate for subpopulations was 0.010 (SYN-DER), 0.036 (PC), and 0.012 (AL) (Supplementary Table S2). The GD values in A, B, and D genomes were 0.319, 0.312, and 0.263, respectively (Supplementary Table S2). At the subpopulation level, SYN-DER (0.294) and AL (0.294) had the highest GD while PC (0.281) had the lowest GD (Supplementary Table S2). The average PIC values varied from 0.09 to 0.375, with an average of 0.245 in the WP (Figure 1D and Supplementary Table S2). The numbers of GBS-SNPs with PIC values ranging from 0.2 to 0.4 were 12,301 (66%), 12,056 (65%), and 12,977 (70%) for SYN-DER, PC, and AL, respectively (Figure 1D).

Two types of GBS-SNPs were determined according to nucleotide substitution analysis: (1) transitions (A/G and C/T) and (2) transversions (A/T, A/C, G/T, and C/G). Transition-type SNPs (73.24%) were more frequent than the transversions (26.76%), and transition/transversion (Ts/Tv) ratio was 2.73 (Table 2). The C/T transition type (36.80%) had the highest frequency, while the A/T transversion type (3.67%) had the lowest frequency among all six SNP types of nucleotide substitution (Table 2). The frequencies were almost similar between A/C and G/T. To better understand the patterns of nucleotide variations within a population, we also evaluated the genome-wide nucleotide diversity (also known as π) and Tajima's D statistics based on 18,589 high-quality GBS-SNPs in the WP and among three subpopulations (Figure 2). The mean nucleotide diversity across windows for the WP was estimated at 8.28E-07 (Supplementary Table S3). Nucleotide diversity was high in the telomeric regions than in the pericentromeric regions of all 21 chromosomes (Figure 2A). Furthermore, among the three genomes, D-genome had the lowest nucleotide diversity π and Tajima's D statistics than the A- and B-genomes (Figure 2). The AL (8.32E-07) subpopulation showed the highest nucleotide diversity, followed by SYN-DER (8.12E-07) and PC (7.93E-07) Supplementary Table S3). The average Tajima's D statistics across windows was positive (i.e., 1.42) for the WP (Figure 2B and Supplementary Table S3). On the other hand, the average Tajima's D statistics in SYN-DER, AL, and PC were 1.13, 1.21, and 1.04, respectively. The mean Tajima's D values were positive for the WP and all subpopulations (Supplementary Table S3), reflecting populations may have gone through balancing selection.

To investigate possible population structure and genetic relationships among the 422 wheat accessions, NJ-tree, PCA, and kinship analyses were conducted (Figure 3). Based on the NJ-tree analysis, we found that the three subpopulations (SYN-DER, PC, and AL) were separated with some admixture (Figure 3A). The AL and SYN-DER were more scattered over PC-1, while PC were more separated along the PC-2. We also inferred the genetic structure and relatedness among the WP by PCA analysis (Figure 3B). In PCA, the first and second principal components explained 9.12 and 4.97% of the total variation, respectively. The weak population structure was detected by both NJ-tree and PCA in the diversity panel, as revealed by Figures 3A,B. The SYN-DER and AL subpopulations were relatively more scattered than the PC subpopulation, indicating that there exists broad genetic divergence in the present collection (Figure 3B). A low level of population structure was also supported by the VanRaden kinship analysis (Figure 3C), which was in accordance with NJ-tree and PCA analyses. The kinship coefficient between pairs of 422 accessions ranged from 0.00 to 3.42 (Figure 3C).

To gain further insights into genetic relationships, we calculated the frequency distribution of kinship coefficients and genetic distances for the WP and the three subpopulations (Supplementary Tables S4, S5). Kinship coefficients near zero indicate no relationship, while those near 2.0 indicate a closer relationship (Supplementary Table S4). For the WP, 60% of the kinship coefficients were equal to 0, 39% varied between 0.01 and 0.8, and the remaining 1% fell between 1.2 and 3.4 (Supplementary Figure S1A). The proportion of kinship coefficients equal to 0 in SYN-DER, PC, and AL were 65, 63, and 60%, respectively (Supplementary Figure S1A). The pairwise genetic distances among the 422 accessions varied from 0 to 0.68 (Supplementary Figure S1B and Supplementary Table S5), with an average of 0.53. The genetic distance between pairs of accessions fell in the range of 0.50 to 0.70 were 82.23% of the WP, 90.06% of the SYN-DER, 90.06% of the PC, and 86.57% of the AL (Supplementary Figure S1B). Most accessions had an estimate between 0.50 and 0.60, regardless of the subpopulations.

Genetic differentiation of predefined subpopulations (i.e., SYN-DER, PC and AL) was assessed using AMOVA analysis (Table 3). AMOVA results showed that 3.41% of the total variation was attributable to the differences among subpopulations, whereas 90.74% was within subpopulations (Table 3). Furthermore, pairwise F_ST analysis was computed to investigate subpopulation divergences and presented in Table 4. The F_ST coefficient among subpopulations varied from 0.0492 to 0.075. The F_ST coefficients showed that the divergence between the SYN-DER and AL (0.0492) was lowest, while the divergence between PC and AL was highest (0.075). Results suggest a low level of genetic differences and in accordance with the NJ-tree and PCA analyses.

The summary statistics results for each chromosome in each genome of LD between adjacent GBS-SNPs were computed in the three subpopulations and the WP (Supplementary Table S6). The average r² values ranged from 0.06 (5D) to 0.33 (4B). The average r² for WP was found to be 0.19. The averaged r² reached the lowest in the AL subpopulation (0.06), and the highest in the AL subpopulation (0.35) (Supplementary Table S6). The averaged r² was decreased with an increase in distances of the genome for all the subpopulations, suggesting that the probability of LD was low between widely separated SNP pairs (Figure 4 and Supplementary Figure S2). The LD decays at 8.52, 5.79, 8.34, and 6.25 Mb for SYN-DER, PC, AL, and WP at r² of 0.1, respectively (Figure 4 and Supplementary Figure S2).

To identify genomic regions that have undergone directional selection (the so-called “selection signature”) during wheat improvement, EigenGWAS was conducted based on positive and negative coordinates of the 422 wheat accessions from Pakistan on the corresponding EV, and their selection differentiations were quantified by F_ST (Table 5). The average genetic relatedness among 422 Pakistan wheat accessions was −0.0046, suggesting that the effective sample size of the WP was 218.60. The effective number of genome segments was 14.33. The largest eigenvalue was 74.45, while the 10th eigenvalue was 15.61 (Supplementary Table S7). The genomic inflation factor, namely, the λ_GC computed from EigenGWAS, which is commonly used in adjusting population stratification for GWAS ranged from 29.98 to 7.20 (Supplementary Table S7). To facilitate the comparisons, EigenGWAS and F_ST analyses were performed on the WP, and two results were drawn as the Miami plot for each of the 10 EVs (Figure 5). Generally, the peaks from –log₁₀(P_GC) and F_ST fairly mirrored each other, indicating reasonable grouping as defined by F_ST. Overall, EigenGWAS detected selection signatures on all 21 chromosomes (Figure 5 and Supplementary Table S8), while 83 significant GBS-SNPs were identified on 6 of the 10 EVs. To exclude the effect of LD, significant GBS-SNPs overlapping with each other were merged within the 5 Mb genomic window, and highly significant GBS-SNP [i.e., SNP with largest –log₁₀(P_GC)] within one region was used to declare as representative. In total, therefore, 38 selection regions were identified and are shown in Table 5. The total length of the selection regions was 418.97 Mb (Table 5). The distribution selection regions across different chromosomes varied considerably, except for 1B, 1D, 4A, 4B, and 4D. Chromosome 2B (6) and 7B (6) had the highest number of selection regions while 2A (1), 3A (1), 5A (1), 5B (1), 5D (1), 6A (1), 6D, and 7A (1) had the lowest across EVs. Besides, significant selection regions were only identified under EV1 (5), EV2 (24), EV4 (2), EV8 (1), EV9 (5), and EV10 (1) (Table 5). The number of selection regions identified in the B-genome were 19, followed by the D-genome (11), and only 8 selection regions were located in the A-genome. The largest selection region (i.e., 1A_45252081) was identified on chromosome 1A under EV2 spanning roughly over 30 Mb (Table 5). In contrast, a region (i.e., 2B_769950981) on chromosome 2B spanned approximately 5 Mb was the smallest selection region detected under EV8. To understand the biological background of the identified selection regions, we particularly aligned previously reported genes, marker-trait associations, and biparental QTL described for grain yield and yield-related traits, baking quality, disease resistance, adaptation, and flowering-time–related traits (Table 5). Results revealed that 22 (i.e., 57%) out of 38 selection regions were falling within proximity of known functional genes and/or QTL with meaningful agronomic implications as existing support. Results suggest that the 2D region (11.96–21.96 Mb) could be involved in dwarfism in wheat (Rht8) and flowering time (Ppd-D1a.1) (Table 5). The region 3D (567.09–577.09 Mb) consists of grain color gene Myb10-D1, which controls the red pigment of wheat grain. In this region, two QTL for grain yield and kernel width were also reported. The region 1A (0–10.49 Mb) includes one gene and five QTL, which includes a low–molecular weight glutenin subunit Glu-A3 controlling gluten quality of the wheat, while QTL were associated with phenology-related traits such as grain yield, biological yield, flag leaf length, and kernel width (Table 5). The 5A region spans from 36.22 to 46.22 Mb, which includes photo-period responsive gene (Ppd-A1) controlling flowering time in wheat. The 5A region (558.36–569 Mb) consists of QTL associated with spike-related traits in wheat such as spike number, awn length, and spike length (Table 5). Similarly, two regions 5D (436.03–446.03 Mb) and 6B (0–9.52 Mb) also included QTL associated with awn length. A region 7D (42.66 −52.66 Mb) encompasses the Lr34 that is known to be associated with leaf rust resistance (Table 5). Notably, the association of selection regions with previously known genes/QTL is speculative; however, further pieces of evidence are required to validate the present results. Allele frequencies of selected regions across three subgroups are presented in Supplementary Table S9.

Functional annotation was carried out to evaluate the genome composition (e.g., intergenic, exon, intron, UTRs) using the whole-genome SNPs (i.e., 18,589 GBS-SNPs) and significant SNPs identified by EigenGWAS (i.e., 38 selection regions) (Figure 6). Of the whole-genome GBS-SNPs, over one-third were located in the intergenic region; more than 15% were in the regions of transcript (i.e., 19%), downstream (18%), and upstream (16%), respectively (Figure 6A). A similar proportion of genome composition could be observed from the gene annotation of 38 selected regions (Figure 6B). Functional enrichment analysis based on genes (IDs recognized by Panther Classification System) within the selected region was performed to identify possible biological pathways associated with the differentially expressed genes (DEGs). Of the 7,263 GO terms annotations, 4,010 GO terms were in the biological function, 2,358 GO terms were in the molecular function, and 895 GO terms were in the cellular component. The distribution of most significantly enriched GO terms revealed several important processes as catalytic activity (GO:0003824), Adenyl ribonucleotide binding (GO:0032559), iron ion binding (GO:0005506), molecular function (GO:0003674), cellular respiration (GO:0045333, GO:0009060), phosphotransferase activity, alcohol group as acceptor (GO:0016773), transferase activity (GO:0016740), ATP binding (GO:0005524), carbohydrate derivative binding (GO:0097367), ribonucleotide binding (GO:0032553), nucleotide binding (GO:0032559, GO:0030554, GO:0032553, GO:0017076, GO:0032555, and GO:0000166), nucleoside binding (GO:1901265 and GO:0035639), NADH dehydrogenase (ubiquinone) activity (GO:0008137), oxidoreductase activity, acting on NAD(P)H (GO:0016651), and so on (Supplementary Figure S3). Significant SNPs identified by EigenGWAS were also subjected to GO enrichment analysis and shown in Supplementary Figure S4. Annotation of DEGs revealed that they were involved in chloroplast organization (GO:0009658), plastid organization (GO:0009657), response to far-red and red light (GO:0010218, GO:0071489, GO:0010017, GO:2000030, GO:0010114), regulation of photomorphogenesis (GO:0010099), lipid modification (GO:0030258), cellular response to light stimulus (GO:0071482), cellular response to radiation (GO:0071478), cellular response to abiotic and environment stimulus (GO:0071214 and GO:0104004), response to salicylic acid and gibberellin (GO:0009751 and GO:00009739), and so on (Supplementary Figure S4).

In this study, the 18,589 high-quality GBS-SNPs were identified across the three wheat genomes (A, B, and D) using stringent filtering criteria, and used for downstream analysis (Table 1). In concordance with previous studies based on different types of molecular markers such as GBS-SNPs, 90K SNP array, RFLP, SSR, AFLP, and DArT markers (Liu and Tsunewaki, 1991; Röder et al., 1998; Peng et al., 2000; Chao et al., 2009; Nielsen et al., 2014; Voss-Fels et al., 2015; Alipour et al., 2017; Eltaher et al., 2018), we found that a high level of GBS-SNPs were located in the B-genome, while low levels were located in the D-genome, indicating that D-genome is the least diverse wheat genome. Furthermore, Dubcovsky and Dvorak (2007) concluded that a large proportion of natural gene diversity in hexaploid wheat came from the polyploid nature of its tetraploid ancestor (AABB) than the diversity found in Ae. tauschii (DD) during domestication. This conclusion could be a good explanation of the high levels of GBS-SNPs in the A- and B-genomes observed in this study (Table 1). The filtered markers spanned a physical distance of 14,053.07 Mb, with an average marker density of 1.28 Mb/marker for the WP, which was slightly lower than previous reports on wheat (Liu et al., 2019). The present study observed more transitions-type SNPs than transversion-type SNPs across three wheat genomes (Table 2), which is in agreement with several previous studies on hexaploid wheat (Alipour et al., 2017; Kumar et al., 2020). The abundance of the transition-type SNPs was due to the mutation of methylcytosine to uracil and then into thymine (Alipour et al., 2017). The hexaploid wheat genome is highly methylated because it arose from two polyploidization events, which may explain transition-type SNP abundance in wheat. Moreover, various studies support the fact that transition-type SNPs are preferred over transversion-type SNPs, in addition to InDels or multiple allelic SNPs for SNP array development (Bianco et al., 2016; Clarke et al., 2016). The higher Ts/Tv ratio improves the accuracy of SNP prediction with a greater level of confidence.

The molecular characterization of genetic resources remains the most promising option for efficient conservation and sustainable use of their diversity in crop breeding (Alipour et al., 2017; Liu et al., 2019). Genetic variability in Pakistan wheat panel revealed by GD and PIC reflected genetic diversity at the nucleotide level of a genetic population which is a key to understanding the effect of past selective forces on germplasm resources. The average GD (0.30) and PIC (0.25) were estimated in the WP in this study, which is fairly similar to GD and PIC values in previous investigations on wheat (Eltaher et al., 2018; Mourad et al., 2020). On the other hand, Kumar et al. (2020) reported higher GD and PIC values for a set of 483 spring wheat genotypes from India genotyped with 35K Axiom Wheat Breeder's Array. In the present study, AL and SYN-DER subpopulations had higher GD than the PC subpopulation, possibly reflecting recent breeding progress in the diversification of germplasm resources (Supplementary Table S2). Similarly, the higher PIC value was also noted for the AL subpopulation, followed by SYN-DER and PC subpopulations (Supplementary Table S2). Moreover, considerable variation was also noted within the different subpopulations for diversity among the three wheat genomes. As expected, the D-genome showed the lowest genetic diversity for all three subpopulations (Liu et al., 2019). These observations were further supported by nucleotide diversity π and Tajima's D analysis (Figure 2). The differences in genetic diversity among AL, SYN-DER, and PC subpopulations indicated that AL and SYN-DER subpopulations were relatively more diverse. This might be because AL and SYN-DERs subpopulations were developed by crossing Pakistan and exotic parental genotypes as parents (Supplementary Table S1), and have been selected in the field for agronomic superiority (Afzal et al., 2019). It is also possible that artificial selection has fixed targeted regions and resulted in genetically homogenous individuals. Consequently, the genotypes present in AL and SYN-DER subpopulations can be used to enhance genetic variation for selection and to accelerate wheat improvement.

The main challenges associated with the analysis of any genetic data are (1) to explore whether the studied population is genetically homogeneous or contains distinct subgroups, and (2) to find quantitative evidence that supports the presence of these subgroups (Patterson et al., 2006). In this study, NJ-tree, PCA, and pairwise kinship analyses were used to assess the population structure of 422 wheat accessions (Figure 3). Interestingly, these different analyses could not differentiate accessions from predefined subpopulations, which raises the possibility of exchanging adapted germplasm in crop improvement activities within the country. It has been widely reported that foreign wheat genotypes (e.g., Mexico and United States) have extensively been used as parents in Pakistan crossbreeding programs that lead to new cultivars (Ain et al., 2015; Rasheed et al., 2016; Liu et al., 2019), which was in general agreement with pedigree information (Supplementary Table S1). Furthermore, the targets of Pakistan crossbreeding programs included improvement of yield potential; resistance leaf and stripe rust; and tolerance to drought, salinity, and terminal heat stress (Rasheed et al., 2016), which could be another reason for overlap between accession from predefined subpopulations. The present results are generally consistent with several other studies (Rasheed et al., 2016; Afzal et al., 2019) which reported overlap between bread wheat cultivars (or landraces) and synthetic derivatives based on genotypic data. However, it was also noted that advanced lines derived from synthetic wheat were separated from non-synthetic wheat.

The resolution in terms of similarity, the coefficient of kinship matrix is dependent on the number of genotypes and markers used in a study. The low numbers will restrict the exploitation purposes in deciphering novel alleles for economic traits and will exhaust in the short term. Large numbers of both genotypes and markers will increase the dissimilarity coefficients, and this can give a possible overview of the collection in use. It is recommended for long-term breeding goals to explore genetic relatedness and divergence among genotypes and to subject for high-density genotyping (Kumar et al., 2020). Regarding kinship, 60% of the pairwise kinship estimates were equal to zero, indicating that these accessions were distantly related (Supplementary Figure S1A). The kinship estimates identified in the present study could be useful to avoid inbreeding. The average Rogers' genetic distance was larger for SYN-DER as compared with AL and PC (Supplementary Table S4). Approximately, 82% of pairwise comparisons of genetic distance among 422 accessions were in the range of 0.50–0.70 (Supplementary Figure S1B). Consequently, these results provide evidence of a very low degree of genetic redundancy with this diversity panel and support our conclusion that the AL and SYN-DER subpopulations are genetically diverse (Figure 3). Furthermore, AMOVA suggested a high degree of genetic diversity within subpopulations and a low degree of variation among populations (Table 3). These variations were highly significant according to the partition value (p < 0.001). The selection for agronomic traits in the Pakistan crossbreeding programs was considered the main reason for this high variation within subpopulations. The low degree of diversity among subpopulations could be due to high genetic exchange or gene flow (Eltaher et al., 2018; Kumar et al., 2020). Therefore, breeders can select genotypes as parents in crossbreeding for improving economic traits, from the same subpopulation than selecting from different subpopulations. Pairwise F_ST among subpopulations revealed moderate genetic differentiation (Table 4), which is in agreement with population structure analysis (Figure 3). In the present study, the low level of F_ST was found between AL and SD, indicating a low-level genetic differentiation between these two subpopulations (Table 4). This coincided with the AMOVA result (Table 3), where a large proportion of genetic diversity was accounted for within subpopulations. LD reflects the degree of linkage between loci, referring to the non-random association of two or more loci in the genome, and influences the genetic forces that structure a population (Morrell et al., 2012). LD decay is one of the most important factors in evaluating the marker coverage to determine the resolution of GWAS results. It is well reported that different populations and different genomic regions of chromosomes always show varied LD, in accordance with the results of the present study and with previous reports (Morrell et al., 2012; Afzal et al., 2019; Liu et al., 2019). In the present collection, the overall genome-wide LD decay was shorter than that reported for other investigations on Pakistan wheat germplasms using different classes of molecular markers (Afzal et al., 2019; Liu et al., 2019). Some researchers have reported low or null decay in diversity in different germplasm resources (e.g., landraces and modern cultivars), although they observed the impact of breeding on LD patterns and allele frequency (Taranto et al., 2020). The distances of LD decays in the SYN-DERs and AL were higher than in the PC due to SYN-DERs and AL germplasm under high directional selection pressure than in the PC. In all three subgroups, the mean r² value was higher than for the entire population (Supplementary Table S6), indicating that more alleles are in LD with a weak population genetic structure.

The genetic bottleneck is an important challenge in crop breeding and artificial selection, which also eliminates standing variation of surrounding genomic regions. Identification of genomic regions for artificial selection is a basic step in understanding breeding history (Liu et al., 2019; Li J. et al., 2020). The trait-associated genes/QTL identified within selection regions should have undergone selection during wheat breeding activities and could be valuable for marker-assisted selection of traits useful for agriculture and assist the use of germplasm. In the present study, 38 genomic regions were found in the present diversity panel and were compared with previously known genes and reported QTL in different wheat populations (Table 5). From population structure analysis, there is no clear separation among the three subpopulations (Figure 3), and in agreement with pairwise F_ST calculated among the three subpopulations (Table 4). Therefore, EigenGWAS approach was chosen because it does not require predefining the subpopulations (Chen et al., 2016). An eigenvalue reflects the mean genetic variation captured and was used as the phenotype in EigenGWAS. Whereas λ_GC of eigenvector in EigenGWAS represents median of variation (Supplementary Table S6). Moreover, the difference between eigenvalue and λ_GC is equivalent to the difference between the mean and a median of a population, implicating the existence of strong selection that could be due to natural or artificial selection during domestication or breeding (Li J. et al., 2019; Liu et al., 2019). In recent studies, the EigenGWAS combined with F_ST analysis has been deployed to identify wheat selection regions (Afzal et al., 2019; Liu et al., 2019). Liu et al. (2019) detected genomic regions by wheat 90 K SNP array in 687 accessions, mainly collected from Pakistan and China, and found that most of the selected regions were associated with known phenotypes for disease resistance, vernalization, quality, and adaptability traits. Afzal et al. (2019) investigated 240 Pakistan wheat accessions, including 171 accessions for SYN-DERs and 69 accessions for PC and AL with wheat 90 K SNP array, and found 89 selection regions within the proximity of functional genes associated with phenology-related traits such as Vrn-D3 and TaElf3-D1 for flowering time, TaCwi-A1, TaCKX-D1, TaSus1-7A, and TaGS-D1 for grain size and weight. However, the use of 90K SNP array could lead to ascertainment bias because the representative SNPs were discovered mainly from the wheat cultivars from Australia, United States, and Europe. The majority of selected regions identified in the present diversity panel were located in B genome (50%) as compared to D-genome (29%) and A-genome (21%), a finding consistent with previous reports (Afzal et al., 2019; Liu et al., 2019), reflecting that the B-genome has experienced intense selection pressure than the D-genome. A total of 38 selected regions were identified in the present study. Of which 22 selected regions were overlapped with previously reported functional genes or/and QTL for important agronomic traits including yield-related traits (QGy, QGw, QTkw, QKw, QSn, QKl), plant height (Rht-8, QPh), end-use quality (Glu-A3), flag-leaf-related traits (QFla, QFlw, QFll), biotic resistance (Lr34), and vernalization (Vrn-D3 and Ppd-A1) (trait nomenclature is presented in Table 5 legends; Table 5). Similar investigations reported in other crops such as maize, barley, and soybean, also suggested that most of the selected regions are associated with phenology-related traits (Liu et al., 2017; Li J. et al., 2019, 2020; Li Z. et al., 2020). Our findings suggest that the selected regions observed in the Pakistan diversity panel may be (or have been) under direct selection and are plausible because it reflects wheat breeding targets in Pakistan. These selected regions will be of interest to further understand their contribution to crop improvement and adaptation of Pakistan wheat germplasm resources.

Although many selection regions had been identified in Pakistan wheat germplasm before (Afzal et al., 2019; Liu et al., 2019), the function of many genomic regions remains unclear. Thus, it is crucial to gain more information on the architecture of selected regions. The functional annotation results revealed that most of the loci were mapped to intergenic regions than that to coding regions (Figure 6), in agreement with previous reports (Jordan et al., 2015). Several studies reported that intergenic regions are genetically diverse and are associated with phenotypic variations (Mei et al., 2018). Viewing the whole-genome GBS-SNPs, the GO analysis revealed that inferred genes were mainly associated with molecular functions (e.g., catalytic and enzymatic activity), the biological process of protein phosphorylation, cellular respiration, aerobic respiration, signal transduction, and cellular component (e.g., Photosystem II reaction center) (Supplementary Figure S1). Whereas, annotation results of selected regions showed inferred genes mostly encoding chloroplast and plastid organization, lipid oxidation, cellular response to red or far-red light, cellular response to abiotic, environment stimulus, response to gibberellin, response to salicylic acid, and so on (Supplementary Figure S2). For instance, response to salicylic acid (GO:0009751) was a significant GO term, which controls the growth and stress response (e.g., drought) in wheat (Loutfy et al., 2012). Liu et al. (2019) identified drought tolerance genes (NAM-6A, and 1-FEH-w3) within the selection regions in Pakistan wheat germplasm, which supports our observation. Comprehensive knowledge of genetic diversity, population structures, and the identification of selection regions offer the potential to assist plant breeders in better understanding the implications of the selection regions on targeted crop improvement and facilitate the use of germplasm.