Genomic Selection and Genome-Wide Association Studies for Grain Protein Content Stability in a Nested Association Mapping Population of Wheat

Grain protein content (GPC) is controlled by complex genetic systems and their interactions and is an important quality determinant for hard spring wheat as it has a positive effect on bread and pasta quality. GPC is variable among genotypes and strongly influenced by the environment. Thus, understanding the genetic control of wheat GPC and identifying genotypes with improved stability is an important breeding goal. The objectives of this research were to identify genetic backgrounds with less variation for GPC across environments and identify quantitative trait loci (QTLs) controlling the stability of GPC. A spring wheat nested association mapping (NAM) population of 650 recombinant inbred lines (RIL) derived from 26 diverse founder parents crossed to one common parent, ‘Berkut’, was phenotyped over three years of field trials (2014–2016). Genomic selection models were developed and compared based on predictions of GPC and GPC stability. After observing variable genetic control of GPC within the NAM population, seven RIL families displaying reduced marker-by-environment interaction were selected based on a stability index derived from a Finlay–Wilkinson regression. A genome-wide association study identified eighteen significant QTLs for GPC stability with a Bonferroni-adjusted p-value < 0.05 using four different models and out of these eighteen QTLs eight were identified by two or more GWAS models simultaneously. This study also demonstrated that genome-wide prediction of GPC with ridge regression best linear unbiased estimates reached up to r = 0.69. Genomic selection can be used to apply selection pressure for GPC and improve genetic gain for GPC.


Introduction
Grain protein content (GPC) is a high-priority determinant of end-use quality for most cereals [1], including pasta (Triticum turgidum L.) and bread wheat (Triticum aestivum L.), for which higher GPC is preferred. Due to the large dependence on wheat, rice (Oryza sativa L.), and maize (Zea mays L.) as primary sources of carbohydrates, they contribute 80% of dietary calorie requirements for humans [2]. Compared to other agricultural commodities, cereal grains contain a relatively low concentration of protein. In a screening of 12,600 lines from the USDA world wheat collection, GPC varied from 7% to 22%, with the genetic component accounting for only a third of the variation [3]. Breeding efforts to improve GPC have been difficult due to strong environmental influences and the high variability of GPC across years and locations [4], combined with variable economic value-based trade-offs between starch and protein yields in grains [5].
Phenotypic plasticity refers to the flexibility that allows for changes in a particular trait as a result of environmental variability [6]. Improved understanding of genotypes and environmental interactions could provide new strategies for breeding crop varieties that stably perform between changing environments [7]. One of the key challenges to identifying genes controlling environmental stability lies in the quantification of stability as a trait [8]. There are various approaches for measuring environmental stability, including the use of variance components of individuals across environments (represented by a coefficient of variation), the comparison of mean responses of genotypes to the overall mean of individuals in the trial (calculated as a coefficient of regression on an environment index) [9,10], and expected change in the performance of a genotype as a function of environmental effect [11]. Compared to unstructured genotype by environment interaction models, the Finlay-Wilkinson (FW) regression fits every level of genotype and environment and reveals genotype performance across environments [12]. GPC stability can potentially help select stable genotypes across increasingly unpredictable environments [13]. Various biparental and association mapping studies have identified QTLs controlling GPC, but none of them has focused on the stability of GPC [14,15].
Nested association mapping (NAM) is a multi-parental design with higher allelic variation than biparental populations and stronger statistical power than association mapping populations, combining the advantages of each approach [16,17]. A universal parent is crossed to multiple genotypes, followed by inbreeding to make a combination of both full-sib and half-sib recombinant inbred lines (RILs) [18]. NAM has proven its success in mapping complex traits in barley (Hordeum vulgare L.) [19], maize [20], rice [21], and arabidopsis (Arabidopsis thaliana) [22]. The NAM population design mirrors the structure of many breeding programs, where multiple superior, diverse, or exotic lines are crossed with a few elite breeding lines for population development in elite and pre-breeding pipelines. The cross of exotic germplasm with elite breeding lines for population development aids in normalizing the genetic background and selecting segregating alleles in favor of the adapted parental alleles [23]. The resolution and power of NAM populations allow for the assessment of complex traits like GPC stability in structured germplasm via genome-wide association studies (GWAS) [24].
Joint linkage association mapping is applied in multi-parental mapping populations where QTL terms are nested within families [25] instead of testing marker effects across the families, as in GWAS [26]. Larger allelic classes and more balanced allele frequencies increase the power to detect QTLs with small effects. NAM combines better resolution by targeting historical recombination events in parents and also includes more causative events that are likely to segregate in biparental progenies [27]. Thus, NAM is an effective tool for identifying major and minor QTLs associated with a particular trait [28].
In plant breeding, the accumulation of minor QTLs is a major constraint given the dozens of loci and resource-limited population sizes. Hence, genomic selection (GS) is used to sum the effects of genome-wide markers to predict genomic estimated breeding values (GEBVs) [29]. In GS, a training population is genotyped and phenotyped for the traits of interest to estimate the genetic effect of each marker. The success of GS depends upon the predictive accuracy of the GS models, which is measured as the correlation between estimated breeding values and the observed phenotypic values of the selected populations [30]. GS can translate into higher genetic gains by reducing the number of progeny and cycles needed and also by improving selection intensity and reducing cycle time [31,32].
Several studies have reported the application of GS for wheat yield and disease resistance [33,34] but not thoroughly for quality traits, especially GPC and GPC stability. GS has been tested in soft winter wheat for end-use quality. Heffner et al. found that end-use quality and processing traits are more predictive than grain yield [35]. The prediction accuracy for flour yield and flour protein was 0.56 and 0.39, respectively, using a ridge regression model. GS accuracy for different agronomic traits and their stability was predicted in a winter wheat population of 273 lines, and it was observed that GS accuracy varied from 0.33 to 0.67 for yield, with yield stability having a higher accuracy [36].
Although several studies have focused on genomic selection for yield stability, none of them has investigated GPC stability. The main objectives of this study were to (1) detect marker-trait associations for GPC stability and GPC and (2) identify the ability of GS models to predict GPC and GPC stability under cross-and independent validations.

Plant Material and Trait Measurement
Thirty-two spring wheat accessions from the USDA-ARS National Small Grains Collection were chosen as parental lines for the creation of the NAM population. These parents were crossed with the common cultivar 'Berkut' ('Irena'/'Babax'/'Pastor'; released in 2002) to create 32 half-sib families [37]. Berkut was included in this study because it is a broadly adapted photoperiod insensitive and semi-dwarf cultivar developed by the International Maize and Wheat Improvement Center (CIMMYT), Mexico. More information about the population development and field experiment is described in Sandhu et al. [37,38].
Twenty-six families with genotype data provided by Kansas State University [39,40] were then selected, and 25 random RILs from each family (650 total RIL; named NAM 650 ) were planted between 2014 and 2016 at the Spillman Agronomy Farm in Pullman, WA, under rainfed conditions. A modified augmented design was used in each trial with three check cultivars (Berkut, 'Thatcher', and 'McNeal' [41]). Planting was completed on 5 May 2014, 8 April 2015, and 10 April 2016 based on field conditions. The agronomic practices, including nitrogen fertilizer goals, were uniform in each environment, assuming a 4.3 t/ha average grain yield goal. The percentage of protein content in the grain was measured using a Perten DA 7000 NIR analyzer (Perkin Elmer, Springfield, Sweden). High-performance liquid chromatography (HPLC) is an instrument used to measure GPC, but in our study NIR was used as a measure. This is a commonly used tool to measure GPC and approved by the American Association of Cereal Chemists (Method 39-25.01). This method is also approved by the United States Department of Agriculture to meet all domestic and export requirements and specifications. Our NIR machine was calibrated with over 2000 samples generated from the USDA Western Wheat Quality Laboratory, and the calibration model was shown to have a 99.3% accuracy in predicting GPC. Since the late 1990s, NIR has become a standard method for measuring GPC [42,43]. Furthermore, NIR has been used in multiple studies to determine GPC. Most notable are those which used NIR to fine map and identify the impact of the GPC-B1 locus in wheat [44,45]. Grain yield was obtained using grain weight per plot with a Wintersteiger Nursery Master combine (Ried im Innkreis, Austria).

Statistical Analysis
Best linear unbiased estimates (BLUEs) were calculated with the augmented complete block design (ACBD) in the R program for individual environments [46,47] using the model: where Y ij is the grain protein content of an individual line, µ is the mean effect, Block i is the fixed effect of the ith block, Gen j represents the fixed effect of unreplicated genotypes, Check j is the effect of the replicated checks within each block, and e ij is the standard normal error. All the effects were considered fixed in BLUE calculations.
The significance of differences in GPC were analyzed across years and between families using analysis of variance (ANOVA). Standard deviation and coefficients of variation were calculated within the different families to identify those with lower variation for GPC across years. The GPC variation within families in each environment and across years was used to select families that had less variation for GPC. Broad sense heritability for GPC was obtained as: where σ 2 g and σ 2 e are the genotypic and error variance components, respectively.

Stability Analysis
The stability of each RIL was assessed using a Finlay-Wilkinson (FW) program implemented in the FW package in R [12]. The FW package jointly estimates the parameters of the FWR equation: where Y ij is the GPC of the ith RIL from the j th environment, µ is the mean effect, t j is the main effect of the j th environment, g i is the main effect of the ith RIL, (1 + b i ) is the change in expected performance of the ith RIL per unit change in the environmental effect (Stability index), and e ij is an error term assumed to be independent and identically distributed with mean zero and variance σ e 2 . All parameters are treated as random effects with distributions: g~N(0, Aσ 2 g ), b~N(0, Aσ 2 b ), and t~N(0, Tσ 2 t ), where A is an N × N kinship matrix for all the RILs calculated from the complete set of genetic markers, and T is an M X M Pearson variance-covariance matrix, describing the relationship of phenotypic values among M environments.
The stability of each RIL was calculated using GPC values for a selected 175 RILs (NAM 175 ) in each environment as dependent variables. These RIL families were selected due to less variation in GPC across the environments, and hence were used for further analysis instead of the full 650 RILs. The selected RILs represent seven RIL families that demonstrated higher GPC stability across environments and showed segregation for stability within the family. The FW package returned the environmental stable genotypic effect (g i ), environmental effect (t j ), and the stability index of each RIL (b i ). The stability index for each RIL provided an idea of plasticity across the environments [7]. A stability index of 1 and −1 means that a genotype is highly plastic and responds according to environmental changes whereas a stability index of 0 suggests that the genotype performs stably under different environments. Loci associated with GPC stability were identified by GWAS of the stability index absolute values. The parents for the selected 175 RILs were 'Dharwar Dry', 'PI210945', 'CItr15144', 'PI92569', 'PI92569', 'CItr4174', and 'PI43355' (Tables S1 and S2).
Furthermore, grain protein content deviation (GPD) was obtained using grain yield and GPC information across the environments for the NAM 175 population [48]. GPD represents the relationship between GPC and grain yield, which helps a plant breeder to make a selection. A linear regression model was fitted on grain yield and GPC to derive the GPD for the selected population using residuals form the model [49]. The GPD distribution in this study varied from −1.7 to 1.32, where values above zero mean higher deviation of GPC with higher grain yield. Furthermore, GPD information was used to identify the loci controlling this trait in the NAM175 population.

Genotyping
The NAM population genotyping, curation methods, and population maps were previously reported [33,37,38,40]. The initial genotypic data used in this study comprised 73,345 high-quality markers which were anchored to the Chinese Spring RefSeqv1 reference map [50] and NAM 175 selected for stable GPC were used. Individual RILs with missing GPC data were removed before filtering. Markers with more than 20% of data missing were removed for further analysis. Individual RILs missing more than 10% allele data were removed before culling based on markers with less than 5% minor allele frequency (MAF) to lower the probability of Type I errors during GWAS. After filtering, 175 RILs from NAM 175 and 38,588 markers were used for GWAS analysis.

Population Structure and Genome-Wide Association Studies
Kinship matrix and structure parameters were calculated by GAPIT [51]. The Van-Raden algorithm was used to derive the kinship parameter from marker genotypes [52]. Population structure was analyzed using principal component analysis (PCA) implemented with the R function 'prcomp' in the software GAPIT [53]. Genetic relatedness among NAM 175 and NAM 650 was calculated using a population differentiation coefficient (F st ) through the "Population Measures" function in JMP [54]. The GWAS was conducted in the GAPIT R package using a mixed linear model (MLM) (Q + K model) [55], a compressed mixed linear model (CMLM) [56], a fixed and random model circulating probability unification (FarmCPU) [57], and the Bayesian information and Linkage-disequilibrium Iteratively Nested Keyway (BLINK) model [58]. Multiple models were used for the GWAS analysis as this would help us to reduce the number of false positives. A general description of all the above models is provided below.
MLM includes the kinship matrix as a random effect in the mixed effect model, and MLM can be represented as: where Y is a matrix of phenotypic information, SNP represents the matrix of markers, Q represents the population structure, and Kinship represents the relationship matrix between the individuals included in the model. SNP and Q are set as fixed effects, while kinship is a random effect in the model [55].
MLM was computationally very intensive because computational time varies with the third power of the number of individuals in the random effect model. Furthermore, there were confounding issues between testing marker, structure, and kinship matrix, as the same set of markers were double counted. CMLM clustered the individuals into different groups, resulting in a reduction of the effective size of the random effect model [56]. CMLM obtains the kinship among the groups and is computationally more efficient than MLM. CMLM can be represented as: where Kinship is the relationship matrix among the groups and other terms are the same as in the MLM, described above.
Both of the above models were single locus models which are not the best for handling complex traits; hence, we used FarmCPU and BLINK for comparison. The working principle of FarmCPU is divided into fixed and random effect models. The fixed effect model tests single markers at a time with multiple associated markers as a covariate to control for false positives. Furthermore, model overfitting is avoided in the random effect model by obtaining kinship using multiple associated markers. The p-value of each tested and associated marker is unified at each iteration. The FarmCPU model can be represented as: This is the fixed effect component of the FarmCPU model, with individual markers tested one at a time and other terms of the equation as described previously.
This is the random effect component of the FarmCPU model, and all terms of the equation are as described previously [57].
FarmCPU has an efficient fixed model, but it has a computationally expensive random effect model. Furthermore, QTNs in random effect models were selected based on their even distribution over the genome. In order to increase computational efficiency, the random effect model was replaced with a fixed effect model using Bayesian information criteria. This new method is known as the Bayesian information and linkage disequilibrium iteratively nested keyway (BLINK), as QTNs are selected based on linkage disequilibrium information [58].
Association studies were performed on the NAM 175 using a GPC stability index, environmentally stable phenotypic GPC values were derived from the FW regression, and grain protein content deviation (GPD) was combined across all environments. Model selection function from the GAPIT manual and the Q-Q plot results were used to decide the number of PCs for inclusion in the GWAS models. PCA groups were included as covariates in the GWAS model to account for population structure. We also performed a pairwise correlation among the GWAS models results and its value varied between 0.99 (MLM and CMLM), 0.93 (MLM and FarmCPU), and 0.95 (MLM and BLINK). The Bonferroni correction with a stringent α = 0.05 was used to identify highly significant associations. Bonferroni corrected p-values (Bonferroni adjusted p-values) were calculated using alpha/number of tests performed and alpha was set to 0.05. This is one of the stringent tests and allowed us to reduce the number of significant markers [59].

Genomic Selection
Genome-wide marker effects for GPC and GPC stability were estimated using ridge regression best linear unbiased prediction (rrBLUP) [60], according to the model: where y is an N × 1 vector of BLUEs for GPC or GPC stability for each RIL, µ is the overall mean, Z is an N × M matrix of markers, u is a vector of marker effects, and e is a vector of residuals. GS was performed with five-fold cross-validation, including 80% of the samples in the training population and predicting the GEBVs of the remaining 20% of the samples under each environmental condition. For accuracy assessment, 250 replication sets were performed, each replication consisting of five model iterations.
Genomic selection models were developed for the NAM population of 650 RILs (NAM 650 ) and separately for the 175 RILs (NAM 175 ) from seven families selected based on lower variation for GPC across environments. Furthermore, independent validations were performed using both sets of RILs for predicting GEBVs for GPC. During independent validations, GS models were trained on the previous year's data set, and predictions were made for the upcoming year. Models trained on 2014 GPC data were used for predictions in 2015 and 2016. Similarly, the 2015 GPC training model was used for 2016 predictions.

Variation of Grain Protein Content across Environments
The GPC values of the NAM 650 population ranged from 11.  Tables S1 and S2. Using this data, seven families totaling 175 RILs (NAM 175 ) were selected for low variation in GPC to identify loci controlling stability of the GPC (Table S3).

Stability Analysis
Stability index (b), environmental effect (t), and genotypic effect (g) values were obtained for the NAM 175 population. Absolute values of the stability index for NAM 175 RILs ranged from 0.00 to 2.504 and were normally distributed. Thirty RILs were identified which had no significant difference of stability index from 0 using individual t-test (p < 0.05). Environment 2014 was observed to be the most favorable for high GPC. The NAM 175 population was divided into five categories based on the GPC and stability index, the trends of which are depicted by solid lines in Figure 2. Categories 2, 3, and 5 include the stable GPC lines, with GPC in the range of 15-17%, 14-15%, and 9-14%, respectively (

Stability Analysis
Stability index (b), environmental effect (t), and genotypic effect (g) values were obtained for the NAM175 population. Absolute values of the stability index for NAM175 RILs ranged from 0.00 to 2.504 and were normally distributed. Thirty RILs were identified which had no significant difference of stability index from 0 using individual t-test (p < 0.05). Environment 2014 was observed to be the most favorable for high GPC. The NAM175 population was divided into five categories based on the GPC and stability index, the trends of which are depicted by solid lines in Figure 2. Categories 2, 3, and 5 include the stable GPC lines, with GPC in the range of 15-17%, 14-15%, and 9-14%, respectively (Figure 2). Categories 1 and 4 represent the plastic lines, with reversal of GPC when moving to other environments. Categories 1, 2, 3, 4, and 5 include 81, 15, 9, 64, and 6 RILs, respectively.

Population Structure Analysis
Population structure analyzed with PCA separated the NAM175 population into seven separate groups based on their different parents, in addition to the common parent 'Berkut' (Figure 3). The PC1 accounted for 7.0% of the variation, whereas the PC2 explained 5.0% of the genetic variation ( Supplementary Figures S2 and S3). Kinship plots obtained from VanRaden algorithms in GAPIT also separated the population into seven different groups (data not shown). Inclusion of a different number of PCs as covariates in the GWAS model demonstrated that the first three PCs best control the false positives and false negatives, as evident from the Q-Q plots (Supplementary Figure S4), and, similarly, the model selection function based on BIC showed that first PCs should be used in the GWAS models. The Fst coefficient was 0.11 for the NAM175 and 0.18 for that NAM650, and this provided us with information about the genetic relatedness among the individuals within each population.

Population Structure Analysis
Population structure analyzed with PCA separated the NAM 175 population into seven separate groups based on their different parents, in addition to the common parent 'Berkut' (Figure 3). The PC1 accounted for 7.0% of the variation, whereas the PC2 explained 5.0% of the genetic variation ( Supplementary Figures S2 and S3). Kinship plots obtained from VanRaden algorithms in GAPIT also separated the population into seven different groups (data not shown). Inclusion of a different number of PCs as covariates in the GWAS model demonstrated that the first three PCs best control the false positives and false negatives, as evident from the Q-Q plots (Supplementary Figure S4), and, similarly, the model selection function based on BIC showed that first PCs should be used in the GWAS models. The F st coefficient was 0.11 for the NAM 175 and 0.18 for that NAM 650 , and this provided us with information about the genetic relatedness among the individuals within each population.

Marker-Trait Associations for the Stability of Grain Protein Content
Eighteen QTLs controlling GPC stability were identified with the help of four different GWAS models using a stringent Bonferroni correction of α = 0.05 (Table 1; Figure S4). The variation explained by each locus ranged from −4.19 to 3.98%. Out of these eighteen QTLs, eight were identified by two or more GWAS models simultaneously and will be referred to as high confidence QTLs (to be discussed later) (Table 1). Cumulatively, these loci explained 35.40% of the phenotypic variation. Out of those eight high confidence QTLs, three were located on chromosome 3B and chromosomes 1A, 2A, 4D, 5B, and 7D had one significant association. The QTLs on 1A, 2A, 3B, 4D, and 7D had a positive effect on increasing GPC stability, while three other OTLs had a negative effect on GPC stability. The parents of origin for associated alleles are provided in Table 1. Removal of those alleles from a breeding program by selecting for alternative alleles at these loci will favor the development of lines having increased GPC stability.

Marker-Trait Associations for the Stability of Grain Protein Content
Eighteen QTLs controlling GPC stability were identified with the help of four different GWAS models using a stringent Bonferroni correction of α = 0.05 (Table 1; Figure S4). The variation explained by each locus ranged from −4.19 to 3.98%. Out of these eighteen QTLs, eight were identified by two or more GWAS models simultaneously and will be referred to as high confidence QTLs (to be discussed later) (Table 1). Cumulatively, these loci explained 35.40% of the phenotypic variation. Out of those eight high confidence QTLs, three were located on chromosome 3B and chromosomes 1A, 2A, 4D, 5B, and 7D had one significant association. The QTLs on 1A, 2A, 3B, 4D, and 7D had a positive effect on increasing GPC stability, while three other OTLs had a negative effect on GPC stability. The parents of origin for associated alleles are provided in Table 1. Removal of those alleles from a breeding program by selecting for alternative alleles at these loci will favor the development of lines having increased GPC stability.
Eight significant QTLs were identified using the GPC values obtained across the environments as the response trait (Table 2; Supplementary Figures S5 and S6). Out of these eight QTLs, three were identified by two or more GWAS models simultaneously and will be referred to as high confidence QTLs (to be discussed later) ( Table 2). Cumulatively, these loci explained 24.80% of the phenotypic variation. Out of those three high confidence QTLs, each of them was located on chromosomes 2B, 7A, and 7B. The QTL on 2B and 7B positively affects GPC, while the QTL on 7A negatively affected GPC such that removing those alleles would increase GPC. As there are multiple alleles with favorable effects identified from various different parental lines, pre-breeding efforts will be required to introgress all the QTLs identified. As some QTLs were identified from landraces, selection during this pre-breeding process may help reduce linkage drag. Twelve significant QTLs were identified for grain protein content deviation obtained across the environments as the response trait ( Table 3). Out of these twelve QTLs, three were identified by two or more GWAS models simultaneously and will be referred to as high confidence QTLs (to be discussed later) (Table 3). Cumulatively, these loci explained 27.63% of the phenotypic variation. Loci controlling GPC stability, GPC, and GPD are different, demonstrating separate genetic architectures for each trait which could be selected simultaneously.

Prediction for Grain Protein Content and Stability
The prediction accuracy for GPC in the NAM 650 population was r = 0.50 in 2014, r = 0.55 in 2015, and r = 0.53 in 2016. Overall, GPC stability was less predictable with prediction accuracy values between r = 0.34 and 0.44 and a mean of 0.40. There was a significant difference in prediction accuracy for GPC and GPC stability, which was assessed using Tukey's test (p-value < 0.05; F statistics = 9.9). The prediction accuracy of GPC in the selected NAM 175 subset ranged between r = 0.56 to 0.69. The maximum prediction accuracy of r = 0.69 was achieved again for the 2015 environment, while the lowest prediction accuracy of r = 0.56 was obtained for the 2014 environment. Comparison of prediction accuracies for each environment using the two different sets of populations is presented in Figure 4. Prediction accuracy is generally high when using the NAM 175 , increasing by 15% as compared to the NAM 650 .     GS accuracy was significantly lower for the independent validation of each population set and under different environments compared to cross-validation GS accuracies assessed using Tukey's test (p-value < 0.05; F statistic = 9.2) ( Figure 5). Independent prediction accuracies for the NAM 650 population ranged between r = 0.30 to 0.39 and ranged between r = 0.35 to 0.43 for the NAM 175 for GPC. The independent prediction accuracies were higher for the NAM 175 , and the same results were observed during the cross-validation prediction scenario. The highest independent prediction accuracy was obtained by a training model on 2015 GPC to 2016 GPC ( Figure 5).

Stability of Genotypes across Environments
A primary goal of plant breeding programs is to select germplasms with superior adaptation to the targeted environments. The performance of genotypes may range from those that are very well adapted to a narrower set of environments and perform below average in others to genotypes that perform consistently relative to others across a wider range of environments and are considered to have greater stability. Herein, we used a static stability concept, targeting a predetermined value of GPC across all the environments. There are numerous statistical tools for analyzing static stability but here we applied FW regression analysis as it is capable of summarizing the interactions in comprehensible ways [9].
Crossover and non-crossover interactions were observed in our FW regression analysis which included reversals in rank and scale effects for GPC due to environmental effects and G*E interactions. Identification of only thirty lines with stable GPC across the environments out of 650 lines highlights the challenge of maintaining an adequate population size when selecting for a complex quantitative trait in genetically diverse germplasm. Evaluating this population in a greater number of environments would provide more insight into the genetic control of this trait. Additionally, the founder parents of the NAM population are primarily landraces that have not undergone routine selection for GPC [37], as has been performed with contemporary germplasms for at least the past 60 years. The stable lines varied in their average GPC, but lines having high and stable GPC (Category 2) demonstrate the ability to select genotypes for these traits. The lines identified in this study having high GPC and stability are particularly useful for introgression into modern breeding germplasm to expand available allelic variation for this important trait.

Genomic Regions Controlling Stability of GPC
The QTLs associated with GPC stability and GPC per se performance were not detected in the same genomic regions, suggesting that GPC stability and GPC are under the control of different genes. Accounting for the different genetic architectures of these two traits could aid the indexed selection for the desired GPC along with stability. Different QTLs controlling yield stability and yield per se in wheat have similarly been documented [61]. In maize, yield and yield stability were observed to be independent, demonstrating the potential for simultaneous selection for both traits [62]. Loci controlling yield stability were located in the same regions as QTLs for yield and yield-component traits in a barley mapping effort [63]. Critical factors that are not captured or quantifiable in each of these yield-related studies are the precise abiotic or biotic factors limiting yield potential in different environments. Thus, stability is an important aspect to consider, and our study suggests that breeding programs may select for both GPC and GPC stability by treating them as separate traits for developing varieties grown in climates with environmental variables that are difficult to predict.
Genomic regions controlling GPC stability have not been investigated in wheat based on the available literature. In the present study, eight high confidence QTLs controlling GPC stability distributed over six chromosomes were identified that individually explained −4.19 to 3.98% of variation for GPC stability. The QTLs on 1A, 2A, 3B, 4D, and 7D had a positive effect on GPC stability, while three other OTLs had a negative effect on GPC stability. Similar results were obtained by Sehgal et al. when mapping genomic regions for yield stability, where 11 QTLs associated with yield stability were distributed on seven chromosomes [64]. The amount of phenotypic variation explained by their yield stability QTLs varied from 3.2 to 8.1%. Thus, although stability is an important consideration when selecting for complex quantitative traits, such as yield and GPC, appropriate index weighting in genomic selection approaches is most likely an improved alternative in selection schemes.
Marker-trait associations have reported QTLs for GPC on all 21 chromosomes of wheat [65][66][67][68][69]. Loci controlling GPC were mapped to chromosomes 1A, 2B, 3B, 4A, 6B, 7A, and 7B in this study. The region linked to SpringWheatNAM_tag_190170 on 1A was previously identified in a DH population using composite interval mapping [70] and Groos et al. [66] used an F 7 RIL population grown in five environments. The favorable allele for this locus was identified in Berkut and Dharwar Dry, both of which are modern cultivars. The 7A GPC locus was discovered in the same region as another investigation [70,71]. These studies also reported that this locus had a negative effect on total GPC. Loci on chromosomes 3B and 4A have not been previously reported for GPC in wheat [14,67,72], which may be due to no prior utilization of these landraces in mapping studies [37]. These two loci had a negative effect on GPC, with associated alleles identified in the landraces. The utilization of diverse landraces in this study provides information about different genomic regions that are absent in present-day cultivars [73].
Similar to GPC stability loci, QTLs controlling GPC per se explained a small portion of the variance (24.80%). Given the number of small-effect loci controlling each trait, the GS approach proposed in this study should improve the ability to select for these traits simultaneously during the breeding process. The utilization of indexed genomic selection would assist in selecting simultaneously for GPC stability and GPC [74]. Rapp et al. [71] demonstrated the use of phenotypic and genomic selection indices to select durum wheat lines for high GPC and grain yield. Similarly, lines having high GPC and grain yield were selected using index selection in a multi-variate GS model for wheat [75]. These studies suggested that utilization of index selection in multi-variate GS models would potentially aid in selecting lines having stable GPC in addition to high GPC and grain yield.

Accuracy for Predicting GPC and GPC Stability
Prediction accuracy for GPC ranged between 0.50 and 0.69, which is moderately higher than prediction accuracies for grain yield [38,39]. Heritability of GPC, the effect of the NAM population, and population sizes used for training the GS models would each affect prediction accuracy [76,77]. Our results are consistent with previous studies where high heritability has resulted in better prediction accuracies in cereals [78]. The heritability of GPC is usually higher than yield, which ultimately resulted in better predictions for GPC [79]. In a genomic prediction experiment for grain yield in oats (Avena sativa L.), lower prediction accuracy was obtained because the experiment was planted under diverse environmental conditions, resulting in reduced genetic variance as compared to G*E interactions and thereby reduced heritability [80]. The NAM population in the current study was investigated in more homogenous target environments for wheat production in the PNW. This should result in a relatively lower G*E variance relative to genetic variance, leading to a higher heritability estimate and an increase in genomic selection accuracy.
Prior studies are not available to compare the accuracy of GS for GPC stability with stability index values obtained from FW regression in wheat. GS accuracies for GPC stability were significantly (p < 0.05) lower than GPC, suggesting a more complex architecture for GPC stability. Huang et al. [36] conducted GS for grain yield, test weight, and flour protein content stability using an additive main effect and multiplicative interaction (AMMI) model in wheat. They observed an accuracy from 0.14 to 0.31 for the stability of flour protein using four different GS models, namely, Bayesian ridge regression (0.14), elastic net (0.27), rrBLUP (0.31), and RKHS (0.31). Their study also demonstrated that the stability index for quality traits has less prediction accuracy than the trait itself. Our results, coupled with those of Huang et al. [36], suggest that GS can be used for predicting GPC stability. Furthermore, the prediction accuracy for GPC stability could be improved by obtaining a stability index from a larger number of field trials, as a large number of environments are useful for more reliable estimations of stability [81]. The GS models could be retrained by incorporating data from additional environments in subsequent years. High prediction accuracies for stability could better predict the performance of genotypes at multiple locations during breeding cycles [82].
There was an improvement of prediction accuracy when the model was trained on NAM 175 compared to the NAM 650 population. This is in contrast to other studies where it was observed that prediction accuracies improve when the number of individuals in the training population increases [83][84][85]. This argument was strengthened by the population differentiation coefficient (F st ), suggesting that lines in the NAM 175 were more genetically related compared to the NAM 650 . It has been observed that prediction accuracies increase when training and testing populations are more genetically related to each other, as was the case in the NAM 175 population set [86,87]. These results suggest that selection for lines having GPC stability also aids in the improvement of GS accuracy for GPC. This study opens up an avenue for the utilization of GS in a spring wheat breeding program selecting lines having GPC stability in addition to high GPC.

Conclusions
Selection for GPC is often secondary to grain yield in terms of breeding objectives for spring wheat, although it is a primary trait that producers consider when selecting varieties. We report the first large-scale study of nested association mapping and evaluation of GS for GPC stability in wheat. This study identified wheat lines having less variation for GPC and mapped QTLs controlling GPC stability, an important and often overlooked trait. The identification of stable genotypes with high GPC could help in developing cultivars that can perform similarly in multiple environments. The QTLs identified in this study explained a small amount of phenotypic variation, demonstrating the complexity of the trait, which suggests a GS approach could best address breeding for GPC and GPC stability. Prediction accuracy is sufficiently high for the implementation of GS for GPC and GPC stability in this study. With the implementation of GS for these traits, predictions can be available while evaluating grain yield, enabling selection for GPC along with agronomic and yield traits.