The spring wheat panel consisting of 111 accessions from the Genebank of IPK in Gatersleben (http://www.ipk-gatersleben.de/en/genebank/) was described in Muqaddasi et al. (2016). The winter wheat panel of 96 accessions was previously described in Neumann et al. (2011). The whole set was field-grown in Gatersleben, Germany (51°49′N 11°16′E, 112 m asl) during 3 years (SP; 2013, 2014, 2015:WP; 2013, 2015, 2016), arranging the accessions as a randomized complete block design (four replicates in 2013 and 2014, and one in 2015 and 2016). The individual plot size was 1 × 1.5 m, and each plot was split into four rows spaced 0.20 m apart. Standard agronomic wheat management practices were applied.

Anther retention (AR; number of non-extruded anthers) was scored on field in both panels 5–10 days post-anthesis for the years 2013, 2014, and 2015 by observing the anthers retained inside the four pairs of primary and secondary florets sampled from the central portion of four spikes per plot. For years 2015 and 2016, AR was scored by taking 10 spikes per plot, harvested 5–10 days post-anthesis and held at −20°C until analysis in the laboratory. Field and laboratory based AR data (2015) were compared to check the usefulness of both methods. Anther extrusion was calculated by subtracting the number of retained (non-extruded) anthers from 24 (since the total number of anthers housed by eight florets is 24) as described in Muqaddasi et al. (2016).

The mean AE in each growing year was used to calculate the best linear unbiased estimates (BLUEs), assuming fixed genotypic effects. The calculations were performed using GenStat v16 software (VSN International, Hemel Hempstead, Hertfordshire, UK) applying the “Mixed models REML” module and the “Linear mixed models” function. A Pearson product moment correlation among the growing years and BLUEs was also calculated. The repeatability among the replicates (R) and heritability (H²) among the years were calculated from the individual variance components as R=σG2σG2+(σe2nR) and H2=σG2σG2+(σG×E2nE)+(σe2nE×nR) where, σG2 represents the variance of genotype, σG×E2 is the variance of the interaction of genotype by environment, σe2 the variance of error, nE is the number of environments and nR the number of replicates (Nyquist and Baker, 1991).

All 207 wheat genotypes were analyzed using a 15K Infinium SNP array, which is an optimized and reduced version of the 90K iSELECT SNP-chip described by Wang et al. (2014). The development of the 15K SNP-chip and genotyping was performed by TraitGenetics GmbH (http://www.traitgenetics.com). The initial number of polymorphic loci was 13,006 for both panels. The SNP-chip data is available at digital object identifier (DOI) http://dx.doi.org/10.5447/IPK/2017/4. This DOI was created with e!DAL (Arend et al., 2014). A minor allele frequency threshold (MAF) of 5% was applied prior to the determination of marker-AE associations which reduced the size of the data matrices to 12, 066 × 111 for SP and 12, 191 × 96 for WP.

Genetic relationships among the SP and WP accessions were investigated graphically via principal component analysis (PCA) based on SNP genotypes in software R (Team, 2016). The first two principal components of each panel were depicted in two-dimensional space to show the clustering of accessions.

Assuming that the accession share some degree of relatedness, a genetic variance-covariance kinship matrix (K) as coefficient of co-ancestry (Lynch and Walsh, 1998) was derived from all SNP markers (MAF > 5%) between all pairs of accessions in each panel individually by using GenStat v16.

Genome-wide association scan was performed by following Yu et al. (2006). Marker trait associations (MTAs) involving AE and the SNP markers were identified using the “QTL analysis” module and the “Single trait association analysis” function implemented in GenStat v16, applying a variance-covariance kinship matrix (VCOV(Gi)=2Kσg2) for the purpose of correction for population stratification (Lynch and Walsh, 1998).

Assuming the genetic relatedness among accessions, the following mixed model for MTA analysis on BLUEs dataset with SNP markers was used: yi=μ+αxi+Gi+(Gi+ei), where, y_i is the AE's BLUE value for i^th genotype, μ is the overall intercept, α is the marker effect, x_i is the genetic predictor for i^th genotype, G_i is the effect of i^th genotype and e_i is the error for genotype i. The residual effect (G_i + e_i) was assumed to arise from a normal distribution as ~N(0, σ²) which imposed a genetic variance-covariance structure. Genetic relatedness (structure) among accessions was considered by assuming the random genotype effect Gi~N(0,2Kσg2). The threshold of |log₁₀(P)| > 3.0 was set in GenStat v16 for the detection of significant MTAs. The MTA analysis was performed separately for each panel based on BLUEs. Observed vs. expected |log₁₀(P)| value distributions were displayed in the form of quantile-quantile (qq) plots. For the set of markers unlinked to an AE-QTL, the |log₁₀(P)| values would be expected to be uniformly distributed along the diagonal, however, extremely large deviations are indicative of spurious associations (Yu et al., 2006; Kang et al., 2008). The additive contribution of each linked marker was calculated with GenStat v16: negative effects reflect a decrease in AE (unfavorable alleles), and positive ones an increase (favorable alleles). For the bi-allelic SNP markers the effects were related to the respective allele with minor allele frequency at each locus. TASSEL v3.0 (Bradbury et al., 2007) was used to calculate the phenotypic variances (R²) explained by individual markers.

By using a set of mapped SNP markers (SP = 10,578; WP = 10,649) genome-wide linkage disequilibrium (LD) was studied in both panels. Intra-chromosomal and genome-wide quantification and graphical representation of LD heat maps and intra-chromosomal LD decay were accomplished with GenStat v16 and R (Team, 2016) by plotting the squared correlation coefficients (r²) vs. map distances in centiMorgans (cM).

The 15K SNP Infinium array provided genome-wide polymorphic markers in high density. We performed the pairwise marker LD analysis for each chromosome individually for both panels. We set the significant LD regions by assuming significant marker to be in LD with other markers on the same chromosome by setting the threshold of r² ≥ 0.30. Thus, we extracted the flanking sequences of all the markers from wheat 90K database (Wang et al., 2014) associated with AE in significant LD regions (r² ≥ 0.30) and blasted (BLASTN; Altschul et al., 1990) against the wheat genome assembly (IWGSC1 + popseq) (Consortium, 2014) to find their corresponding genes and transcripts. The functional descriptions of these genes were retrieved from (ftp://ftpmips.helmholtz-muenchen.de/plants/wheat/IWGSC/genePrediction_v2.2/). All the markers sequences were blasted in Geneious v9.1.5 (Kearse et al., 2012).

Additionally, the orthologs of the corresponding genes of the close relatives of wheat i.e., Brachypodium distachyon, Sorghum bicolor, rice, and barley were retrieved from Wang et al. (2014) and barley sequence assembly RefSeq v1.0 (International Barley Sequencing Consortium pre-publication access). Furthermore, the wheat transcript IDs of all the genes in significant LD regions were given in http://www.wheat-expression.com to retrieve their expression profiles in above-ground plant tissues at log2 transcripts per million (tpm) scale (Borrill et al., 2016). Heat maps of expression profiles were created by using the R package gplots (Warnes et al., 2009) in R (Team, 2016). We portioned the expression data profiles of above ground tissues in two groups i.e., flowering and non-flowering tissue types and calculated the mean expression for all transcripts. One way repeated-measures (RM) analysis of variance (ANOVA) was used to determine significant difference between expression profiles of flowering and non-flowering tissue types. Differences in the means were determined by Holm-Sidak post-hoc test (Holm, 1979) by setting a statistical significance of P < 0.05. Statistical tests were done in SigmaPlot v13.0.

Anther extrusion (AE) was scored for 207 wheat accessions (111 spring and 96 winter) and showed a wide range of variation among accessions. The habit of growth, country of origin and AE scores of all 207 accessions are given in Table S1. For the year 2015, two methods to measure the AE in both spring (SP) and winter wheat (WP) panels were adopted, i.e., (1) direct scoring of the anthers on the field and (2) collecting and freezing the spikes and determining the AE in the laboratory. Both methods strongly correlated with one another (Pearson′sr(SP) = 0.81, P < 0.0001; r(WP) = 0.85, P < 0.0001) (Figure S1). The laboratory-based method was preferred for the subsequent genetic analysis, because it avoided difficulties related with the on field AE scoring and lodging of the non-dwarf accessions and because it allowed the processing of large numbers of accessions. The BLUEs of both SP and WP approximated a normal distribution and ranged from 2.94 to 20.92 (mean = 13.24) for SP and from 5.55 to 22.14 (mean = 14.44) for WP (Table S1). The observed range in AE was comparable in the 2013 and 2014 trials, while the 2015 and 2016 experiments produced higher mean values in both panels (Figure 1). The extent of the AE Pearson's correlation among growing years and BLUEs ranged from r = +0.53−0.91 in SP and r = +0.38−0.80 in WP (Figure 1). The repeatability levels in all years were >0.95 in both panels. The broad sense heritability was 0.84 across years for SP and 0.65 for WP. Six SP and 12 WP accessions showed AE >80% based on their BLUEs ranging from 19.33 to 20.92 (SP) and 19.65 to 22.14 (WP) (Table S1), and the trait was stable across years (Figure 2).

The imposition of minor allele frequency (MAF) threshold of 5% reduced the SNP marker set to 12,066 markers for SP and 12,191 for WP, of which 10,578, and 10,649 had known map locations, respectively (Wang et al., 2014), spreading unevenly across all three constituent genomes; however, each of the 21 chromosomes was represented (Figure S2). Almost half of the markers mapped to B-genomes in both panels (SP = 49.42%; WP = 49.33%) followed by A-genomes (SP = 39.08%; WP = 38.87%), while D-genomes had the lowest marker coverage (SP = 11.50%; WP = 11.80%).

Intra-chromosomal and genome-wide linkage disequilibrium (LD) assessment was performed using adjacent (pair-wise) SNP loci which revealed that LD decayed rapidly with increasing genetic map distances. The LD analyses of the chromosomes harboring the markers significantly associated with the AE (Figures S3, S4) showed a considerable variability in chromosome-wise and genome-wide analysis of LD. This suggested that the power of GWAS can be improved by increasing the polymorphic loci density, especially on the D-genome which is least represented among all genomes.

A PCA based test for the existence of clustering among accessions for both panels is illustrated in Figure 3. In order to simplify the PCA plots, we highlighted the accessions with different colors based on their continent of origin. A clear separation was found between the European and Asian accessions in SP, while there was absence of pronounced clustering for the American ones. For WP, no distinct grouping of sub-populations was observed.

Marker-trait associations (MTAs) were identified for AE by using a mixed linear model (Yu et al., 2006) correcting for population stratification by using a kinship matrix (Lynch and Walsh, 1998) among pairs of accessions. In total, 23 (SP = 11; WP = 12) significant MTAs (|log₁₀(P)| > 3.0) were identified in GWAS analyses (Figures 4, 5). Of the SP, nine significant markers had favorable additive effects, while two were unfavorable and the significant markers explained phenotypic variances (R²) ranging from 9.75 to 14.27%. The significant MTAs of SP were detected on chromosomes 2B, 3A, 3B, 6A, 7A, and 7B (Table 1 and Table S2). Of the WP significant markers, three markers had favorable while nine had unfavorable effects on AE and the markers explained phenotypic variances ranging from 9.44 to 16.98%. The WP significant MTAs were detected on chromosomes 1A, 1B, 2B, 6A, and 6B (Table 1 and Table S3). Chromosome 6A harbored MTAs in the genomic region (135–141 cM) shared by both panels. The quantile-quantile (qq) plots, which compare observed vs. expected |log₁₀(P)| values, indicated that the model we used was sufficiently stringent to control for false positives (Figures 4, 5).

Analysis of all the significant markers and the markers in significant LD regions (r² ≥ 0.30) revealed 26 and 52 distinct markers in SP and WP, respectively. The corresponding orthologous genes/regions to these SNP markers were retrieved from wheat, barley, Brachypodium distachyon, sorghum and rice. Tables S2, S3 document the transcripts found in the different species for SP and WP, respectively.

We extracted the wheat transcript IDs in significant LD regions and retrieved the gene expression data of above-ground tissues from http://www.wheat-expression.com. The expression data showed that some of these genes were expressed relatively higher in the tissues active at the anthesis stage i.e., pistils, stamens, spikes, and spikelets (Figures 6, 7). One way RM-ANOVA and all pair-wise multiple comparisons (Holm-Sidak method) on expression data (flowering vs. non-flowering tissues) suggested that the differences in the mean expression values between both groups are significant (P < 0.05) (Tables S4A,B). The genes and annotated functions of all the markers in significant LD region (r² ≥ 0.30) revealed putative genes which could be directly or indirectly involved in AE in both panels.

Anther extrusion (AE) was scored for a total of 207 wheat accessions including 111 spring (SP) and 96 winter (WP) over the period of 3 growing years which resulted in broad sense heritability values of 0.84 for SP and 0.65 for WP. Anther extrusion showed a strong phenotypic variation among accessions in both panels. The repeatability and broad sense heritability values for AE in both SP and WP complemented previous reports based on various wheat panels and different environments (Singh et al., 2007; Skinnes et al., 2010; Lu et al., 2013; Langer et al., 2014; Buerstmayr and Buerstmayr, 2015; Boeven et al., 2016). The continuous frequency distribution of AE suggested the possible contribution of several genes which was confirmed in the subsequent genetic analyses. These findings suggest that selection for AE could be performed to improve AE in wheat.

For efficient QTL detection, a high marker allele density increases the chances for the detection of linked QTL. We used 15K SNP array which guarantees the best polymorphic genetic markers taken from the 90K SNP chip reported by Wang et al. (2014). In our marker analysis, sub-genome-wise coverage was maximum for B-genome and least for D-genome in both spring and winter panels. A potential reason for the lowest marker coverage on D-genome is the lack of polymorphism reported for this genome (Akhunov et al., 2010). One reason may be the relative younger age (whereby less mutations could have accumulated) of the D-genome (ca. 1–2 million years) compared to A and B-genomes which diverged ca. 7 million years ago and which were potentially involved in its homoploid formation (Marcussen et al., 2014, but see also Li et al., 2015). Moreover, the D-genome hybrid origin itself would have been an intense bottleneck reducing genetic diversity. This least D-genome coverage is also seen by using other marker types, i.e., DArT and GBS markers (Muqaddasi et al., 2016, 2017). By using these resources, we detected several AE-QTL with moderate significances and effects. This result is in accordance with the studies about QTL mapping of AE (Skinnes et al., 2010; Lu et al., 2013; Buerstmayr and Buerstmayr, 2015; Muqaddasi et al., 2016). It implies that no single major gene is governing the floral architecture of wheat, but that the co-evolution of several genes resulted in the observed phenotypes. Nevertheless, one could imagine the pyramiding of several favorable AE-QTL in a variety.

Our linkage disequilibrium (LD) analysis resulted in extremely low r² values between the adjacent pairs of marker loci. This suggested that due to the absence of closely linked SNPs especially on the D-genome, some QTLs may possibly have not been detected. As fast LD decay is a precondition for genome-wide association mapping, the increase of marker density may improve the detection of QTL for AE.

Anther extrusion in our study behaved as a quantitative trait in both SP and WP; controlled by the concerted action of several loci and explaining a minor to modest level of phenotypic variance (R²(SP) = 9.75−14.24%; R²(WP) = 9.44−16.98%) for individual SNPs.

Our AM analysis revealed 11 significant SNPs in SP and 12 in WP, defining six (on chromosomes 2B, 3A, 3B, 6A, 7B, and 7D) and five (on chromosomes 1A, 1B, 2B, 6A, and 6B) QTL regions, respectively. Most studies published on AE are linkage-mapping studies (Skinnes et al., 2010; Lu et al., 2013; Buerstmayr and Buerstmayr, 2015) and used amplified fragment length polymorphism (AFLP), diversity array technology (DArT) and simple sequence repeat (SSR) makers. To our knowledge, AM studies published on the subject comprise of AM in spring wheat using mainly DArT and GBS markers (Muqaddasi et al., 2016, 2017) and SNP markers in winter wheat (Boeven et al., 2016). The use of dissimilar marker systems and/or background of mapping populations both in linkage and association mapping studies on AE make it challenging to compare the QTL intervals. However, comparison of QTL with previously published key QTL (imparting more than 10% of phenotypic variance) in the similar regions of chromosomes can give a better hint about the nature of AE and novelty of QTL.

For example, in our study chromosome 1A harbored a QTL (113.19–115.72 cM), represented by significant SNP “JD_c285_1438” in WP and explaining 14.51% of phenotypic variation (R²). Skinnes et al. (2010) detected a major QTL for AE on chromosome 1A in a bi-parental doubled haploid population which was located in the confidence interval between 92 and 102 cM (R² = 18.3%). Similarly a WP-QTL on chromosome 6B (71.18–72.29 cM) imparted an average of 11.73% of phenotypic variation. Buerstmayr and Buerstmayr (2015) detected QTL for anther retention, an equivalent of AE, in a bi-parental mapping population on chromosome 6BL, with an interval covering a distance of 25 cM (~70–90 cM) and explaining a mean phenotypic variance of 11.3%. Based on the similar QTL regions in these two examples we assume they might constitute similar QTL, although the marker systems and populations are different. To our knowledge, no other AE-QTL explaining more than 10% of phenotypic variance has been reported in the similar regions to our QTLs and therefore these QTL can be considered novel.

Nevertheless, it is important to keep in mind that all recent reports suggest an extremely complex nature of AE which warrants that the AE is controlled by many loci, mostly small effect, contrary to what was initially reported that AE is controlled by very small number of genes (Sage and De Isturiz, 1974). Moreover, as most QTL are non-overlapping in earlier studies as well, it becomes apparent that AE greatly depends on the mapping population background and the marker system used. Moreover, our GWAS analyses yielded 70.76 and 63.94% of the total phenotypic variance imparted by SP and WP QTL, respectively. These results are analogous to all previous reports where a large part of the phenotypic variance remained unexplained which shows that marker density could be improved to capture the small effect loci and unexplained phenotypic variance.

We found several orthologous genes in close wheat relatives i.e., barley, Brachypodium distachyon, rice and sorghum. The annotated functions of several genes suggested that the proteins they encode can be important components of flowering and reproduction related traits, for example, D111/G-patch domain-containing protein (Quiapim et al., 2009), MADS-box (Theissen et al., 2000), nuclear pore anchor (Xu et al., 2007). Likewise, the expression profiles of most of these genes suggested a significantly higher expression in flowering tissues as compared to other non-flowering above-ground tissues. Based on annotated function and expression profiles, we suggest that these genes could either be the genes-of-interest or linked to the true gene-of-interest. However, further investigation and validation of the roles of these genes in wheat is still needed. These AE-associated genes can be considered as possible candidate genes and therefore provide a good resource for future map based cloning. Furthermore, additional genes in LD which were not present on the 15K SNP-chip may be the causal genes for AE-QTL.

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.