Genomic Prediction and Genome-Wide Association Studies of Flour Yield and Alveograph Quality Traits Using Advanced Winter Wheat Breeding Material

Use of genetic markers and genomic prediction might improve genetic gain for quality traits in wheat breeding programs. Here, flour yield and Alveograph quality traits were inspected in 635 F6 winter wheat breeding lines from two breeding cycles. Genome-wide association studies revealed single nucleotide polymorphisms (SNPs) on chromosome 5D significantly associated with flour yield, Alveograph P (dough tenacity), and Alveograph W (dough strength). Additionally, SNPs on chromosome 1D were associated with Alveograph P and W, SNPs on chromosome 1B were associated with Alveograph P, and SNPs on chromosome 4A were associated with Alveograph L (dough extensibility). Predictive abilities based on genomic best linear unbiased prediction (GBLUP) models ranged from 0.50 for flour yield to 0.79 for Alveograph W based on a leave-one-out cross-validation strategy. Predictive abilities were negatively affected by smaller training set sizes, lower genetic relationship between lines in training and validation sets, and by genotype–environment (G×E) interactions. Bayesian Power Lasso models and genomic feature models resulted in similar or slightly improved predictions compared to GBLUP models. SNPs with the largest effects can be used for screening large numbers of lines in early generations in breeding programs to select lines that potentially have good quality traits. In later generations, genomic predictions might be used for a more accurate selection of high quality wheat lines.


Introduction
Baking quality of wheat is a complex trait controlled by many genes with minor effects and few genes with larger effects [1,2]. The amount and the composition of gluten proteins have large effects on baking quality of wheat. The major gluten loci are the high molecular weight glutenins (HMWGs) Glu-A1, Glu-B1, and Glu-D1 and the low molecular weight glutenins (LMWGs) Glu-A3, Glu-B3, and Glu-D3 [3,4]. Milling quality and water absorption are affected by the hardness of the grain. Grain hardness is, to a large extent, controlled by the Hardness locus on chromosome 5D, consisting of the genes Pina-D1, Pinb-D1, and Gsp-1 [5,6].
Breeding for improved wheat quality is challenging, because phenotyping of most quality traits requires laborious analyses of relatively large amounts of grain using expensive equipment. Baking tests can be used for evaluating the quality of wheat lines by determining bread loaf volume and texture. However, breeding programs typically do not have the resources to perform baking tests with

Plant Material
In total, 635 F 6 winter wheat lines from two breeding cycles in the Danish plant breeding company Nordic Seed A/S (Holeby, Denmark) were used in this study. The 321 lines of the first breeding cycle (set2014) were harvested in 2014, and the 314 lines of the second breeding cycle (set2015) were harvested in 2015. Six out of 96 crossing parents were used for both sets (years), while the remaining 90 crossing parents were used for one of the sets only. Each line was grown in an unreplicated 9.9 m 2 plot at Lolland, Denmark, following standard Danish agricultural practices. Approximately 180 kg of nitrogen were applied per hectare during the growth season, and no irrigation was used.

Phenotyping
Phenotyping was done at the Wheat Chemistry and Quality Laboratory at International Maize and Wheat Improvement Center (CIMMYT), Mexico. Grain samples were conditioned to 13.5% moisture content and then milled one time using a Brabender Quadrumat Jr. (Brabender GmbH & Co. KG, Duisburg, Germany). Flour yield was measured as the percentage of refined flour obtained from each grain sample after the bran fraction was sieved away through a 75 µm mesh sieve. A Chopin Alveograph (Tripette and Renaud, Villeneuve-la-Garenne, France) was used to obtain the Alveograph traits P, L, and W (Alveo P is dough tenacity, Alveo L is dough extensibility, and Alveo W is dough strength) using a modified version of the American Association of Cereal Chemists (AACC) method 54-30A [8,43]. Flour was mixed with a saltwater solution to form a dough, which was cut into discs. After resting 20 min at 25 • C, the dough discs were inflated with air, thus the dough expanded as a bubble. During inflation of each disc, a curve was recorded of pressure inside the bubble until it burst. Alveo P was the maximum height of the curve, Alveo L was the length of the curve, and Alveo W was the area under the curve. The coefficient of variation was calculated for each trait by dividing the standard deviation of the raw phenotypes with the mean.

Genotyping
DNA extraction was performed with a modified cetyl trimethylammonium bromide (CTAB) method [44] using leaves of three bulked, two-week-old seedlings for each line. Genotyping was done by TraitGenetics (Gatersleben, Germany) with the 15K Illumina Infinium iSelect HD Custom Genotyping BeadChip technology. The 13,006 called SNP markers were edited for minor allele frequency (MAF) > 1% and missing values < 10%, and the remaining 10,802 SNPs were used for the analyses. For each line, at least 90% of the SNPs were successfully genotyped.

Statistical Analysis
Genome-wide associations were studied using single marker regression. The following model was run for each of the 10,802 SNPs: where y is a vector of observed phenotypes, X and Z 1 are design matrices, b is a vector of fixed effects (mean and year/set), w i is the vector of genotypes of the i th SNP coded as 1, 0, −1, a i is the additive genetic effect of the i th SNP, u is a vector of additive genetic effects of the lines (u~N(0,G 1 σ g 2 ), where G 1 is a G-matrix (genomic relationship matrix) and σ g 2 is additive genetic variance), and e is a vector of random residual effects (e~N(0,Iσ e 2 ), where I is an identity matrix and σ e 2 is the residual variance).
Model effects and variance components were estimated by restricted maximum likelihood using the software package DMU [45]. For each chromosome, a G-matrix was calculated based only on the SNPs mapped to the remaining chromosomes and then used to correct for structure when analyzing the SNPs mapped to the excluded chromosome. This was done so that the SNP effect was not included twice in the model. G-matrices were calculated using the first method proposed by Van Raden [27]: where p i is the MAF of i th marker, Z 2 = M − P, M is a matrix with the marker alleles coded as 1, 0, -1, and P is a matrix where the i th column contains the MAF of SNP i calculated as 2(p i -0.5). Missing genotypes were set to 0 in matrix Z 2 . Genomic inflation factors, λ IF , were calculated for each trait by dividing the observed median value of the chi-squared statistic for the SNPs with the expected median value [46]. The inflation factor λ IF was used to correct the p-values for inflation by dividing the chi-squared statistic with λ IF and then re-calculating the p-values. The significance threshold was set using a Bonferroni correction: 5% divided by number of SNPs (0.05/10,802 = 4.6 × 10 −6 ).
DNA sequences surrounding the significantly associated SNPs were blasted against the annotated reference genome of the bread wheat variety Chinese Spring, IWGSC RefSeq v1.0 [47], using the BLAST tool of EnsemblPlants [48].
A Bayesian Power Lasso model where all SNPs were fitted simultaneously was also used for genome-wide association analyses in the Bayz software [49]: where y is a vector of observed phenotypes, b is a vector of the mean + year/set effect with design matrix X, Z 3 is a matrix of the alleles of the SNPs coded as 0, 1, 2, u is a vector of additive genetic SNP effects, and e is a vector of residual effects. The prior distribution of SNP effects was specified as an exponential power distribution [50]: where m is the number of markers, and β is shape parameter to control the sparsity, which affects the shrinkage of the SNP effects. Setting β to 1 makes the model equivalent to a standard Bayesian Lasso model. If β is set to less than 1, the difference between large and small marker effects can be increased [50]. Models with β of 0.2, 0.4, 0.8, and 1.0 were run, and the Deviance Information Criterion was used to determine the optimal β for each trait [51]. Residual effects were assigned a normal prior distribution. The residual variance, the mean, the year/set effect, and the rate parameter, λ RP , were assigned flat prior distributions. Model parameters were estimated using Markov Chain Monte Carlo with a length of 100,000 with 30,000 cycles as burn-in. The tool pbayz supplied with Bayz was used to compute posterior means, and the R package CODA was used to check for convergence to the posterior distribution [52]. Genomic predictions based on all 10,802 SNPs were conducted using the Bayesian Power Lasso model (3) and using a GBLUP model: where y is a vector of observed phenotypes, X and Z 4 are design matrices, b is a vector of fixed effect (mean and year/set), u is a vector of additive genetic effects (u~N(0 G 2 σ 2 g ), where G 2 is a G-matrix computed as above (2) using all SNPs, σ 2 g is additive genetic variance), and e is a vector of random residual effects (e~N(0,Iσ e 2 )).
Model effects and variance components for the GBLUP and Bayesian Power Lasso models were estimated by DMU and Bayz packages, respectively. For the GBLUP models, the narrow-sense genomic heritability corresponding to records of single plots was calculated as: where d(G 2 ) is the average diagonal element of the G-matrix (calculated using all SNPs), σ 2 g is additive genetic variance, and σ 2 e is residual variance. The following cross-validations (CVs) were used to study the effectiveness of possible strategies for implementing genomic selection in breeding programs: Leave-one-out (LOO): The GEBV of each line was predicted from the rest of the lines. The training set used in the LOO strategy was as large as possible (634 lines), and the genetic relationship between lines in the training and validation set was higher compared to the other CV strategies.
Leave-family-out (LFO): The GEBVs of lines in each half-sib family were predicted from lines of the remaining families. The average size of the half-sib families was 46 lines. Using this strategy, the effect of the genetic relationship between the lines in training and validation sets was studied.
Leave-set-out (LSO): The GEBVs of lines in each set were predicted from lines from the other set. The training set sizes were 314 or 321 lines. The LSO CV strategy was used for studying the predictive ability when GEBVs of lines from one breeding cycle were predicted from lines from another breeding cycle.
k-fold: The lines were randomly divided into k folds (2, 5, or 10) of equal size. The GEBVs of lines in each fold were predicted from lines in the other folds. The training set sizes were approximately 318 lines for the 2-fold, 508 lines for the 5-fold, and 572 lines for the 10-fold. Approximately half of the lines in the training sets were from set2014 and half from set2015. The k-fold CV strategy was used for studying the effect of the training set size.
Furthermore, the effect of the training set size was studied by selecting from 10% to 90% of the 635 lines as a training set for genomic predictions using LOO CV. The lines were randomly selected, and selection and predictions were repeated 100 times for each 10% interval.
Correlations between observed phenotypes corrected for fixed effects and GEBVs were calculated to determine predictive abilities of the models and were compared with the maximum correlation (the square root of the narrow-sense genomic heritability). Biases of the genomic predictions were calculated as the deviation from the expectation of the slope (1.0) of the regression line of the corrected phenotypes on the GEBVs.
Genomic feature models with two G-matrices were tested in order to possibly utilize the most significant SNPs better [53,54]. The lines were randomly divided in two folds. The lines of one fold were used for genome-wide association studies (GWAS), and the remaining lines were used for LOO genomic predictions. The SNPs used for computing the two G-matrices were selected for each trait based on their p-value in the GWAS. The most significant SNPs were used for one G-matrix, and the remaining SNPs were used for the other. The following thresholds for number of SNPs to include in the group of most significant were tested: 5, 10, 50, 100, 500, 1000, 3000, 5000, 7000, 10,000, and all 10,802 SNPs. The following model was used for the genomic predictions: where yis a vector of observed phenotypes, X, Z 5 , and Z 6 are design matrices, bis a vector of fixed effect (mean and year/set), sand nare vectors of additive genetic effects (s~N(0,G s σ 2 g s ) and n~N(0,G n σ 2 g n ), where G s and G n are G-matrices computed as above (2) using significant (G s ) or nonsignificant (G n ) SNPs, σ 2 g s and σ 2 g n are additive genetic variances, and eis a vector of random residual effects (e~N(0,Iσ e 2 )).

Phenotyping and Genotyping
A total of 635 F 6 winter wheat lines from two different breeding cycles (set2014 and set2015) were phenotyped for the quality traits flour yield and Alveos P, L, and W (Table S1). The phenotypic distribution for each trait is shown in Figure 1. Phenotypic variation was higher for the Alveograph traits than for flour yield ( Table 1). The wheat lines were genotyped for 10,802 SNPs (Table S2). A G-matrix was computed using all SNPs. A dendrogram based on the G-matrix showed that the lines were genetically related both within and between sets ( Figure 2). In total, the two sets consisted of 159 full-sib families with an average of four full-sibs per family.   For all traits, additive genetic variance was observed. Variance components were estimated based on GBLUP models and were used for estimation of narrow-sense genomic heritabilities. Heritabilities ranged from 0.38 for flour yield to 0.72 for Alveo W ( Table 2).      For all traits, additive genetic variance was observed. Variance components were estimated based on GBLUP models and were used for estimation of narrow-sense genomic heritabilities. Heritabilities ranged from 0.38 for flour yield to 0.72 for Alveo W ( Table 2).  For all traits, additive genetic variance was observed. Variance components were estimated based on GBLUP models and were used for estimation of narrow-sense genomic heritabilities. Heritabilities ranged from 0.38 for flour yield to 0.72 for Alveo W ( Table 2).

GWAS
Single marker regression was performed for each of the 10,802 SNPs ( Figure 3). Two linked SNPs on chromosome 5DS were significantly associated with flour yield, Alveo P, and Alveo W. On chromosome 1DL, significantly associated SNPs were identified for Alveos P and W. Additionally, SNPs associated with Alveo W were identified both on the short arm and on the long arm of chromosome 1B. A region on chromosome 4AL was significantly associated with Alveo L. The frequencies of the SNP alleles that were positively associated with each trait ranged from 13% to 64% (Table 3). A large difference was observed in Alveo P and in Alveo W for lines with the positive alleles of all significant SNPs compared to lines with the negative allele of one or more of the SNPs (Tables 4 and 5).

GWAS
Single marker regression was performed for each of the 10,802 SNPs ( Figure 3). Two linked SNPs on chromosome 5DS were significantly associated with flour yield, Alveo P, and Alveo W. On chromosome 1DL, significantly associated SNPs were identified for Alveos P and W. Additionally, SNPs associated with Alveo W were identified both on the short arm and on the long arm of chromosome 1B. A region on chromosome 4AL was significantly associated with Alveo L. The frequencies of the SNP alleles that were positively associated with each trait ranged from 13% to 64% (Table 3). A large difference was observed in Alveo P and in Alveo W for lines with the positive alleles of all significant SNPs compared to lines with the negative allele of one or more of the SNPs (Tables  4 and 5).    GWAS were also performed using Bayesian Power Lasso models to fit all 10,802 SNPs simultaneously [50]. The optimal value for β was 0.4 for flour yield and for Alveo W, 0.6 for Alveo P, and 0.8 for Alveo L, respectively. The SNPs most significantly associated with flour yield, Alveo P, and Alveo W according to the single marker regressions were also the SNPs with the highest genetic effects according to the Bayesian Power Lasso models ( Figure 4). For Alveo L, each SNP had very low genetic effect (less than 0.2).

Genomic Predictions
Genomic predictions based on all 10,802 SNPs using a GBLUP model were evaluated using different CV strategies. Predictive abilities of the models were determined as the correlations between observed phenotypes corrected for fixed effects and the GEBVs. The predictive abilities were intermediate to high, ranging from 0.50 for flour yield to 0.79 for Alveo W based on the LOO CV ( Figure 5). The LFO and the LSO CVs resulted in lower predictive abilities. The lowest predictive ability was 0.3 for flour yield based on the LSO CV. The predictive abilities of the k-fold CVs increased slightly when using a higher number of folds, and they were very close to the LOO when using 10 folds. The predictions were unbiased for the LOO and the k-fold CV and only slightly biased for the LFO and the LSO (Table 6).

Genomic Predictions
Genomic predictions based on all 10,802 SNPs using a GBLUP model were evaluated using different CV strategies. Predictive abilities of the models were determined as the correlations between observed phenotypes corrected for fixed effects and the GEBVs. The predictive abilities were intermediate to high, ranging from 0.50 for flour yield to 0.79 for Alveo W based on the LOO CV ( Figure 5). The LFO and the LSO CVs resulted in lower predictive abilities. The lowest predictive ability was 0.3 for flour yield based on the LSO CV. The predictive abilities of the k-fold CVs increased slightly when using a higher number of folds, and they were very close to the LOO when using 10 folds. The predictions were unbiased for the LOO and the k-fold CV and only slightly biased for the LFO and the LSO (Table 6).   For every CV strategy, the predictive abilities of the Bayesian Power Lasso model were a little better compared to the GBLUP for flour yield, Alveo P, and Alveo W, but not for Alveo L ( Figure 6).    For every CV strategy, the predictive abilities of the Bayesian Power Lasso model were a little better compared to the GBLUP for flour yield, Alveo P, and Alveo W, but not for Alveo L ( Figure 6).    For every CV strategy, the predictive abilities of the Bayesian Power Lasso model were a little better compared to the GBLUP for flour yield, Alveo P, and Alveo W, but not for Alveo L ( Figure 6).  The effects of training population sizes ranging from 10% to 90% of the 635 lines were studied (Figure 7). Predictive abilities increased from 0.33 at 10% (64 lines) to 0.50 at 90% (572 lines) for flour yield. For all traits, the predictive abilities increased when increasing the size of the training set, but the increases were smaller at larger sizes. Additionally, the variation around the mean of the predictive abilities decreased when increasing the size of the training set. The effects of training population sizes ranging from 10% to 90% of the 635 lines were studied (Figure 7). Predictive abilities increased from 0.33 at 10% (64 lines) to 0.50 at 90% (572 lines) for flour yield. For all traits, the predictive abilities increased when increasing the size of the training set, but the increases were smaller at larger sizes. Additionally, the variation around the mean of the predictive abilities decreased when increasing the size of the training set. Genomic features models with two G-matrices were also tested (Figure 8). One G-matrix was calculated from significant SNPs, and the other G-matrix was calculated from nonsignificant SNPs. The number of SNPs to include as significant ranged from five to all 10,802 SNPs. For flour yield, Alveo P, and Alveo W, predictive abilities were highest (0.52, 0.73, and 0.79, respectively) when fewer than 1000 SNPs were considered significant (10, 500, and 100 SNPs, respectively). For Alveo L, predictive abilities were highest (0.61) when 3000 SNPs were considered significant. For all traits, predictive abilities were higher when using the optimal number of significant SNPs for one G-matrix and the remaining SNPs for another G-matrix compared to using all SNPs for one G-matrix: 0.  Genomic features models with two G-matrices were also tested (Figure 8). One G-matrix was calculated from significant SNPs, and the other G-matrix was calculated from nonsignificant SNPs. The number of SNPs to include as significant ranged from five to all 10,802 SNPs. For flour yield, Alveo P, and Alveo W, predictive abilities were highest (0.52, 0.73, and 0.79, respectively) when fewer than 1000 SNPs were considered significant (10, 500, and 100 SNPs, respectively). For Alveo L, predictive abilities were highest (0.61) when 3000 SNPs were considered significant. For all traits, predictive abilities were higher when using the optimal number of significant SNPs for one G-matrix and the remaining SNPs for another G-matrix compared to using all SNPs for one G-matrix: 0.

Discussion
Wheat quality traits typically have intermediate or high heritabilities, although the traits can be considerably affected by environmental effects and genotype-environment (G×E) interactions [36][37][38]. Thus, additive genetic variation across environments also affects the traits. Here, narrow-sense genomic heritabilities ranged from 0.38 for flour yield to 0.72 for Alveo W (Table 2). Therefore, breeding for improved wheat quality traits should be possible.
Advanced breeding material was used in the present study. The breeding program has, until now, focused more on increasing yield rather than improving baking quality due to restrictions in application of nitrogen fertilization to the fields in Denmark, which have made it challenging to grow high quality bread wheat lines. Nevertheless, variation was observed for the quality traits, indicating that both high and low quality wheat lines were used as crossing parents for the studied lines, and that genetic variation was maintained throughout the breeding program ( Figure 1, Table 2). Six of the 96 crossing parents were used in both breeding sets. The dendrogram of the wheat lines indicated genetic relationships within each of the two sets ( Figure 2). However, the lines from each set were not clearly divided in the dendrogram, indicating that the lines were also related between the two breeding sets.
QTL for wheat quality traits were identified across the genome [1,2,17]. However, many QTL were only identified in certain environments or populations [55,56]. In the present study, two closely linked SNPs on chromosome 5D were significantly associated with flour yield, Alveo P, and Alveo W (Figure 3). The puroindoline genes Pina-D1 and Pinb-D1 have a large effect on grain hardness. These genes are located on chromosome 5D and can affect several quality traits [1,12,57]. Several Pinb-D1 alleles with a positive effect on wheat quality were identified [58,59]. However, relatively few of the markers in the present study were located on chromosome 5D [37], thus additional markers or sequencing would be needed to distinguish between each of the alleles.
For Alveo P and Alveo W, significant SNPs were also identified on chromosome 1D. The HMWG loci Glu-D1 is located on chromosome 1D and has a large influence on wheat baking quality [4]. Similarly, the glutenin loci Glu-B1 and Glu-B3, which are located on chromosome 1B, can affect

Discussion
Wheat quality traits typically have intermediate or high heritabilities, although the traits can be considerably affected by environmental effects and genotype-environment (G×E) interactions [36][37][38]. Thus, additive genetic variation across environments also affects the traits. Here, narrow-sense genomic heritabilities ranged from 0.38 for flour yield to 0.72 for Alveo W (Table 2). Therefore, breeding for improved wheat quality traits should be possible.
Advanced breeding material was used in the present study. The breeding program has, until now, focused more on increasing yield rather than improving baking quality due to restrictions in application of nitrogen fertilization to the fields in Denmark, which have made it challenging to grow high quality bread wheat lines. Nevertheless, variation was observed for the quality traits, indicating that both high and low quality wheat lines were used as crossing parents for the studied lines, and that genetic variation was maintained throughout the breeding program ( Figure 1, Table 2). Six of the 96 crossing parents were used in both breeding sets. The dendrogram of the wheat lines indicated genetic relationships within each of the two sets ( Figure 2). However, the lines from each set were not clearly divided in the dendrogram, indicating that the lines were also related between the two breeding sets.
QTL for wheat quality traits were identified across the genome [1,2,17]. However, many QTL were only identified in certain environments or populations [55,56]. In the present study, two closely linked SNPs on chromosome 5D were significantly associated with flour yield, Alveo P, and Alveo W (Figure 3). The puroindoline genes Pina-D1 and Pinb-D1 have a large effect on grain hardness. These genes are located on chromosome 5D and can affect several quality traits [1,12,57]. Several Pinb-D1 alleles with a positive effect on wheat quality were identified [58,59]. However, relatively few of the markers in the present study were located on chromosome 5D [37], thus additional markers or sequencing would be needed to distinguish between each of the alleles.
For Alveo P and Alveo W, significant SNPs were also identified on chromosome 1D. The HMWG loci Glu-D1 is located on chromosome 1D and has a large influence on wheat baking quality [4].
Similarly, the glutenin loci Glu-B1 and Glu-B3, which are located on chromosome 1B, can affect quality. No SNPs on chromosome 1A were significantly associated with the Alveograph traits. The loci Glu-A1 and Glu-A3 are also known to affect wheat quality traits [4]. Thus, these loci might not be polymorphic in the studied material, their effects might be too low to detect, or the genotyped SNPs might not be located in or close enough to the loci. Several markers may be required in order to distinguish between the alleles of each of the Glu loci [60]. Therefore, the SNP chip array used to genotype the studied wheat lines can perhaps not capture all the genetic variants. Characterization of Gluand Pin loci in the breeding material could be useful for more accurate estimation of their effects and consequently more accurate predictions of the baking quality traits [4,38,58].
The GWAS indicated that only few QTL with large or intermediate effects controlled the quality traits. The identified QTL only explained a relatively small proportion of the total genetic variance. This indicates that the traits were also controlled by many QTL with small effects. Identification of such minor QTL is challenging, especially if they have a low MAF or if they are located near major QTL. Lines with the positive alleles of the four largest QTL for Alveo W had considerably higher dough strength than lines with negative alleles of any of the QTL (Table 5). However, the four QTL together explained only 26.3% of the additive genetic variance based on the effects estimated from the single marker regressions, and these effects were possibly overestimated due to the Beavis effect [61,62]. The effects were lower when estimated using the Bayesian Power Lasso (Figure 4). Here, all SNP effects were estimated simultaneously, thus the effects were shrunk towards zero, and each QTL effect might have been distributed across several SNPs.
Since only a relatively small proportion of the genetic variance could be explained by the identified QTL, genomic predictions based on a large number of genome-wide markers could be useful. The different CV strategies showed that the predictive abilities were affected by size of the training set, by the genetic relationship between lines in training and validation sets, and by the G×E interactions ( Figure 5). The LSO strategy represents one way of implementing genomic predictions in breeding programs. Predicting GEBVs of new lines based on lines from previous years could possibly enable selection before any phenotypic information is available for the new lines. However, the LSO strategy resulted in lower and more biased predictive abilities than the other CV strategies ( Figure 5, Table 6). Possible reasons for the lower predictive abilities could be the size of the training set, the genetic relationship between lines, and the G×E interactions. Including lines from the same year and the same families in the training set improved the predictive abilities considerably. Reducing the size of training sets had a negative impact on the predictive abilities ( Figure 7). However, the decrease for very small training sets seen in the present study was not as drastic as in other studies [29,32], possibly due to the high heritabilities of the traits included in this study. Thus, a few hundred lines might be enough for the training set, if they are highly related to the validation set. The G×E interactions and the genetic relationship between lines might be partly confounded, since the lines were only tested from one location. These effects had a larger impact on predictive abilities than the size of the training set and the model used for the predictions. G×E interactions can be accounted for in, for example, reaction norm models if phenotypic data are available for lines replicated across several locations or years. Additionally, data about climatic conditions and soil types might be included to obtain higher predictive abilities [63]. Such models could also be used for selecting lines for target environments or lines that are performing well over many environments. Predictive abilities can also be affected by the number of markers used for the predictions [29,32]. If the markers are selected based on GWAS, the number can possibly be reduced to a few hundred markers without affecting the predictive abilities. However, the markers that are selected would not be the same for each trait, and the markers could change after each new breeding cycle [29].
The predictive abilities of the Bayesian Power Lasso models were slightly higher compared to the GBLUP models for flour yield, Alveo P, and Alveo W ( Figure 6). In the Bayesian Power Lasso model, the difference between small and large QTL effects can be bigger than in GBLUP [50]. Thus, the largest QTL effects might be shrunk too much in the GBLUP models. For Alveo L, no improvements were observed when using the Bayesian models, indicating that no major QTL were present for this trait ( Figure 6). Previous studies have also shown that different types of Bayesian models can, in some cases, give slightly more accurate predictions compared to GBLUP models, especially for traits influenced by major QTL and for populations with low genetic relationships between training and validation sets [26,37,50].
Other studies of genomic selection for wheat quality traits using breeding material have reported predictive abilities in similar ranges as in the present study [36,37,64]. Thus, genomic prediction appears to be a promising strategy for improving wheat quality in breeding programs. Nevertheless, the use of information from GWAS or from known major QTL in genomic predictions might be useful for modeling marker effects more accurately. Markers for major QTL can be specified as having fixed effects, while remaining markers have random effects. Bernardo (2014) [39] recommended to specify fixed effects for markers that explain more than 10% genetic variance for traits controlled by fewer than 10 major genes. Zhao et al. (2014) [26] used W-BLUP (weighted BLUP) to give few functional markers a larger weight than other markers. This improved prediction accuracies of heading time and plant height in hybrid wheat compared to marker-assisted selection (MAS), ridge regression BLUP (equivalent to GBLUP) and BayesCπ. Improved prediction accuracies were also reported by Arruda et al. (2016) [34] when using fixed effects for QTL associated with Fusarium head blight resistance traits compared to ridge regression BLUP, and similarly for pre-harvest sprouting tolerance by Moore et al. (2017) [40]. However, including QTL identified using the same lines that were also used for the genomic predictions could lead to inflation of the accuracies [34]. In a recent study by Michel et al. (2018) [38], prediction accuracies for several wheat quality traits were higher based on ridge regression BLUP compared to specifying fixed effects for associated markers identified in independent populations. However, including one or more of the three Glu-1 loci as fixed effects resulted in improved accuracies. For each trait, Glu-1 loci were included if the locus explained more than 5% genetic variance [38].
In the present study, an alternative approach using one G-matrix computed from the most significant SNPs and another G-matrix computed from the remaining SNPs could slightly improve predictive abilities (Figure 8). However, the optimal significance threshold depended on the trait. For the traits controlled by few QTL with large effects, a conservative significant threshold seemed to be best, while a loose threshold seemed to be better for traits controlled only by QTL with small effects. To avoid inflation of the predictive abilities, half of the lines were used for the GWAS to select significant SNPs, and the other half were used for the genomic predictions. Thus, larger datasets might be necessary to study the predictive abilities of the genomic feature models and the optimal number of SNPs more thoroughly.
Implementation of genetic markers and genomic predictions in breeding programs could likely lead to increased genetic gains for the wheat quality traits. Resources needed for phenotyping could be reduced, and the selection intensity could be increased. The identified SNPs with large effects might be used for screening large numbers of lines early in the breeding program and selecting lines that potentially have good quality traits. Genomic predictions could be used for a more accurate selection of lines with good quality in later generations.

Conclusions
SNPs significantly associated with flour yield, Alveo P, and Alveo W were identified on chromosome 5D. For Alveos P and W, associated SNPs were also identified on chromosome 1D. Likely candidate genes could be the Pina-D1 or the Pinb-D1 on chromosome 5D and the Glu-D1 loci on chromosome 1D. Furthermore, SNPs associated with Alveo W were identified on chromosome 1B, and SNPs associated with Alveo L were identified on chromosome 4A. Additive genetic variance explained by a single SNP was up to 13.3% (SNP on 5D for flour yield). The identified SNPs can be used in early generations of breeding programs to screen large numbers of lines. In later generations, it would be advantageous to use a large number of SNPs to ensure accurate prediction of breeding values. Predictive abilities of GBLUP models were 0.50 for flour yield, 0.75 for Alveo P, 0.79 for Alveo W, and 0.64 for Alveo L based on the LOO CV. Predictive abilities were lower when using smaller training sets but were still moderate when using only 10% of the 635 lines as the training set. Furthermore, predictive abilities were significantly lower when using LSO and LFO CV strategies because of reduced genetic relationship between lines in training and validation sets and because of G×E interactions. Predictive abilities were similar or slightly higher based on Bayesian Power Lasso and genomic feature models. Thus, GBLUP models could be used for genomic prediction of wheat quality traits with moderate to high predictive ability. Other models might give slightly higher predictive abilities for traits where major QTL are present in the breeding material.