The plant material comprised 97 accessions including 48 Iranian bread wheat (T. aestivum L.) cultivars and 49 breeding lines provided by the Seed and Plant Improvement Institute, Karaj, Iran (Supplementary Table 1). All accessions were planted at experimental farm of Faculty of Agriculture, Urmia University, Urmia (Iran), during the 2012–2013 growing season. The experiment was conducted on the basis of randomized complete block design (RCBD) with three replications. The experimental plot consisted of a 1-m-long single row with a spacing of 30 cm. During the growing season, the data were recorded on the basis of five individual plants for each plot. Eighteen agronomical traits including plant height (cm), peduncle length (cm), flag leaf length (cm), flag leaf width (cm), flag leaf area (cm²), flag leaf weight (g), stem diameter (mm), number of node, days to flowering, spike length (cm), awn length (cm), spike density, biological yield in row (g), number of spikelets spike⁻¹, number of grains spike⁻¹, grain yield in row (g), thousand grain weight (g) and harvest index (%) were considered in this study based on IBPGR descriptors [16]. Since we were not interested in the variation within accessions due to the highly selfing nature of wheat, the data for each trait were averaged prior to analysis and expressed as plot value per accession. Phenotypic means calculated for each accession were used for subsequent association analysis. Descriptive statistics such as minimum, maximum, mean, standard deviation, environmental coefficient of variation and broad sense heritability were determined for each phenotypic trait.

Genomic DNA was extracted from young leaves of 20-day-old seedlings of five plants from each genotype using the method described by Ausubel et al. [17], with minor modifications.

The concentration and quality of genomic DNA was determined by a spectrophotometer and by 0.8% agarose gel electrophoresis. Primer combinations, PCR conditions, electrophoresis and visualization of PCR products for ISSR and retrotransposon-based markers (IRAP and REMAP) have been described by Dashchi et al. [14] and Nasri et al. [15], respectively. In total, 23 ISSR, 9 IRAP, and 20 REMAP primers were chosen to test for associations between markers and 18 quantitative traits in the current study.

We used the program STRUCTURE version 2.3.4 [12] to investigate the genetic structure (Q) among wheat accessions on the basis of ISSR, IRAP and REMAP markers. The number of groups (K) was set from 1 to 20, and 10 independent runs were performed for each K. A length of 100,000 burn-in and 100,000 MCMC (Markov Chain Monte Carlo) replications were used for the analysis of the genetic structure. The admixture model with correlated allele frequencies was considered. To determine the optimal number of groups, as described by Evanno et al. [18], an ad hoc statistic ΔK was applied. The neighbour-joining clustering algorithm among accessions was run on the basis of the whole number of markers using the TASSEL 3 software [13].

Prior to the analysis of marker-trait associations, all bands with minor allele frequencies less than 0.05 or major frequencies greater than 0.95 were removed from the analysis. Two statistical models were used and compared:

• • the Q model incorporating the Q-matrix as covariate and performed using a general linear model (GLM);
• • the Q + K model incorporating both the Q- and K-matrices using the mixed linear model (MLM) implemented in TASSEL 3.

The Q-matrix was extracted from the STRUCTURE output from the identified optimal group, i.e. K = 4 for the current studied population. The K-matrix is a matrix of pairwise relatedness between accessions calculated from marker data using the TASSEL 3 software. For all traits, the Q + K model performed significantly better than the Q model. Thus, we concluded that MLM approach can significantly reduce the spurious marker–trait associations (type-I error showing false positives). Hence, only the results obtained from MLM were shown and discussed for this study. Significant associations between a marker locus and a trait (i.e. whether a QTL was associated with a marker or not) was indicated by the P-value, and the magnitude of QTL effects was examined by R²-marker. A false discovery rate (FDR) of 0.01 was used as a threshold for significant association using Bonferroni's correction [19].

Basic descriptive statistics for quantitative traits are presented in Table 1. The data showed a wide range of variation for all studied traits. The maximum value of environmental coefficient of variation (CV) observed for awn length (20.72%) and grain yield in row (20.41%). The other remaining traits recorded moderate to low CV values, ranging from 19.82% for thousand-grain weight to 1.97% for days to flowering. Broad sense heritability (h²) is a suitable indicator of the trait repeatability and the proportion of trait variation accounted for by genetic factors. The minimum amount of heritability obtained for number of node (17.70%). Twelve out of the 18 studied traits showed a value of broad-sense heritability of more than 65% (Table 1).

Genome of studied individuals was surveyed using 9 IRAP, 20 REMAP and 23 ISSR primers. A total of 658 clear and reproducible bands were scored, of which 86% (572) were polymorphic among genotypes. No less common loci with frequency lower than 50% were amplified. The number of total loci per primer ranged from 7 (UBC844) to 18 (wilma-bagy2), with an average of 12.65. The number of polymorphic loci per primer ranged from 5 to 18. The size of the amplified fragments ranged from 75 to 2000 bp.

The results of model-based clustering analyses showed the highest peak at K = 4 (Fig. 1), suggesting that the analyzed wheat germplasm can be divided into four genetically distinct groups or populations (Fig. 2). In general, most of the lines located in population A and B, while most of the cultivars were in populations C and D. Individuals with a cluster membership coefficient higher than 0.8 in a population were assigned to that population, while those with a cluster membership coefficient less than 0.8 were considered as mixed individuals. Based on this criterion, three mixed genotypes were identified in each population of A and B. In populations C and D, 8 and 1 mixed samples were detected, respectively. Overall, 15 mixed individuals (15.5%) were found at K = 4. Cluster analysis using the neighbor-joining method also placed the studied germplasm into four main groups, which was mostly in concordance with the results of the model-based clustering method (Fig. 3).

The MLM method was used to identify associations between markers and studied traits. Significant associations between markers and morphological traits along with their corresponding R² and P-value are presented in Table 2. All markers reported in this table were significant at a comparison-wise error rate of 0.01. Ninety-four loci were identified to be associated with studied traits, consist of 41 ISSRs, 10 IRAPs and 43 REMAPs. The highest number of associated markers (11) was found for the number of spikelet spike⁻¹, while the lowest was for the grain yield in row (2). Locus ISSR16-9 (A12) was associated with three traits (peduncle length, flag leaf area, and flag leaf width). Out of the associated IRAP loci, IRAP2-7 (flag leaf area and flag leaf width) and IRAP5-8 (number of spikelets spike⁻¹ and awn length) was linked each with two traits. REMAP13-10 was linked with three traits (spike density, days to flowering, and spike length). Each associated marker explained a maximum regression from 2% (number of grains spike⁻¹) to 17% (spike density) of the total available variation for the individual associated traits.