Insights into the Genetic Architecture of Phenotypic Stability Traits in Winter Wheat

Examining the architecture of traits through genomics is necessary to gain a better understanding of the genetic loci affecting important traits to facilitate improvement. Genomewide association study (GWAS) and genomic selection (GS) were implemented for grain yield, heading date, and plant height to gain insights into the genetic complexity of phenotypic stability of traits in a diverse population of US Pacific Northwest winter wheat. Analysis of variance using the Additive Main Effect and Multiplicative Interaction (AMMI) approach revealed significant genotype and genotype by environment interactions. GWAS identified 12 SNP markers distributed across 10 chromosomes affecting variation for both trait and phenotypic stability, indicating potential pleiotropic effects and signifying that similar genetic loci could be associated with different aspects of stability. The lack of stable and major effect loci affecting phenotypic variation supports the complexity of stability of traits. Accuracy of GS was low to moderate, between 0.14 and 0.66, indicating that phenotypic stability is under genetic control. The moderate to high correlation between trait and trait stability suggests the potential of simultaneous selection for trait and trait stability. Our results demonstrate the complex genetic architecture of trait stability and show the potential for improving stability in winter wheat using genomic-assisted approaches.


Introduction
The complex genetic architecture of traits with agronomic and economic importance, such as grain yield, impose challenges in breeding for trait improvement. Complex traits express continuous variation, are affected by many loci with small effects, and are highly influenced by genotype-by-environment interactions (GEI), which refers to the differences in the ability of a genotype to exhibit changes in a particular trait across different environments [1]. A thorough understanding of GEI is a key requirement for progress in any breeding program [2]. When performing phenotypic selection to improve genetic gain, it is necessary to evaluate trait stability and GEI, as such would facilitate breeding of crop varieties that may better cope with changing environments while securing phenotypic stability for quantitative traits such as yield [1,3]. Evaluating trait stability will likely become more relevant as changes in the global climate increase environmental variation [2].
One approach which can facilitate a more accurate estimation for the improvement of trait stability is through examining the genetic architecture underlying stability and implementing genomic-assisted breeding [4,5]. Genetic mapping through association studies and genomic selection has been widely used to dissect the architecture of important traits in wheat (Triticum aestivum L.) [6][7][8][9]. A genomewide association study (GWAS) identifies genetic loci linked to phenotypic values of trait variation through

Experimental Material
A diverse winter wheat panel consisting of 456 lines adapted to Pacific Northwest growing conditions of the US were used in this study. These represent advanced winter wheat breeding lines from public state university breeding programs in the region, including Oregon State University, University of Idaho, and Washington State University, USDA-ARS, and private companies. The majority of these lines were common wheat (284; 62%), whereas the remaining lines were club wheat (172; 38%). This population has been characterized previously for different traits such as grain yield [35], snow mold tolerance [11], eyespot resistance [36], and end-use quality traits [37].

Collection and Analysis of Phenotypic Data
The winter wheat panel was planted at the Washington State University Dryland Research Station in Lind (LND) and the Spillman Agronomy Farm near Pullman (PUL), Washington, between the 2015 and 2019 growing seasons under an augmented design, described previously by Lozada and Carter [35]. Briefly, each block consisted of un-replicated genotypes and repeated checks "Eltan" [38] and "Madsen" [39]. A combination of location and year was regarded to be an environment. In 2016, significant soil crusting delayed emergence of lines in LND, and thus no phenotypic data was collected from the winter wheat panel in this environment.
Grain yield (GY) was evaluated by harvesting whole field plots using a Zurn 150 Plot Combine (Zurn Harvesting, Baden-Württemberg, Germany). Heading date (HD) was recorded as the day when 50% of the heads in a plot emerged from the boot, in Julian Days, and plant height (PH; in cm) was measured as the distance from the ground to the tip of the spike, excluding the awn when present. Phenotypic values of trait adjustments were implemented using the Augmented Complete Block Design in R (ACBD-R) program [40], where best linear unbiased predictors (BLUP) values for combined analysis across environments from raw data were derived based on augmented design [35].
The following model was used for the calculation of BLUE values for individual environments: where Y ij is the phenotype for the i th block and j th individual, µ is the overall mean, B i is the block effect of the i th block; ID, G, and C are the effects of check identifier, genotype, and checks, respectively; ε ij is the residual. For the estimation of BLUP values across all environments, the following model was used: where Y ijk is the phenotypic values of trait value for the i th environment of the j th genotype and k th block, µ is the mean effect; ID, G, C, and , where σ 2 G is the variance due to genotype, n is the number of environments; and σ 2 e is the residual variance. H 2 values for trait stability was calculated from variance components derived from an AMMI approach using a formula previously described by Wang et al. [41].

Genotype by Environment Interactions and Stability Indices
Various parameters for estimating trait stability for GY, HD, and PH were evaluated. These indices included AMMI stability index (ASI), AMMI stability value (ASV), and the Finlay-Wilkinson regression coefficient (FW) implementing a Gibbs sampler where independence of lines and environments was assumed. Two additional stability indices, namely, rank sum (RS) and yield stability index (YSI), were evaluated for GY.
Analysis of genotype by environment interaction through an AMMI approach was conducted using the "agricolae" package [42] in R [43] using the adjusted BLUE values for GY, HD, and PH. AMMI model was in the form where Y ge is the yield of genotype g in environment e; µ is the overall mean; α g is the mean deviation of the genotype; β e is the mean deviation of the environment; λ n is the eigenvalue for the principal component (PC) axis, n; N is the number of PC retained in the model; Υ g and η e are the genotype and environment PC scores for the PC axis, n and θ ge is the residual [44]. The ratio between the sum of squares (SS) due to GE signal (GE signal ) and of genotype (G) was calculated following the methods of Gauch [27] to determine whether implementing an AMMI analysis was relevant for the evaluated traits. AMMI stability index, ASI (D i ) was represented as the sum of the squares scores for the significant IPC and calculated using the following equation from Zhang et al. [31]: the score for the genotype i in IPC, and c is the number of significant IPC. ASV was calculated as follows, as previously described by Purchase et al. [30]: From Equation (4), ASV could be interpreted as the distance between the IPC1 scores against the IPC2 scores from the origin in a two-dimensional plot and could be used to evaluate trait stability after the reduction of noise from the GEI effects [45]. The square root of the product of the ratio between the IPC1 and IPC2 sum of squares and IPC1 score and the IPC2 scores are considered in the calculation. The higher the absolute value for IPC score, the more specifically adapted a genotype is to specific environments, whereas lower values for ASV indicate a more stable genotype across environments [32].
Finlay-Wilkinson (FW) regression coefficients were calculated using the R package "FW" assuming independence of lines and of environments and drawing samples from the model posterior distribution using a Gibbs I sampler [34]. Regression coefficients were calculated using 50,000 iterations and 5000 burn-ins. FW coefficients were derived from the following reaction norm model: where y ijk is the phenotype for the k th replicate of the i th variety in j th environment; µ is the mean; g i is the main effect of the i th variety; b i + 1 is the expected change in the variety performance per unit change of the environment main effect; h j is the main effect of the j th environment; and ε ijk is the residual. In this scenario, the prior distributions for g, b, and h were all multivariate normal [34], and the coefficient b i , i.e., the sensitivity of genotype i to the environment, was regarded as the FW measure of trait stability. Average computation time for FW regression of the 456 wheat lines using an Intel Core i5 3.10 GHz processor and 8.00 GB of RAM was estimated through the "proc.time" function in R.
Rank sum (RS) was calculated as RS = mean rank across environments (R) + standard deviation where R ij is the rank of X ij in the j th environment, R i is the mean rank across environments for the genotype i, and l is the number of environments [32].
Yield stability index (YSI) was measured as YSI = RASV + RY, where RASV is the rank of AMMI stability value, and RY is the rank of the mean grain yield of genotypes across environments [32]. The relationships between trait values, breeding values (GEBV), and stability indices were evaluated using Spearman's rank correlation coefficient (ρ), which is less sensitive to outliers compared with the Pearson correlation coefficient [5].

SNP Genotyping and Genomewide Association Study
SNP genotyping was conducted using the Illumina ® 90K SNP [46] through the USDA-ARS Northern Genotyping Laboratory in Fargo, ND. Markers with minor allele frequency >0.05 and <20% missing data were retained for analyses, resulting in 15,229 SNP markers for GWAS. Out of this number, 12,681 markers (83%) were with known mapped positions in the genome. GWAS was conducted using the Fixed and Random Circulating Probability Unification (FarmCPU) package [47] in R using two different models: (1) considering only the Kinship (K model), and (2) Kinship and the first three principal components (PC) derived from GAPIT [48] included as fixed effects (K-PC model).
A K-PC model with PC = 3 was previously observed to be the most reliable in identifying significant marker-trait associations (MTAs) for yield and agronomic traits in winter wheat [7], and hence, a total of three PC were included in the GWAS model. The number of significant MTAs identified in either or on both GWAS models was also compared. A Benjamini-Hochberg false discovery rate (FDR; [49]) value of 0.05 was used in declaring significant marker-trait associations for trait values and stability indices. Percent variation explained for each significant SNP was estimated by fitting phenotypic and genotypic data in a stepwise regression model and computing the difference in the variance when all the markers are present and when a significant SNP is removed from the model in JMP Pro v.11 [50]. BLUP values of trait values for Pullman (PUL), Lind (LND), and across all locations (ALL) together with the stability indices ASI, ASV, FW, RS, and YSI were used for association mapping. The significant MTAs identified for trait and phenotypic stability were then compared with previous studies in wheat and other crops in terms of phenotypic variation explained and chromosomal locations.

Genomic Selection
Genomic selection was implemented using the ridge regression best linear unbiased prediction (RRBLUP) model [51] under a five-fold cross validation and 10-iteration scheme through the Intelligent Prediction and Association Tool (iPAT) program [52]. RRBLUP uses the function "mixed.solve" and fits the prediction model: where X is a full-rank design matrix for the fixed effects, β; Z is the design matrix for the random effects u; K is a positive semidefinite covariance matrix, obtained from markers using an additive relation matrix function "A.mat"; residuals have a mean of zero and are normal, with constant variance; u and ε being statistically independent [51]. Prediction accuracy for the model was reported as the Pearson correlation between GEBV and BLUP values of trait and trait stability.

Genotype-by-Environment Interactions, Trait Heritability, and Correlations
Analysis of variance using AMMI revealed significant (p < 0.0001) genotype (G) and genotype by environment (GE) interaction effects ( Table 1). The first interaction principal component (IPC1) was responsible for 29% of the variation for GY, whereas it was related to 34% and 33% of the variance for HD and PH, respectively. The first two IPC explained between 47% and 56% of the variation for GY, HD, and PH. There were five significant (p < 0.0001) IPC for HD and PH, and seven for GY (Table  S1). The ratios between the SS of GE signal and G were 1.15 (PH), 1.73 (HD), and 2.55 (GY), indicating that there were substantial G and GE signal effects for the traits, and performing AMMI analysis was appropriate for our datasets. PC biplots of IPC1 from AMMI against adjusted values for GY showed that PUL had higher yield overall than LND ( Figure 1). There was also a clear separation based on trial locations, where the LND and the PUL environments grouped together in different clusters. Winter wheat line entries 119 ("ID581"), 129 ("J010039-1"), 140 ("J950086-002"), and 400 ("X960731") were among the least stable in terms of yield ( Figure 1a) as they were located farther from the origin of the AMMI biplot. J010039-1 (entry 129) had the highest ASV for GY among the lines (6.86), indicating that it has the least stability. Conversely, the most stable genotypes were "SSD02071", "V/W-32", "X960411-1", "A00181", and "ARS010704-1L", based on the ASV derived from AMMI analysis (Table S2). The calculated phenotypic stability indices for the GY, HD, and PH are presented in Table S3. PUL18 and PUL19 environments tend to classify genotypes better for GY, as indicated by the length of their vectors. Genotypes had earlier heading dates in LND compared with PUL environments. The PUL17 and PUL19 environments classified wheat lines better for HD. In contrast to GY and HD, there was no clear separation between trial location for PH based on the PC plot; nonetheless, wheat genotypes were observed to be taller in the PUL environments.
Agronomy 2020, 10, x FOR PEER REVIEW 6 of 16 PC biplots of IPC1 from AMMI against adjusted values for GY showed that PUL had higher yield overall than LND ( Figure 1). There was also a clear separation based on trial locations, where the LND and the PUL environments grouped together in different clusters. Winter wheat line entries 119 ("ID581"), 129 ("J010039-1"), 140 ("J950086-002"), and 400 ("X960731") were among the least stable in terms of yield ( Figure 1a) as they were located farther from the origin of the AMMI biplot. J010039-1 (entry 129) had the highest ASV for GY among the lines (6.86), indicating that it has the least stability. Conversely, the most stable genotypes were "SSD02071", "V/W-32", "X960411-1", "A00181", and "ARS010704-1L", based on the ASV derived from AMMI analysis (Table S2). The calculated phenotypic stability indices for the GY, HD, and PH are presented in Table S3. PUL18 and PUL19 environments tend to classify genotypes better for GY, as indicated by the length of their vectors. Genotypes had earlier heading dates in LND compared with PUL environments. The PUL17 and PUL19 environments classified wheat lines better for HD. In contrast to GY and HD, there was no clear separation between trial location for PH based on the PC plot; nonetheless, wheat genotypes were observed to be taller in the PUL environments. Heritability values for the evaluated traits were moderate to high, where HD was the most heritable trait (H 2 = 0.83), followed by PH (H 2 = 0.79) and GY (H 2 = 0.56) (Table 2). Similarly, H 2 values for trait stability based on ANOVA using an AMMI approach (Table 1) showed moderate to high heritability (between 0.55 and 0.70), where PH was the most heritable trait. Spearman's Rank correlation coefficient (ρ) between BLUP values of trait and breeding values and stability indices were low to high, ranging between −0.92 and 0.84 (Table 3). Negative significant ρ values were observed for grain yield BLUP values and GEBV, and the stability indices, except for FW. Significant correlations were also observed for HD trait values and ASI, ASV, and FW indices. Similarly, for PH, highly significant correlations (p < 0.0001) between BLUP values of trait and breeding values with the stability indices ASI and FW were observed, whereas non-significant correlations were observed for ASV. FW regression coefficient showed a normal distribution, with a center of ~0.25, whereas all the other GY stability indices did not follow this distribution (Figure 2). Several outliers for ASI, ASV, and FW indices were also observed in the boxplots. Computing time for the datasets (456 wheat lines, Heritability values for the evaluated traits were moderate to high, where HD was the most heritable trait (H 2 = 0.83), followed by PH (H 2 = 0.79) and GY (H 2 = 0.56) (Table 2). Similarly, H 2 values for trait stability based on ANOVA using an AMMI approach (Table 1) showed moderate to high heritability (between 0.55 and 0.70), where PH was the most heritable trait. Spearman's Rank correlation coefficient (ρ) between BLUP values of trait and breeding values and stability indices were low to high, ranging between −0.92 and 0.84 (Table 3). Negative significant ρ values were observed for grain yield BLUP values and GEBV, and the stability indices, except for FW. Significant correlations were also observed for HD trait values and ASI, ASV, and FW indices. Similarly, for PH, highly significant correlations (p < 0.0001) between BLUP values of trait and breeding values with the stability indices ASI and FW were observed, whereas non-significant correlations were observed for ASV. FW regression coefficient showed a normal distribution, with a center of~0.25, whereas all the other GY stability indices did not follow this distribution (Figure 2). Several outliers for ASI, ASV, and FW indices were also observed in the boxplots. Computing time for the datasets (456 wheat lines, 4104 observations) was approximately 2 min, on average, comparable to the recorded computation time by Lian and de Los Campos [34] using a Gibbs I sampler.

Marker-Trait Associations for Trait Values and Stability Indices
A GWAS approach identified 146 and 110 significant MTAs for GY, HD, PH, and trait stability indices, respectively (Table S4). Out of these, 12 multi-trait loci controlling both BLUP values and trait stability distributed across 10 chromosomes (Table 4) were identified. These loci were responsible for 0.69% to 11.44% of phenotypic values of trait variation. IWB72787 (Tdurum_contig61410_542) on chromosome 1B was associated with GY, RS, and YSI and related to 2.0% to 8.0% of the phenotypic values of trait variation. Another locus on 7B, IWB7147 (BS00022542_51), was associated with GY, RS, and YSI and explained 3.0% to 10% of the variation for the traits. Four significant SNP markers on chromosomes 2B, 3B, 6A, and 7A were detected for the yield stability traits FW and YSI in both K and K-PC models and were associated with an average of 4.33% variation for the traits. There were 10 and 20 significant SNP markers for the BLUP values of trait and stability indices identified for both GWAS models, respectively. There were 20 SNP markers associated with ASI, 18 for ASV, and 46 for FW for GY, PH, and HD, whereas there were 19 and 12 SNP markers associated with RS and YSI for GY, respectively. Twenty SNP markers associated with either BLUP values of trait or stability indices were with no known mapped position in the genome.

Prediction Accuracy for Trait Values and Stability Indices
Mean genomic prediction ability for BLUP values of trait using an RRBLUP model was 19% higher than for the stability indices (0.43 vs. 0.36; Figure 3). Overall, predicting HD for both BLUP values of trait and stability parameters resulted in the highest mean prediction (0.52), followed by GY (0.35), and PH (0.31). Within the BLUP values of trait, HD_ALL had the highest mean prediction ability (0.70), followed by HD_PUL (0.69), and PH_ALL (0.55). Using BLUP values of trait across all environments resulted in a 60% advantage, on average, compared with using the LND and PUL environments alone for prediction. The FW regression coefficients were the most predictable among the stability indices, with prediction ability of 0.66, 0.49, and 0.46 for HD_FW, GY_FW, and PH_FW, respectively. RS (0.46), on average, was more predictable than YSI for GY (0.41), whereas ASV had higher prediction ability values compared with ASI (0.26 vs. 0.21).

Discussion
Association mapping and GS approaches were used to gain insights into the genetic architecture of trait and stability indices for 456 diverse winter wheat lines evaluated in nine US Pacific Northwest environments. Analysis of GEI using an AMMI approach was found to be suitable for GY, HD, and PH. Genomic regions associated with trait values and stability indices were identified using GWAS, whereas GS showed a low to moderate prediction accuracy for the stability indices. Altogether, results showed a complex genetic architecture for trait stability in winter wheat.
Performing an AMMI to evaluate GEI was appropriate for our datasets, based on a method earlier proposed by Gauch [27], in determining whether performing such analysis to multienvironment trials is advisable. Accordingly, the SS for the signal caused by GEI should be at least as large as that of G, such that there should be a substantial GEsignal, and most of the effects are not buried in noise or residual. In our case, the ratio between the SS of GEsignal and G was between 1.15 (PH) and 2.55 (GY), suggesting that a considerable amount of GEI was due to signal and not to error. Furthermore, significant variation was caused by most, if not all, of the IPC where there were five and seven significant IPC (out of a total of eight) for HD and PH, and GY, respectively. This observation further indicates that GEI was due to signal and not due to error or noise.

Discussion
Association mapping and GS approaches were used to gain insights into the genetic architecture of trait and stability indices for 456 diverse winter wheat lines evaluated in nine US Pacific Northwest environments. Analysis of GEI using an AMMI approach was found to be suitable for GY, HD, and PH. Genomic regions associated with trait values and stability indices were identified using GWAS, whereas GS showed a low to moderate prediction accuracy for the stability indices. Altogether, results showed a complex genetic architecture for trait stability in winter wheat.
Performing an AMMI to evaluate GEI was appropriate for our datasets, based on a method earlier proposed by Gauch [27], in determining whether performing such analysis to multi-environment trials is advisable. Accordingly, the SS for the signal caused by GEI should be at least as large as that of G, such that there should be a substantial GE signal , and most of the effects are not buried in noise or residual. In our case, the ratio between the SS of GE signal and G was between 1.15 (PH) and 2.55 (GY), suggesting that a considerable amount of GEI was due to signal and not to error. Furthermore, significant variation was caused by most, if not all, of the IPC where there were five and seven significant IPC (out of a total of eight) for HD and PH, and GY, respectively. This observation further indicates that GEI was due to signal and not due to error or noise.
Genomewide association studies capture genetic loci linked to significant variation for traits of interest [53,54]. In this study, genomic regions controlling BLUP values of trait and stability indices were identified on 20 out of 21 chromosomes (all excluding 4D), including some that were not mapped to any chromosome. Most, if not all, were minor effect loci, with an average phenotypic variation of 2.70% explained across all BLUP values of trait and stability parameters. Moreover, the majority (89%) of identified loci for trait values and stability were not detected in both K and K-PC GWAS models, indicating the lack of stable loci. Our results, together with previous observations in wheat [3,45] and in rye (Secale cereale) [5,55] demonstrate the complex genetic architecture for trait and phenotypic stability. We identified 12 SNP loci on 10 chromosomes linked to both trait value and stability parameters, demonstrating that these traits are not exclusive and independent, and it is possible to select for both trait and stability indices in Pacific Northwest winter wheat. SNP markers, Tdurum_contig61410_542 (1B) and BS00022542_51 (7B), associated with phenotypic variation for GY, RS, and YSI, represent important loci which could be examined further for marker-assisted breeding as they control multiple yield trait and trait stability measures. The important role of chromosome 5B in controlling GY stability was confirmed from a previous report in wheat substitution lines [45], as most of the significant MTAs controlling yield trait stability were detected in this chromosomal region. We further observed some potential pleiotropic effects as shown by the SNP loci, Excalibur_c53772_302 (5B) and Tdurum_contig29357_338 (6A), associated with GY_ALL, HD_LND, and GY_RS, and GY_ASI and PH_FW, respectively. This was consistent with a previous observation in barley (Hordeum vulgare) where the same QTL regions control different measures of phenotypic stability such as phenotypic variances across environments and FW regression slope for yield and component traits [56]. Altogether, our observations demonstrate that the same genetic loci could affect different aspects of phenotypic values of trait stability in winter wheat. Weak negative relationships between ASI and phenotypic values of trait (BLUP and GEBV) were observed, consistent with previous observations by Huang et al. [3] in soft winter wheat. The moderate to high degree of phenotypic correlation observed between trait values and stability indices demonstrate the ability to simultaneously select for both sets of traits in winter wheat breeding. Genetic loci linked to phenotypic stability and yield performance were also identified on chromosomes 1D and 4B (Table S4), consistent with previous observations in wheat substitution lines [57]. Our results provided additional evidence on the genetic complexity of phenotypic trait stability in wheat. The multi-trait loci related to trait and phenotypic stability identified in this study were not previously reported elsewhere, and hence could represent novel genomic regions affecting variation for the stability of traits in winter wheat.
Genomic-assisted breeding strategies offer the potential for training prediction models enabling the genetic improvement of phenotypic stability traits in different crops [5]. In the current study, moderate to high mean prediction ability values for trait and stability indices were observed. BLUP values of trait were generally more predictable than stability indices, suggesting a more complex architecture for the stability parameters. There was a weak and non-significant correlation (r = −0.08) between prediction ability for the stability of traits and heritability, indicating that other factors aside from H 2 were the main drivers of the observed predictability for the stability indices. Our observations are consistent with another study using the same population of winter wheat lines evaluated for single and multiple trait partial least square regression (PLSR) prediction models, where a weak correlation (r = 0.06) was observed between prediction ability and heritability under cross-validations [35]. In contrast, higher heritability was related to improved predictability for the BLUP values of trait (r = 0.57; p < 0.0001), demonstrating the influence of heritability in the prediction accuracy for the trait values. Recently, it was observed that using multiple secondary correlated traits collected from high-throughput field phenotyping was related to improved accuracy for grain yield in winter wheat [35], and hence could also be used for the indirect selection of lines with high yield potential [58]. Prediction accuracies across environments were also lower compared to predictions within the same environments for yield in winter wheat [35], indicating the relevance of using similar environments to achieve optimal prediction accuracies [3].
The FW regression consistently showed the highest prediction ability (mean of 0.54), whereas AMMI stability parameters (ASI and ASV) had the lowest predictability among the indices for the evaluated traits. This suggests that most of the variation observed for the ASI and ASV may not be under genetic control and that the pattern of GEI was not adequately captured by all IPC used in the model for the AMMI indices [3], even if most of these IPC had significant effects. Only an average of 5% of the total variance was attributed to the IPC, and this could have resulted in lower accuracies for AMMI stability. The FW model measures the stability and/or the responsiveness of a genotype to a specific environment [59]. Previously, a large-scale comparison of 25 genetic improvement studies in barley, cotton (Gossypium hirsutum L.), maize (Zea mays), soybean (Glycine max), and wheat indicated that newer cultivars yielded more than older varieties and had higher FW values. This indicates that breeders have been successful in simultaneously selecting for genotypes that are both higher-yielding and are more specifically adapted to favorable environments across multiple years of selection in major crops [60]. Yield potential and stability among maize hybrid lines were further observed to not be "mutually exclusive" [61], demonstrating the potential for simultaneously selecting for both traits. Improving both stability and average performance of lines nonetheless would warrant selecting for lines that consistently show high performance across all environments [59]. The relatively high predictability values for the FW coefficients in this study demonstrate that phenotypic stability and response of the winter wheat lines are under genetic control, and there is a potential of selecting for stable and responsive genotypes using GS approaches. Furthermore, the significant correlation between yield GEBV and FW coefficients indicates that selecting for lines predicted to yield well could also lead to selecting for cultivars that possess above-average stability. The normal distribution of FW values for the winter wheat panel was in accordance with previous observations for yield in a soybean nested association mapping population (NAM) [62]. Overall, our results suggest that genomewide selection strategies should help facilitate the process of selecting for stable and high-yielding genotypes in winter wheat.

Conclusions
Association mapping and genomic prediction strategies were used to gain insights into the genetic architecture of stability traits for GY, HD, and PH in winter wheat. Evaluating GEI using an AMMI approach was appropriate for the datasets, as most of the effects were not buried in error or noise. Genomewide scans revealed the complexity of stability indices for yield and agronomic traits. Multiple genomic regions on 10 chromosomes were found to control both trait values and stability indices, demonstrating pleiotropic effects or tight linkage and the possibility of simultaneously selecting for both sets of traits. Prediction ability for trait values and stability indices were moderate to high, indicating the potential of selecting for stable and top-performing lines using genomic selection approaches. Altogether, we were able to identify regions controlling trait and trait stability in Pacific Northwest winter wheat, which rendered insights into their complex genetic architecture.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4395/10/3/368/s1, Table S1. Variance, sum of squares, and mean square for the interaction principal components for an AMMI analysis for grain yield, heading date, and plant height in US Pacific Northwest winter wheat.; Table S2. AMMI stability values (ASV) and ASV rank for the US Pacific Northwest winter wheat lines; Table S3. Phenotypic stability indices for grain yield, heading date, and plant height for the US Pacific Northwest winter wheat lines; Table S4. SNP markers associated with trait and trait stability in US Pacific Northwest winter wheat.