Genetic Architecture of Anther Extrusion in Spring and Winter Wheat

Hybrid wheat breeding is gaining prominence worldwide because it ensures higher and more static yield than conventionally bred varieties. The cleistogamous floral architecture of wheat (Triticum aestivum L.) impedes anthers inside the floret, making it largely an inbreeder. For hybrid seed production, high anther extrusion is needed to promote cross pollination and to ensure a high level of pollen availability for the seed plant. This study, therefore, aimed at the genetic dissection of anther extrusion (AE) in panels of spring (SP), and winter wheat (WP) accessions by genome wide association studies (GWAS). We performed GWAS to identify the SNP markers potentially linked with AE in each panel separately. Phenotypic data were collected for 3 years for each panel. The average levels of Pearson's correlation (r) among all years and their best linear unbiased estimates (BLUEs) within both panels were high (r(SP) = 0.75, P < 0.0001;r(WP) = 0.72, P < 0.0001). Genotypic data (with minimum of 0.05 minor allele frequency applied) included 12,066 and 12,191 SNP markers for SP and WP, respectively. Both genotypes and environment influenced the magnitude of AE. In total, 23 significant (|log₁₀(P)| > 3.0) marker trait associations (MTAs) were detected (SP = 11; WP = 12). Anther extrusion behaved as a complex trait with significant markers having either favorable or unfavorable additive effects and imparting minor to moderate levels of phenotypic variance (R² (SP) = 9.75−14.24%; R² (WP) = 9.44−16.98%). All mapped significant markers as well as the markers within their significant linkage disequilibrium (r² ≥ 0.30) regions were blasted against wheat genome assembly (IWGSC1+popseq) to find the corresponding genes and their high confidence descriptions were retrieved. These genes and their orthologs in Hordeum vulgare, Brachypodium distachyon, Oryza sativa, and Sorghum bicolor revealed syntenic genomic regions potentially involved in flowering-related traits. Moreover, the expression data of these genes suggested potential candidates for AE. Our results suggest that the use of significant markers can help to introduce AE in high yielding varieties to increase cross fertilization rates and improve hybrid-seed production in wheat.

Introduction

In many crops, hybrids have almost completely replaced the conventionally bred varieties because of their potential of higher and more stable yields even in marginal environments (Schnable and Springer, 2013). Heterosis in F1 hybrids ensures greater economic yield and static yield stability (Longin et al., 2012; Muhleisen et al., 2014; Zhao et al., 2015). Although, hybrid wheat breeding has a long history, still the wheat hybrid market is very small. This slow progress in hybrid wheat seed production is due to the hermaphroditic cleistogamous architecture of wheat florets which impedes proper anther extrusion to allow sufficient pollen shedding outside the floret. To guarantee a sound hybrid wheat seed production program, the cleistogamous nature of wheat needs an alteration to open flowering (Whitford et al., 2013). This floral modification will allow anthers to extrude during anthesis for sufficient pollen shedding outside the florets, so they may become available for the female lines via wind pollination. Moreover, this floral adjustment will not only profit male lines but female lines as well, since they ought to have open flowers to receive sufficient amount of pollen with receptive stigmas. Subsequently, the existing insignificant cross fertilization rates in wheat can be increased by ensuring a higher level of anther extrusion for appropriate pollen shedding outside of florets.

Male sterility is the most efficient way to ensure the cross pollination in inbreeding crops and this can be achieved in wheat by several means, e.g., cytoplasmic male sterility, genic male sterility and the use of chemical hybridization agents (CHAs) (Perez-Prat and van Lookeren Campagne, 2002; Chen and Liu, 2014; Kempe et al., 2014). The absence of an efficient cytoplasmic male sterility system and toxic effects coupled with CHAs are major reasons for less hybrid wheat seed production (Longin et al., 2012; Whitford et al., 2013). However, currently, the most widely used industrial scale hybrid wheat production system is the induction of male sterility by means of CHAs (Longin et al., 2012). CHAs, applied exclusively on female-parents, render the plants male-sterile. These male-sterile lines (females) receive pollen from neighboring male lines. Nevertheless, besides achieving an efficient male-sterility system and the determination of heterotic groups (Zhao et al., 2015), successful development of hybrid wheat varieties depends upon the modification of floral architecture to ensure high cross-pollination capacity in wheat (Whitford et al., 2013). This can be guaranteed by adequate volume of anther extrusion and viable pollen-shed outside the florets.

Anther extrusion is a complex trait controlled by many genes. The genetic basis of AE has mostly been unveiled via conventional mapping studies (Skinnes et al., 2010; Lu et al., 2013; Buerstmayr and Buerstmayr, 2015). Association mapping (AM) is a more recent method to determine the quantitative trait loci (QTLs) in a non-related germplasm by marker-trait associations (Risch and Merikangas, 1996; Zhu et al., 2008). This method is advantageous over the conventional QTL mapping technique because it profits from the strength of linkage disequilibrium (LD) and higher recombination frequencies available in diverse germplasm and offers an increased mapping resolution, a smaller investment in research time and cost, and an opportunity to survey a wider range of alleles at any given locus (Zhu et al., 2008; Hamblin et al., 2011). Recently, AM studies to unveil the genetic architecture of AE in spring wheat by using mainly DArT and GBS markers (Muqaddasi et al., 2016, 2017) and winter wheat (Boeven et al., 2016) by SNP markers, also suggested a highly quantitative nature of AE.

In this investigation, genome wide association studies were carried out to identify the genetic basis of AE in diverse panels of spring and winter wheat accessions. The 15K SNP marker platform was chosen for its high density and highly informative and genome-wide polymorphic loci—a perquisite for efficient AM studies.

Materials and Methods

Plant Material and Phenotyping of Anther Extrusion

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).

Statistical Analysis

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=\frac{{\sigma }_G^2}{{\sigma }_G^2+\left(\frac{{\sigma }_e^2}{nR}\right)}$ and $H^2=\frac{{\sigma }_G^2}{{\sigma }_G^2+(\frac{{\sigma }_{G\times E}^2}{nE})+(\frac{{\sigma }_e^2}{nE\times nR})}$ where, ${\sigma }_G^2$ represents the variance of genotype, ${\sigma }_{G\times E}^2$ is the variance of the interaction of genotype by environment, ${\sigma }_e^2$ the variance of error, nE is the number of environments and nR the number of replicates (Nyquist and Baker, 1991).

Genotyping, Population Structure, and Kinship Estimates

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.

Marker Trait Association and Linkage Disequilibrium Assessment

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\left(G_i\right)=2K{\sigma }_g^2$) 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:

$y_i=\mu +\alpha x_i+G_i+(G_i+e_i),$

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 $G_i~N(0,2K{\sigma }_g^2)$. 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).

Analysis of Significant LD Regions and Synteny with Related Grasses

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.

Results

Phenotypic Data Analyses

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).

FIGURE 1

Figure 1. Boxplots and Pearson's product moment correlations among growing years and their Best Linear Unbiased Estimates (BLUEs) in spring and winter wheat panels. Numbers on the Y-axis represent the number of extruded anthers. xs indicate the outlier accessions.

FIGURE 2

Figure 2. Anther extrusion (AE) performance of the top anther extruding accessions based on Best Linear Unbiased Estimates (BLUEs) ±SE. Six spring and 12 winter wheat accessions showed AE >80%. The growth habit of the accessions is indicated within the bars. Numbers on the Y-axis represent the number of extruded anthers.

Marker Distribution, Linkage Disequilibrium, and Population Structure

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.

FIGURE 3

Figure 3. Principal component analysis (PCA) of spring and winter wheat panels based on SNP genotypes (MAF >5%). The spring wheat accessions bred in Europe and Asia form recognizable clusters, but those bred in the Americas do not. There is no pronounced clustering among winter wheat accessions. Color code is given in the figure.

Marker-Trait Associations for Anther Extrusion

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).

FIGURE 4

Figure 4. Summary of genome-wide association studies results of anther extrusion (AE) in spring wheat panel (SP). (A) Distribution of Best Linear Unbiased Estimates of AE (B) Manhattan plots based on mixed linear model using kinship matrix (K) to correct for population stratification. Horizontal dashed line corresponds to the threshold (|log₁₀(P)| > 3.0) for MTA estimation. (C) Quantile-quantile plots depicting expected vs. observed |log₁₀(P)| values.

FIGURE 5

Figure 5. Summary of genome-wide association studies results for anther extrusion (AE) in winter wheat panel (WP). (A) Distribution of Best Linear Unbiased Estimates of AE (B) Manhattan plots based on mixed linear model using kinship matrix (K) to correct for population stratification. Horizontal dashed line corresponds to the threshold (|log₁₀(P)| > 3.0) for MTA estimation. (C) Quantile-quantile plots depicting expected vs. observed |log₁₀(P)| values.

TABLE 1

Table 1. List of significant markers (−log₁₀(P) > 3.0) for anther extrusion in spring and winter wheat panels based on their BLUEs values.

Study of Genes in Significant LD Regions

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.

FIGURE 6

Figure 6. Expression of genes from spring wheat panel in above ground tissues. Expression values are given at log2 scale of transcripts per million (tpm): red, high expression; yellow, moderate expression; blue, low expression. Each individual gene (transcript_ID) is represented as horizontal row and different tissues are described in vertical columns. Tissues from pistil to spikelets show higher expression values for most genes compared to other tissues. Gene expression data were obtained from http://www.wheat-expression.com.

FIGURE 7

Figure 7. Expression of genes from winter wheat panel in above ground tissues. Expression values are given at log2 scale of transcripts per million (tpm): red, high expression; yellow, moderate expression; blue, low expression. Each individual gene (transcript_ID) is represented as horizontal row and different tissues are described in vertical columns. Tissues from pistil to spikelets show higher expression values for most genes compared to other tissues. Gene expression data were obtained from http://www.wheat-expression.com.

Discussion

Anther Extrusion is a Heritable and Phenotypically Diverse Trait

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.

SNP Density and Extent of LD Are Important Factors for AE Mapping

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 is Controlled by Concerted Action of Several Genetic Loci

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.

Synteny and Expression Data of Significant LD Regions Reveal Putative AE Genes

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.

Conclusions

Hybrid wheat breeding is a promising strategy to break the yield barriers (Muhleisen et al., 2014; Zhao et al., 2015). Limited cross-pollination resulting from the cleistogamous nature of wheat is a major limitation to sufficient seed set for efficient hybrid production. High AE can ensure higher cross-fertilization rates for wheat hybrid formation. Our results show that AE is a stable trait, has high heritability rates and could be improved by breeding. Moreover, our GWAS results demonstrate that AE is a complex trait controlled by the combined action of several genetic loci. The favorable markers could be used to incorporate AE in high yielding varieties via marker assisted breeding to promote hybrid-heterosis in the key crop wheat.

Author Contributions

QM performed the data and genome-wide association analyses and prepared the manuscript. JB participated in data analysis. AB, KP, and MR participated in the experimental design. MR conceived the idea and participated in the interpretation of results and preparation of manuscript. All authors read and approved the final manuscript.

Conflict of Interest Statement

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.

Acknowledgments

QM thanks DAAD (Deutscher Akademischer Austauschdienst) for supporting his Ph.D. candidature. We acknowledge Ellen Weiß for contribution in data collection. We are grateful to Guozheng Liu and Armin Meister for advices in statistical analysis. The presented barley transcript results were obtained by support of the International Barley Genome Sequencing Consortium (IBSC) through providing pre-publication access to the annotated barley reference sequence RefSeq v1.0 (downloadable after registration).

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017.00754/full#supplementary-material

References

Akhunov, E. D., Akhunova, A. R., Anderson, O. D., Anderson, J. A., Blake, N., Clegg, M. T., et al. (2010). Nucleotide diversity maps reveal variation in diversity among wheat genomes and chromosomes. BMC Genomics 11:702. doi: 10.1186/1471-2164-11-702

PubMed Abstract | CrossRef Full Text | Google Scholar

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2

CrossRef Full Text | Google Scholar

Arend, D., Lange, M., Chen, J., Colmsee, C., Flemming, S., Hecht, D., et al. (2014). e! DAL-a framework to store, share and publish research data. BMC Bioinformatics 15:214. doi: 10.1186/1471-2105-15-214

CrossRef Full Text | Google Scholar

Boeven, P. H., Longin, C. F., Leiser, W. L., Kollers, S., Ebmeyer, E., and Würschum, T. (2016). Genetic architecture of male floral traits required for hybrid wheat breeding. Theor. Appl. Genet. 129, 2343–2357. doi: 10.1007/s00122-016-2771-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Borrill, P., Ramirez-Gonzalez, R., and Uauy, C. (2016). expVIP: a customizable RNA-seq data analysis and visualization platform. Plant Physiol. 170, 2172–2186. doi: 10.1104/pp.15.01667

CrossRef Full Text | Google Scholar

Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635. doi: 10.1093/bioinformatics/btm308

CrossRef Full Text | Google Scholar

Buerstmayr, M., and Buerstmayr, H. (2015). Comparative mapping of quantitative trait loci for Fusarium head blight resistance and anther retention in the winter wheat population Capo x Arina. Theor. Appl. Genet. 128, 1519–1530. doi: 10.1007/s00122-015-2527-8

CrossRef Full Text | Google Scholar

Chen, L., and Liu, Y. G. (2014). Male sterility and fertility restoration in crops. Annu. Rev. Plant Biol. 65, 579–606. doi: 10.1146/annurev-arplant-050213-040119

PubMed Abstract | CrossRef Full Text | Google Scholar

Consortium, I. W. G. S. (2014). A chromosome-based draft sequence of the hexaploid bread wheat (Triticum aestivum) genome. Science 345:1251788. doi: 10.1126/science.1251788

CrossRef Full Text | Google Scholar

Hamblin, M. T., Buckler, E. S., and Jannink, J. L. (2011). Population genetics of genomics-based crop improvement methods. Trends Genet. 27, 98–106. doi: 10.1016/j.tig.2010.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scand. J. Stat. 6, 65–70.

Google Scholar

Kang, H. M., Zaitlen, N. A., Wade, C. M., Kirby, A., Heckerman, D., Daly, M. J., et al. (2008). Efficient control of population structure in model organism association mapping. Genetics 178, 1709–1723. doi: 10.1534/genetics.107.080101

PubMed Abstract | CrossRef Full Text | Google Scholar

Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., et al. (2012). Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649. doi: 10.1093/bioinformatics/bts199

PubMed Abstract | CrossRef Full Text | Google Scholar

Kempe, K., Rubtsova, M., and Gils, M. (2014). Split-gene system for hybrid wheat seed production. Proc. Natl. Acad. Sci. U.S.A. 111, 9097–9102. doi: 10.1073/pnas.1402836111

PubMed Abstract | CrossRef Full Text | Google Scholar

Langer, S. M., Longin, C. F. H., and Würschum, T. (2014). Phenotypic evaluation of floral and flowering traits with relevance for hybrid breeding in wheat (Triticum aestivum L.). Plant Breed. 133, 433–441. doi: 10.1111/pbr.12192

CrossRef Full Text | Google Scholar

Li, L. F., Liu, B., Olsen, K. M., and Wendel, J. F. (2015). A re-evaluation of the homoploid hybrid origin of Aegilops tauschii, the donor of the wheat D-subgenome. New Phytol. 208, 4–8. doi: 10.1111/nph.13294

PubMed Abstract | CrossRef Full Text | Google Scholar

Longin, C. F., Mühleisen, J., Maurer, H. P., Zhang, H., Gowda, M., and Reif, J. C. (2012). Hybrid breeding in autogamous cereals. Theor. Appl. Genet. 125, 1087–1096. doi: 10.1007/s00122-012-1967-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Q. X., Lillemo, M., Skinnes, H., He, X. Y., Shi, J. R., Ji, F., et al. (2013). Anther extrusion and plant height are associated with Type I resistance to Fusarium head blight in bread wheat line “Shanghai-3/Catbird”. Theor. Appl. Genet. 126, 317–334. doi: 10.1007/s00122-012-1981-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, M., and Walsh, B. (1998). Genetics and Analysis of Quantitative Traits. Sunderland, MA: Sinauer.

Google Scholar

Marcussen, T., Sandve, S. R., Heier, L., Spannagl, M., Pfeifer, M., Jakobsen, K. S., et al. (2014). Ancient hybridizations among the ancestral genomes of bread wheat. Science 345:1250092. doi: 10.1126/science.1250092

PubMed Abstract | CrossRef Full Text | Google Scholar

Mühleisen, J., Piepho, H. P., Maurer, H. P., Longin, C. F., and Reif, J. C. (2014). Yield stability of hybrids versus lines in wheat, barley, and triticale. Theor. Appl. Genet. 127, 309–316. doi: 10.1007/s00122-013-2219-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Muqaddasi, Q. H., Lohwasser, U., Nagel, M., Börner, A., Pillen, K., and Röder, M. S. (2016). Genome-wide association mapping of anther extrusion in hexaploid spring wheat. PLoS ONE 11:e0155494. doi: 10.1371/journal.pone.0155494

PubMed Abstract | CrossRef Full Text | Google Scholar

Muqaddasi, Q. H., Reif, J. C., Li, Z., Basnet, B. R., Dreisigacker, S., and Röder, M. S. (2017). Genome-wide association mapping and genome-wide prediction of anther extrusion in CIMMYT spring wheat. Euphytica 213:73. doi: 10.1007/s10681-017-1863-y

CrossRef Full Text | Google Scholar

Neumann, K., Kobiljski, B., Denčić, S., Varshney, R., and Börner, A. (2011). Genome-wide association mapping: a case study in bread wheat (Triticum aestivum L.). Mol. Breed. 27, 37–58. doi: 10.1007/s11032-010-9411-7

CrossRef Full Text | Google Scholar

Nyquist, W. E., and Baker, R. (1991). Estimation of heritability and prediction of selection response in plant populations. Crit. Rev. Plant Sci. 10, 235–322. doi: 10.1080/07352689109382313

CrossRef Full Text | Google Scholar

Perez-Prat, E., and van Lookeren Campagne, M. M. (2002). Hybrid seed production and the challenge of propagating male-sterile plants. Trends Plant Sci. 7, 199–203. doi: 10.1016/S1360-1385(02)02252-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Quiapim, A. C., Brito, M. S., Bernardes, L. A., Malavazi, I., DePaoli, H. C., Molfetta-Machado, J. B., et al. (2009). Analysis of the Nicotiana tabacum stigma/style transcriptome reveals gene expression differences between wet and dry stigma species. Plant Physiol. 149, 1211–1230. doi: 10.1104/pp.108.131573

PubMed Abstract | CrossRef Full Text | Google Scholar

Risch, N., and Merikangas, K. (1996). The future of genetic studies of complex human diseases. Science 273, 1516–1517. doi: 10.1126/science.273.5281.1516

PubMed Abstract | CrossRef Full Text | Google Scholar

Sage, G. C., and De Isturiz, M. J. (1974). The inheritance of anther extrusion in two spring wheat varieties. Theor. Appl. Genet. 45, 126–133. doi: 10.1007/BF00291142

PubMed Abstract | CrossRef Full Text | Google Scholar

Schnable, P. S., and Springer, N. M. (2013). Progress toward understanding heterosis in crop plants. Annu. Rev. Plant Biol. 64, 71–88. doi: 10.1146/annurev-arplant-042110-103827

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, S. K., Arun, B., and Joshi, A. K. (2007). Comparative evaluation of exotic and adapted germplasm of spring wheat for floral characteristics in the Indo-Gangetic Plains of northern India. Plant Breed. 126, 559–564. doi: 10.1111/j.1439-0523.2007.01396.x

CrossRef Full Text | Google Scholar

Skinnes, H., Semagn, K., Tarkegne, Y., Maroy, A. G., and Bjornstad, A. (2010). The inheritance of anther extrusion in hexaploid wheat and its relationship to Fusarium head blight resistance and deoxynivalenol content. Plant Breed. 129, 149–155. doi: 10.1111/j.1439-0523.2009.01731.x

CrossRef Full Text | Google Scholar

Team, R. C. (2016). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

Theissen, G., Becker, A., Di Rosa, A., Kanno, A., Kim, J. T., Münster, T., et al. (2000). A short history of MADS-box genes in plants. Plant Mol. Biol. 42, 115–149. doi: 10.1023/A:1006332105728

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S. C., Wong, D. B., Forrest, K., Allen, A., Chao, S. M., Huang, B. E., et al. (2014). Characterization of polyploid wheat genomic diversity using a high-density 90,000 single nucleotide polymorphism array. Plant Biotechnol. J. 12, 787–796. doi: 10.1111/pbi.12183

PubMed Abstract | CrossRef Full Text | Google Scholar

Warnes, G. R., Bolker, B., Bonebakker, L., Gentleman, R., Huber, W., Liaw, A., et al. (2009). Gplots: Various R Programming tools for Plotting Data. R Package Version 2(4).

Whitford, R., Fleury, D., Reif, J. C., Garcia, M., Okada, T., Korzun, V., et al. (2013). Hybrid breeding in wheat: technologies to improve hybrid wheat seed production. J. Exp. Bot. 64, 5411–5428. doi: 10.1093/jxb/ert333

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, X. M., Rose, A., Muthuswamy, S., Jeong, S. Y., Venkatakrishnan, S., Zhao, Q., et al. (2007). NUCLEAR PORE ANCHOR, the Arabidopsis homolog of Tpr/Mlp1/Mlp2/megator, is involved in mRNA export and SUMO homeostasis and affects diverse aspects of plant development. Plant Cell 19, 1537–1548. doi: 10.1105/tpc.106.049239

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Pressoir, G., Briggs, W. H., Vroh Bi, I., Yamasaki, M., Doebley, J. F., et al. (2006). A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat. Genet. 38, 203–208. doi: 10.1038/ng1702

CrossRef Full Text | Google Scholar

Zhao, Y., Li, Z., Liu, G., Jiang, Y., Maurer, H. P., Würschum, T., et al. (2015). Genome-based establishment of a high-yielding heterotic pattern for hybrid wheat breeding. Proc. Natl. Acad. Sci. U.S.A. 112, 15624–15629. doi: 10.1073/pnas.1514547112

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, C. S., Gore, M., Buckler, E. S., and Yu, J. M. (2008). Status and prospects of association mapping in plants. Plant Genome 1, 5–20. doi: 10.3835/plantgenome2008.02.0089

CrossRef Full Text | Google Scholar