Wheat (Triticum aestivum L.) is a staple food crop that feeds about 30% of the world population and provides over 20% of the calories consumed by humans (FAO, 2015). Due to a rapidly growing world population and climate changes, breeders and farmers are facing the challenge of increasing wheat production up to 70% by 2050 to meet future demands (FAO, 2009; Ray et al., 2013; Marcussen et al., 2014), which needs a 2.4% of yield increase yearly. However, the current global average rate of crop yield increase is only 0.9% per year, which is far slower than the desired rate (Ray et al., 2013).

Wheat grain yield can be increased by improvement of both crop management practices and genetic improvement of cultivars for high yield potential (Sener et al., 2009). Genetic diversity is the foundation for such genetic improvement (Nielsen et al., 2014; Govindaraj et al., 2015). Iran is part of the wheat center of origin, and Iranian wheat landraces are important genetic resources for new alleles or genes to be used in breeding for new cultivars (Ciaffi et al., 1992). Wheat landraces are distinct and locally adapted accessions collected and grown by farmers may have a high level of tolerance to biotic and abiotic stresses hence, they are able to provide higher sustainable yields under low input agricultural conditions (Zeven, 1998; Skovmand et al., 2002). For example, Iranian landraces PI 1377397 (Toit, 1989) and PI 626580 (Valdez et al., 2012) have been used as the sources of resistance to Russian wheat aphid [Diuraphis noxia (Kurdjumov)]. ‘Turkey Red,’ a hard red winter wheat from Turkey, has been the foundation for hard winter wheat cultivars in the United States due to its cold tolerance (Olmstead and Rhode, 2002).

After the green revolution in the mid-20th century (Borlaug, 1968), wheat landraces have been widely replaced with modern semi-dwarf cultivars, which significantly narrowed genetic diversity. Although artificial selection performed during domestication and breeding processes has increased the frequency of favorable alleles controlling yield and input responses, some other desirable alleles such as biotic and abiotic stress resistance have been removed from breeding populations and cultivars. Therefore, genetic diversity of bread wheat is diminishing in breeding programs (Sofalian et al., 2008). Conservation of landraces and their wild relatives becomes a critical measure to avoid genetic erosion and meet future need of wheat yield increase (Chen et al., 1994).

Molecular markers have been widely used to study the population structure and genetic diversity of germplasm collections (Huang et al., 2002; Dreisigacker et al., 2005; Hao et al., 2006, 2008; Cavanagh et al., 2013). Single nucleotide polymorphisms (SNPs) are the most abundant type of sequence variations in plant genomes (Batley and Edwards, 2007). They are suitable for analysis of genetic variation, population structure, marker-trait association, genomic selection, QTL mapping, map-based cloning, and other plant breeding applications that need large number of markers to cover entire genomes (Kumar et al., 2012). High-throughput SNP arrays are available, but the high cost per sample limits their application in breeding. Recently, next-generation sequencing (NGS) provides a high-throughput and cost-effective molecular tool for breeding and has been widely used to speed up breeding processes (Poland and Rife, 2012; Edae et al., 2014). Rapid advances in NGS technology have driven the costs of DNA sequencing down to the point that genotyping-by-sequencing (GBS) can now be used for routine breeding screening in any crops (Elshire et al., 2011). GBS can be used for marker discovery and genotyping simultaneously and many samples can be multiplexed to reduce cost per sample (He et al., 2014). GBS uses restriction enzyme digestion to reduce the complexity of genomes, which makes it possible to analyze plant species with large and complex genomes such as wheat, and the low cost per sample makes it feasible for breeders to use it as a routine tool to aid breeding selection. An improved GBS protocol (Poland et al., 2012b) that uses two-enzyme (PstI/MspI) digestion provides a greater degree of complexity reduction and more uniform libraries for sequencing than the original single enzyme protocol (Elshire et al., 2011). The improved GBS protocol has been successfully used in cereal crops such as barley, wheat and oat (Poland et al., 2012a,b; Poland J.A. et al., 2012; Huang et al., 2014).

Iranian wheat landraces provide a rich source of genetic diversity and carry resistance genes to many different biotic stresses such as bunt diseases (Bonman et al., 2015), Russian wheat aphid (Ehdaie and Baker, 1999; Valdez et al., 2012; Bonman et al., 2015), leaf and stripe rusts (Kertho et al., 2015), stem rust (Rouse et al., 2011; Newcomb et al., 2013), and abiotic stresses such as salinity (Jafari-Shabestari et al., 1995), drought and heat (Ehdaie et al., 1988). To date, most of the Iranian germplasm lines have not been characterized and used in modern plant breeding (Hoisington et al., 1999). These germplasm lines not only provide new sources of resistance to biotic and abiotic stresses, but also can enhance the biodiversity of the breeding materials (Huang et al., 2010). Therefore, assessing genetic variation and differentiation of Iranian wheat landraces and cultivars will facilitate the effective use of these valuable genetic resources in future breeding to broaden the genetic diversity of Iranian breeding materials and identify novel alleles that could be used by geneticists and breeders in Iran and other countries. To the best of our knowledge, the present study is the first to directly compare the population diversity of Iranian wheat landraces to a representative pool of Iranian wheat cultivars using SNP markers.

After eight Ion Proton runs of 369 samples, a total of 566,439,207 reads were identified with 458,363,607 (about 81%) unique reads. A total of 133,039 GBS-SNPs were called after filtering out duplicated reads. Among them, 16,506, 38,824, and 56,560 GBS-SNPs have <20, <50, and <80% missing data, respectively. Only the SNPs with <20% missing data were used to evaluate the genetic diversity of the diversity panel. In the BLASTn analysis, a total of 11,758 (∼71.24%) and 14,697 (∼89.04%) SNPs were aligned to CSSS and W7984 reference genomes, respectively. The highest numbers of SNPs were mapped on the B-genome and the lowest on the D-genome in both CSSS and W7984 references. The highest values of transition-type SNPs were identified on the B genome with 3,465 and 4,107 SNPs while the lowest were on the D genome with 1,397 and 1,721 SNPs on the CSSS and W7984 references, respectively (Table 1). A similar chromosome distribution pattern was observed for transversion-type SNPs. More transition-type SNPs (63.63%) were observed than transversion-type SNPs (36.37%) with a transition/transversion (Ts/Tv) SNP ratio of 1.75 (10,503/6,003) over all three genomes. However, the ratio was significantly higher in the A genome than those in the B or D genomes. As expected, more A/G and C/T transitions ware observed than G/A and T/C transitions. On the other hand, more C/G, A/C, G/T, and A/T transversions were observed than G/C, C/A, T/G, and T/A transversions. The SNP numbers in the most of the chromosomes of W7984 assembly were higher than those of CSSS, with only a few exceptions (Figure 1A). Similarly, the average marker density was lower in D-genome than in other two genomes of both references (Figure 1B). The maximum marker density per chromosome of W7984 assembly (≈10.4 SNP/cM) was much higher than that from CSSS assembly (5.5 SNP/cM). The largest marker intervals were observed for chromosome 3D (∼1.77 cM) of CSSS and 4D (∼0.62 cM) of W7984 assembly. Highly significant correlation was observed between the number of SNPs per chromosome and their physical size for both CSSS (Figure 2A) and W7984 (Figure 2B), though this correlation was much higher for SNPs called using W7984 reference than CSS reference. The negative correlation between average marker intervals and chromosome sizes was significant for SNPs called using W7984 reference (Figure 3B), but not significant for SNPs called using the CSSS reference (Figure 3A).

Using the POPSEQ approach, we ordered 11,758 (∼71.24%) and 14,697 (∼89.04%) SNPs to the W7984 and CSSS assemblies, respectively, which were much higher than these reported (22.5–46.3%) in previous studies (Würschum et al., 2013; Shavrukov et al., 2014; Wang et al., 2014). Edae et al. (2015) obtained 33,664 SNPs with up to 80% missing data from W7984 × Opata M85 RIL population out of which, 16,591 (49.3%) and 9709 (28.8%) SNPs were mapped to the W7984 and CSSS reference assemblies, respectively. Distribution of mapped SNPs among the A, B, and D genomes in this study was similar to these in previous reports (Akhunov et al., 2010; Poland et al., 2012a; Würschum et al., 2013; Marcussen et al., 2014; Shavrukov et al., 2014; Edae et al., 2015) with most SNPs mapped to B genome followed by A genome and D genome (Chao et al., 2009; Berkman et al., 2013; Lai, 2015). D genome is the youngest one among the three genomes in wheat evolutionary history. It is likely that older genomes underwent gene duplication and accumulated more mutations that led to sequence polymorphism. Substantial early gene flow could have occurred between T. aestivum and its tetraploid progenitor T. turgidum (AABB) but not between the hexaploid and Aegilops tauschii (DD). This could have resulted in greater sequence diversity in the A and B genomes than in D genome (Talbert et al., 1998; Caldwell et al., 2004; Dvorak et al., 2006; Berkman et al., 2013). In this study, the number of SNPs that mapped to A or B genome were twice as many as those that mapped to the D genome. These results agree with several other studies (Akhunov et al., 2003; Chao et al., 2009, 2010; Wang et al., 2013; Edae et al., 2014; Iehisa et al., 2014) but contradict with some previous reports which indicated that five times more SNPs mapped to A or B genome than to the D genome (Allen et al., 2011, 2013; Cavanagh et al., 2013). This result suggests that Iranian wheat landraces may have relatively higher SNP variation in the D genome than other sources. Higher diversity in the D genome may provide new elite and desirable alleles controlling agriculturally important traits to deal with global climate and environment changes (Trethowan and Mujeeb-Kazi, 2008; Jia et al., 2013).

We observed an average transition/transversion (Ts/Tv) ratio of 1.75 for all three genomes based on both CSSS and W7984 references (1.89 for A genome, 1.76 for B genome and 1.59 on D genome), which reflects the high frequencies of A to G and C to T mutations following methylation. D genome has the smallest Ts/Tv ratio among the three genomes. Higher frequency of transition mutations has been observed in several species including hexaploid wheat (Lorenc et al., 2012; Winfield et al., 2012; Manickavelu et al., 2014) and barley (Turuspekov et al., 2016) with Ts/Tv SNP ratios ranging from 1.59 to 2.12, which agrees with the current study. Transition abundance in many species may result from the mutation of methyl cytosine to uracil and then to thymine (Coulondre et al., 1978). Bread wheat genome is highly methylated due to the two rounds of polyploidy. Therefore, highly repeated sequences and abundant transitions can be considered an ‘evolutionary footprint’ of methylation (Buckler and Holtsford, 1996; Feldman and Levy, 2012; Kalinka et al., 2017). It has been demonstrated that loss of genes occurred more frequently in the A and B genomes than in the D genome (Berkman et al., 2013; Pont et al., 2013). The observed Ts/Tv bias in wheat provides a high level of confidence in SNP prediction accuracy because such a bias was unlikely caused by erroneously called SNPs due to errors in sequencing or mapping (Lorenc et al., 2012).

Chromosome 4D had the least SNPs while chromosome 3B had the highest using both CSSS and W7984 assemblies. On the other hand, a relatively high positive correlation was observed between the number of SNPs mapped to a chromosome and size of the chromosome, which agrees with previous studies.(Poland et al., 2012a; Saintenac et al., 2013; Edae et al., 2015). In addition, other factors such as the time of evolution of a genome also affect the number of SNPs mapped on chromosomes.

Structure analysis placed 369 Iranian wheat accessions into three groups with only 7% of accessions misplaced into an opposite group. The most misplaced accessions were the cultivars that were placed into one of the landrace groups because most of these cultivars were selected from Iranian landraces. Therefore they should belong to landraces. Cluster analysis generated a similar grouping pattern. When pedigrees, geographical regions of cultivation, years of accession released, growth habits and origins of cultivars were analyzed in cluster analysis, we found that accession pedigree is the main factor for separation of the Iranian cultivars. Most cultivars that originated from Iran or with one parent of Iranian wheat were clearly separated from those which originated from CIMMYT, suggesting that Iranian wheat may have a different genetic makeup from the CIMMYT wheat. That may explain why crossing made between Iranian and CIMMYT wheat genotypes produced high yielding cultivars. For instance, Pishgam and Parsi derived from a cross between Iranian and CIMMYT lines are currently the most widely planted cultivars in Iran (Mahfoozi et al., 2009; Amin et al., 2010). For landraces, a high level of genetic diversity was observed among those collected in different years, or from different geographical regions and climate zones. Genetic variation among the landraces collected from north of Iran was higher than the landraces from the south. For example, most of the landraces collected from Gilan and Mazandaran were grouped closely to each other. Landrace collection years can also separate the landraces into two groups based on whether they were collected before and after 1960. This may be attributed to the breeders at University of Tehran who initiated the purification of landraces after 1960. Therefore, genetic variability among the landraces can be affected by the years of collection, anthropogenic impact through dynamic storage practices, diverse needs, end-uses, seed exchange (gene flow) between farmers, and geographic and environmental conditions where wheat grows. It has been documented that genetic diversity of landraces could be related to their adaptability (Cooper et al., 2001) and farmers have played a pivot role in maintaining the genetic diversity (Zeven, 2002).

Genetic diversity for each of the predefined subpopulations was measured using Nei’s and Shannon’s genetic diversity indices, private alleles, the percentage of polymorphic loci and PIC in this study. The Landrace Group I showed higher values of Nei’s and Shannon’s gene diversity indices and the Cultivar Group demonstrated a lower average of private alleles and percentage of polymorphic loci than the two landrace subgroups. This result highlights the strong genetic separation among the groups and specific adaptation of the Iranian cultivars captured in the founder lines of the landraces. But whether the proportion of the genetic differentiation in the subgroups of Iranian landraces and cultivars are due to differentiation during different domestication events or introduction from CIMMYT remains unknown. It is possible that a single domestication event occurred with various hybridizations to endemic genotypes in each region, as has occurred in the domestication of rice (Londo et al., 2006; Huang et al., 2012). Further, the significant genetic differentiation (p < 0.001) among the three groups as illustrated by pairwise F_ST and AMOVA analyses verified the differentiation of three groups in this panel. However, the Landrace Group I demonstrated a higher gene flow and lower F_ST value with the Cultivar Group, suggesting the Landrace Group I might play a bridge role between the Landrace Group II and the Cultivar Group. Most accessions in Landrace Group I originated from cold or moderate and rainy climates, and had been collected before 1960 while those in Landrace Group II were mainly collected from semiarid or hot and dry climates after 1960. Since the landraces were collected from diverse climate regions and altitudes, many of these germplasm lines should be useful sources of genes to be used in breeding to address the challenge of climate change.

Our studies demonstrated that GBS is a powerful tool for investigating population structure and genetic diversity of wheat landraces and cultivars. In this study, Iranian hexaploid wheat landraces and cultivars collected and released from different years, agro-climatic zones with different growth habits were grouped into three distinct groups. Iranian cultivars can be separated into four major groups with one group was mostly originated from Iranian ancestor(s). Large genetic variation was observed among landraces based on the years of collection, geographical and climate zones. We hope that this genetic diversity will help wheat breeders to select parents for crossing to improve wheat under different climate conditions in Iran and other countries.

HA analyzed the data and wrote the paper. GB contributed reagents, Ion Proton sequencing and analysis tools. MB, VM, and SP provided seeds of the genotypes. HA, MB, VM, SP, GB, and GZ contributed to conceive the project and revise the paper. All authors read and approved the final manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.