Genome-wide association analysis on pre-harvest sprouting resistance and grain color in U.S. winter wheat

Plant materials

A set of 185 winter wheat accessions [33] was assembled to include 130 hard winter wheat and 55 soft winter wheat accessions as listed in Additional file 1: Table S1. A mapping population of 155 F₆ recombinant inbred lines (RILs) derived from the cross of Tutoumai A x Siyang 936 [7, 34] was used to validate the SNPs that showed significant associations with the Qphs.hwwgr-4A.

Pre-harvest sprouting evaluation

In the greenhouse experiments, five plants per accession were grown in a 13 by 13 cm Dura-pot (Hummert Int. Topeka, KS) under the growth condition listed in Additional file 2: Table S2 after vernalization for 7 weeks at 4 °C in a cold chamber. The GWAS panel was evaluated for PHS in the greenhouse experiments of fall (August-December) 2011, spring (January-May) and fall 2012, and spring 2013. All experiments were conducted in a randomized complete block design with two replications of five plants.

The GWAS panel was also planted for PHS resistance evaluation in the Kansas State University Rocky Ford Wheat Research Farm, Manhattan, KS and the Agricultural Research Center-Hays, Hays, KS, respectively, in the summers of 2013 and 2014. About 30 seeds per accession were planted in a 1.22-m-long single-row plot, and each experiment had two replications.

When wheat plants reached physiological maturity, similar to Zadoks scale 91 [35], spikes that lost green color [36] were harvested from both greenhouse and field experiments, and evaluated for PHS in the lab. Five spikes per accession that were harvested from each replicate were air-dried for 5 d in a greenhouse, and then stored at -20 °C to maintain dormancy for PHS evaluation. After all accessions had been collected, the greenhouse- harvested spikes were air-dried 9 d and field-harvested spikes for 5 d at room temperature. The additional drying days were determined based on preliminary test results of randomly selected samples from field and greenhouse experiments that maximize phenotypic differences among genotypes. After the dried spikes had been immersed in de-ionized water for 12 h, they were enclosed in a moist chamber at 22 ± 1 °C with an attached humidifier that ran twice daily at 2 h each time to maintain high moisture in the chamber. After 7 d of incubation, the germinated and non-germinated kernels were hand-threshed and counted separately to calculate the percentage of germinated kernels from all five spikes of each replication.

Evaluation of grain color

Grain color was evaluated for grains harvested from one field experiment (2009-2010 Enid Oklahoma) and the fall 2011 greenhouse experiment at Manhattan KS. For each accession, ten seeds were soaked in 1 M sodium hydroxide (NaOH) for 1 h to increase the color contrast. Grain color intensity was determined visually using a scale of 1 to 4, where 1 represents white, 2 light red, 3 red and 4 dark red.

DNA isolation and genotyping

Leaf tissue was collected at the two-leaf stage, and genomic DNA was isolated using a modified cetyltrimethyl ammonium bromide (CTAB) method [33]. A total of 446 polymorphic SSR markers were selected to genotype the association panel based on PCR product quality, chromosome distribution in available genetic maps (http://wheat.pw.usda.gov/GG3/; verified 11 Aug. 2010), and previously reported associations with PHS resistance. One expression sequence tag (EST), ZXQ118 [12] and three gene markers of PM19A1 and PM19A2 [6] were used to determine the association between PHS resistance and Qphs.hwwgr-4A. Five sequence-tagged sites (STS) from three Tamyb10 genes [25] were analyzed to determine QTL for GC. Amplification, separation and scoring of polymorphic chain reaction (PCR) products followed Zhang et al. [33].

The GWAS panel was also genotyped with the Wheat 9K and 90K SNP arrays [37, 38] at USDA-ARS Cereal Crops Research Unit (Fargo, ND). SNPs with less than 5 % minor allele frequency (MAF) or with more than 15 % missing data were removed. A total of 5,921 and 21,600 SNPs were scored from the 9K and 90K SNP arrays, respectively. Association analysis was initially conducted using the 9K genotypic data, and 28 non-redundant SNPs with p < 0.001 were then selected and pooled together with the 90K data. Totally, 21,628 SNPs were used for the final analysis. Also, one SNP in the promoter region and two SNPs in the coding region of the TaPHS1 gene [3, 4] were analyzed using three Kompetitive Allele Specific PCR (KASP) assays. Sequences that harbored significant SNPs and SSR markers were searched against the W7984 reference sequence to estimate their putative chromosome positions.

Population structure and kinship

Population structure was characterized by a set of 1500 SNPs that are evenly distributed on all the 21 wheat chromosomes using the admixture model in STRUCTURE 2.3.4 [39]. K-values ran from 2 to 20 with 10 iterations set for each k-value. The burn-in time and replication number were set at 2 × 10⁵ and 2 × 10⁴, respectively. For each trait, Bayesian information criterion (BIC) [40] was applied to determine the optimum number of subpopulations. Marker-based kinship was estimated to approximate the probability of two individuals being identical by descent through adjusting the average probability of identical in the state between random individuals [41]. Kinship was calculated with the same set of 1,500 SNPs used for structure analysis using SPAGeDi package [42].

Statistical analysis and genome-wide association analysis

Best linear unbiased predictions (BLUPs) were calculated for each accession evaluated in the greenhouse and field experiments using the 'lme4' package in R 3.2.2 [43] with year and location as random effects in the model. Genome-wide association analysis was conducted using two models: the generalized linear model (GLM) with the Q matrix as fixed effects, and the mixed linear model (MLM) with a Q matrix as fixed effects and a kinship matrix as random effects. These two models were applied to each experiment for GC and PHS resistance, and model fitness was determined based on the BIC values. Association analysis of SNP data was conducted using the genome association and prediction integrated tool (GAPIT) implemented in R [44], and association analysis of SSR data was conducted using PROC MIXED procedure in SAS 9.3 (SAS Institute Inc., Cary, NC). A threshold of p < 0.001 was set to claim significant associations between SSR markers and the traits (GC and PHS resistance), and p < 0.0001 was set to claim significant associations between SNPs and the traits. Linkage disequilibrium and haplotype analyses of the significant SNPs were performed with HAPLOVIEW v.4.2 (http://www.broadinstitute.org/scientific-community/science/programs/medical-and-population-genetics/haploview/haploview).

QTL analysis

A linkage map covering the 4A QTL region was constructed for the RIL population of Tutoumai A x Siyang 936 using KASP markers converted from significant SNPs from the association study, and previously mapped SSR markers in Liu et al. [45] and GBS-SNPs in Lin et al. [7] by JoinMap version 4.0 [46]. Recombination fractions were converted into centiMorgans (cM) using the Kosambi function [47]. Interval mapping (IM) using sprouting data from the 2005 and 2006 greenhouse experiments and their combined mean was performed using WinQTLCart 2.5 [48]. LOD thresholds to claim significant QTL for each dataset were determined from 1000 permutations [49].

Phenotypic variations in grain color and pre-harvest sprouting

Twenty-nine accessions were scored as white wheats, and 156 accessions as red wheats. The GC scores were highly consistent between greenhouse and field grown seeds (Fig. 1) with a high correlation coefficient of 0.87 (P < 0.0001), indicating a low genotype-by-environment interaction for GC.

Significant correlations for sprouting rates were observed among most of the eight experiments (Table 1). Cluster analysis showed high similarities in sprouting rates of accessions among all the field experiments, but significant differences between the field and the greenhouse experiments (Fig. 2a). The broad sense heritability across all eight experiments was high (0.83), with 0.62 in the greenhouse experiments and 0.92 in the field experiments. The population could be roughly divided into three subgroups (Fig. 2b), with average sprouting rates of 13.9 % in Group 1, 35.5 % in Group 2 and 60.3 % in Group 3, and average GC scores of 3.0 in Group 1, 2.6 in Group 2 and 1.7 in Group 3, indicating red wheats were more likely to have low sprouting rates. Most of the soft winter wheats were clustered to Group 1, as well as some hard white winter (HWW) wheat accessions from the Regional Germplasm Observation Nursery (RGON). The rest accessions from the RGON were mostly clustered to Group 2, whereas accessions from the Southern Regional Performance Nursery (SRPN) and the Northern Regional Performance Nursery (NRPN) were mainly clustered to Group 2 and Group 3.

Corr Coeff	2011F	2012S	2012F	2013S	2013_MH	2013_Hays	2014_MH	2014_Hays	GH_BLUP
2012S	0.238***
2012F	0.468***	0.313***
2013S	0.167*	0.543***	0.337***
2013_MH	0.276***	0.171*	0.212**	0.153*
2013_Hays	0.407***	0.227**	0.402***	0.237***	0.741***
2014_MH	0.384***	0.312***	0.576***	0.242***	0.686***	0.721***
2014_Hays	0.358***	0.228**	0.431***	0.266***	0.747***	0.755***	0.821***
GH_BLUP	0.559***	0.718***	0.821***	0.710***	0.268***	0.437***	0.551***	0.448***
Field_BLUP	0.399***	0.265***	0.463***	0.253***	0.869***	0.890***	0.909***	0.929***	0.484***

*p < 0.05; **p < 0.01, ***p < 0.001

Genome-wide association studies on grain color

The association mapping population can be classified into eight genotypic groups based on the allele combinations of the gene markers of Tamyb10 genes on group 3 chromosomes. The average GC scores in each group tended to increase as the number of red color alleles increases. However, red wheat accessions T154, LA02-923, MO040192 and NC04-15533 do not contain the red alleles (abd) at any of the three loci, whereas white accessions KS05HW15-2 and OK06848W carry the red allele of Tamyb10-A1 (Abd), and white accessions KS05HW136-3 and CO03W139 carry the red allele of Tamyb10-D1 (abD) (Fig. 3), suggesting that other genes besides Tamyb10 may also contribute to GC, or the markers for Tamyb10 genes may not be diagnostic in some genetic backgrounds.

Genome-wide association studies on pre-harvest sprouting resistance

Relationships between grain color and pre-harvest sprouting resistance

Validation of the significant SNPs for the 4A QTL in a bi-parental population

Seventeen KASP assays were designed based on the sequences of the significant SNPs identified in the 4A QTL region for PHS resistance. Four of the KASP markers (Additional file 4: Table S4) showed co-segregation among the F₆ RILs of 'Tutoumai A' × 'Siyang 936', and were mapped between the two previously reported flanking GBS SNPs (GBS212432, GBS109947) for the QTL [7] at 1.02 cM to GBS212432 and 2.10 cM to GBS109947 (Fig. 4). These four SNPs showed the highest LOD scores in all experiments, and explained up to 31.76 % of the phenotypic variance in the population.

Linkage disequilibrium

Linkage disequilibrium (LD) parameter r ² was calculated to determine the linkage relationship between SNPs from different QTL and link the markers with unknown positions to known QTL. LD was calculated for the 125 SNPs that were significantly linked to nine PHS resistance QTL in at least two experiments. Strong LD was detected for SNPs within each PHS resistance QTL region, but not between different QTL (Fig. 5b), indicating that those QTL for PHS resistance were independent. Pair-wise r ² values were also estimated for the 17 SNPs that were tightly linked to the five GC QTL. Similarly, strong LD was not detected among the SNPs linked to GC QTL in group 3 chromosomes (Fig. 5a).

Genetic positions of most significant SNPs for GC on chromosome 3D and Tamyb10-D1 could not be determined using the W7984 reference sequence, and SNPs significantly related to GC on chromosome 3D, D_GA8KES402JVT1Y_74 and BS00067163_51, were far apart from each other on the chromosome 3D. However, LD analysis suggested that these SNPs were tightly linked to Tamyb10-D1, and thus they linked to the same 3D QTL for GC (Fig. 5a).

QTL for grain color

Wheat GC has been a classic example for dissection of a quantitative trait [50] and three genes on wheat chromosomes group 3 have long been proposed as the genes controlling wheat GC. Several previous studies have mapped the three genes as major QTL as well as some minor QTL on chromosomes 2B, 2D, 5A and 6B for GC [1, 25, 51]. Being the first association study for wheat GC, we not only validated the effects of these three GC genes, Tamyb10-A1, Tamyb10-B1 and Tamyb10-D1, on the long arms of chromosomes group 3, but also identified a new QTL on the chromosome 1B for GC, suggesting that QTL on other chromosomes than these well-known QTL on chromosomes group 3 may also play a role in regulating GC in some wheat germplasm lines.

Groos et al. [1] mapped all group 3 QTL in a bi-parental population, but they did not discuss their effects of each QTL. In this study, a diverse association panel makes it possible to compare the effects of all the three QTL. Among the three genes on the chromosome group 3, Tamyb10-D1 had the largest effect on GC (R² = 0.24) and Tamyb10-A1 the smallest (R² = 0.11) in the association mapping panel whereas their minor allele frequencies (MAF) were similar (Table 2), indicating that the large effect of Tamyb10-D1 was not due to a higher MAF than other two genes. On the other hand, one single gene changed GC from white to red, and adding one or two additional GC genes only slightly increased redness (Fig. 3). Besides, QTL on chromosome 1B also contribute to GC, which was not reported previously, thus it is likely a new QTL for GC. That the red allele of the 1B QTL presents in the four red wheat accessions that do not carry the red alleles (abd) at any of the three Tamyb10 genes supports this assumption. Therefore, when breeding for white wheat cultivars, breeders not only need to remove the three Tamyb10 genes, but also should watch for other genes that may contribute to GC.

In this study, wheat GC was visually scored after increasing color intensities using sodium hydroxide solution. High repeatability in GC between the greenhouse and field experiments (Fig. 1) indicates that the GC scoring method used in the experiments is highly repeatable. All of the four QTL identified for GC were detected in both experiments, which provided genetic evidence that QTL for GC are relatively stable across environments.

QTL for pre-harvest sprouting resistance

QTL for PHS resistance have been mapped on almost all wheat chromosomes in previous bi-parental mapping studies. Although association studies on PHS resistance have been conducted using several types of markers [11, 30], the current study is the first report to use high density SNPs for GWAS on PHS resistance. We identified 12 QTL that were significance in at least two experiments.

For the QTL on 3AS, the causal gene (TaPHS1) has been cloned [3, 4]. One of the reported functional SNPs in the coding region [4] was significant in one greenhouse (fall 2011), whereas the functional SNP in the promoter region [3] was not significant in any of the experiments (Table 3; Additional file 3: Table S3). However, the most significant markers linked to the 3AS PHS resistance QTL were not the functional SNPs (Additional file 3: Table S3), which was probably due to environmental effects on phenotyping [3]. We had similar result for markers in the 4A QTL region. Only one of the candidate gene markers of PM19A1 was significantly associated with PHS resistance in the fall 2012 experiment, although the 4A QTL showed an extremely large effect on PHS in that experiment (Table 3; Additional file 3: Table S3). This was probably due to the fact that the gene expression was affected by environments or the gene markers are not diagnostic.

The QTL identified at the distal end of chromosome 3DS was not reported previously. LD analysis indicated that Qgc.hwwgr-3DS is different but could be a homeologue of Qgc.hwwgr-3AS (Fig. 5b). For the QTL on chromosome 3B, the sequences of the linked SSR markers are not found in the W7984 reference sequence, thus we cannot determine whether or not the significant SSR markers and SNPs on 3B linked to the same QTL. Similarly, we cannot determine the QTL positions on chromosome 7A.

QTL identified on chromosome 1A could be the same QTL reported by Knox et al. [52] in durum because Xwmc183 was located near the QTL region mapped in our study based on the W7984 reference sequence. The QTL on chromosome 2D is the same QTL as QPhs.ccsu-2D.4 [53] because of the common SSR Xgwm539. However, we cannot determine whether the QTL that were identified on chromosomes 1D, 5A, 5B, 6A and 6B were the same QTL reported in previous studies [1, 9, 30, 51, 54] due to the lack of common markers.

Variation of PHS resistance across environments

PHS is a complicated trait affected by many factors, including seed dormancy (SD) [15, 55–57], GC [1, 58], spike morphology, as well as environmental factors such as temperature, moisture and photoperiod after flowering [59, 60]. In the current study, PHS resistance of the tested accessions and QTL effects varied across environments with more variation observed among the greenhouse experiments than that among the field experiments (Fig. 2a). A total of four greenhouse experiments were conducted in the fall greenhouse cycles of 2011 and 2012 with the harvest time in winter, and the spring cycles of 2012 and 2013 with the harvest time in summer. The two seasons were highly different in growing and post-harvesting temperatures, which has been shown to influence PHS resistance [3, 6]. Meanwhile, in the field experiments at Manhattan and Hays, dry hot winds shortened maturity period, which greatly reduced environment effects on wheat PHS resistance. Therefore, PHS resistance was similar in the four field experiments.

Qphs.hwwgr-3AS and Qphs.hwwgr-4A were the major QTL for PHS resistance, and most frequently identified in all experiments. However, Qphs.hwwgr-3AS was detected more frequently in the field experiments, while Qphs.hwwgr-4A was detected more frequently in the greenhouse conditions (Table 3; Additional file 3: Table S3), which might be due to high temperatures in field conditions during late grain maturation that suppressed the expression of Qphs.hwwgr-4A [6].

According to the heat map derived from individual PHS ratings across all the experiments, the population can be roughly divided into three clusters (Fig. 2b). Most of the soft winter wheats had low germination rates, and were clustered to Group 1. Wheat cultivars from RGON were mostly clustered to Group 1 and Group 2, whereas accessions in SRPN and NRPN showed higher germination rates, and were mainly clustered to Group 2 and Group 3. These results indicated that the soft winter wheat accessions grown in the humid climate during harvest season had a higher selective pressure on PHS resistance than the hard winter wheat accessions from the Great Plains that grown under relatively drier climate.

Validation of the markers for the QTL on 4A

In this study, a RIL population from ‘Tutoumai A’ x ‘Siyang 936’ was used to validate the position of significant SNPs for 4A PHS resistance QTL. Four polymorphic SNPs from GWAS were successfully mapped to the QTL region, and they are more closely linked to PHS resistance than previously reported flanking markers, GBS212432 and GBS109947, for this QTL [7]. This result indicates that GWAS provides more power to increase marker density and mapping resolution, whereas bi-parental populations can further validate the positions of new markers. Barrero et al. [6] proposed PM19A1 and PM19A2 as the candidate genes for the 4A QTL and identified causal deletions in PM19A1 and PM19A2. We analyzed the markers developed based on the causal variation in Tutoumai A and Siyang 936, but did not find any polymorphism between the two parents. Therefore, a different gene or different causal SNP in the gene may control the PHS resistance of 4A QTL in this population, which was also supported by the results from the GWAS that the candidate gene markers contributed much lower phenotypic variation for PHS resistance than three other SNP markers (Ex_c66324_1151, wsnp_Ex_c13031_20625900, wsnp_Ex_rep_c66324_64493429) (Table 3; Additional file 3: Table S3).

Effect of grain color QTL on pre-harvest sprouting resistance

GC has been considered as an important factor for PHS resistance, and previous studies showed that seed dormancy level of a white-grained wheat line was improved by the introgression of an R gene [16]. In the current study, GC explained 26 % to 44 % of the phenotypic variance for PHS resistance, and Tamyb10-A1 and Tamyb10-D1 showed significant effects on both GC and PHS resistance, which agree with a previous study [1]. Tamyb10 genes encode R2R3-type MBY transcription factors, which regulate the accumulation of PA in the biosynthesis pathways [25]. Therefore, it is possible that these transcription factors showed pleiotropic effects by regulating more than one metabolism pathway, and had effects on improving wheat PHS resistance. However, the GC gene on 3BL, Tamyb10-B1, did not show any effect on PHS resistance in this study (Table 2; Table 6).

In this study, GC was significantly related to PHS resistance in field experiments, but had barely any effect in the greenhouse experiments (Table 5). Also, the Tamyb10-A1 gene affected PHS resistance in all of the four field experiments, and the Tamyb10-D1 gene only affected PHS resistance in the 2013 experiments. Such results suggested that environmental factors could be important triggers of pleiotropic effects of the GC genes on PHS resistance. That Tamyb10-B1 did not show any effect on PHS resistance might be due to the field environments of this study that could not trigger the expression of pleiotropic effect of the gene.

Although some GC genes contributed to wheat PHS resistance, many QTL for PHS resistance did not affect GC. Therefore, some red wheats can be highly susceptible to PHS, while some white wheats can be highly resistant [61, 62]. Breeding for PHS resistance, attention should be paid to these QTL with a major effect on PHS in most environments without a pleiotropic effect on GC, such as these on 3AS and 4AL. Pyramiding several of these genes in one cultivar should be able to avoid PHS damage in U.S. HWW.