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

Abstract

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.

1. 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 the strength of linkage disequilibrium (i.e., the non-random association of alleles at multiple loci) between markers and functional polymorphism across diverse germplasm [10]. In wheat, GWAS has been conducted in different traits such as grain yield, yield components (thousand kernel weight and kernel number), Fusarium head blight resistance, snow mold tolerance, plant height, and heading date [11,12,13,14,15,16]. Compared with biparental mapping, GWAS results in a higher mapping resolution as it exploits the majority of the recombination histories of the individuals belonging to the diverse population used [17,18]. The main drawback of GWAS, nonetheless, is that it might not able to capture the causative and rare variants (loci) with small effects, and therefore would be a special case of missing heritability, i.e., the portion of genetic variance that cannot be explained by all significant loci [19]. In this case, genomic selection (GS) is seen as a complementary approach to GWAS.

GS uses genomewide markers to estimate the marker effects and calculate genomic estimated breeding values (GEBV), which represents the genetic “worth” or merit of an individual. A high GEBV for grain yield, for example, would mean that a line is predicted to have better grain yield compared with others. In contrast to GWAS, GS does not test for significance but rather uses genomewide marker data to calculate the GEBV. More genetic variation could be captured in GS compared with marker-assisted selection as all markers are used in estimating breeding values [20]. In GS, a prediction model is first trained using a training population that has both genotypic and phenotypic information [21]. The model is eventually used to predict the un-phenotyped validation population and calculate GEBV. The correlation between the observed and predicted values is called prediction accuracy and is affected by several factors, including genetic relatedness between the training and test populations, training population size and composition, number of markers, population structure, and genetic architecture of traits [22,23,24,25].

Among the measures of phenotypic stability that have been widely used are the Additive Main Effect and Multiplicative Interaction (AMMI) and Finlay-Wilkinson (FW) regression, where each parameter measures different aspects of trait stability. The AMMI model uses the additive analysis of variance (ANOVA) to partition variation into genotype (G), environment (E), and GEI effects, and the multiplicative principal components analysis (PCA) of the GEI [26,27], where G and G + GE effects represent wide and specific adaptation, respectively [28,29]. AMMI further allows an understanding of complex GEI, including identifying mega-environments, determining lines to exploit narrow adaptations, and predicting performance to improve cultivar recommendations, selections, and genetic gain [7]. Other phenotypic stability parameters were subsequently derived from the AMMI. This includes the AMMI stability value (ASV) based from the AMMI interaction principal components IPC1 and IPC2 earlier proposed by Purchase et al. [30], and the AMMI stability index (ASI) developed by Zhang et al. [31], based on the distance of the IPC with its origin in space. Additional stability measures such as the rank sum (RS) and the yield stability index (YSI) derived from the ASV have also been previously evaluated [32]. Low values for both RS and YSI indicate genotypes with high mean yield and stability. The Finlay-Wilkinson (FW) coefficient measures phenotypic stability and adaptation by regressing to the overall mean of the trait across different environments [33]. FW evaluates how a genotype’s expected performance changes as a function of the environment [34].

The objectives of this study are to (1) identify stable lines in a diverse population of US Pacific Northwest winter wheat lines evaluated under multi-environment trials for yield, heading date, and plant height based on phenotypic stability indices; (2) determine genomic regions controlling phenotypic stability of traits; (3) evaluate genomic selection accuracy for trait values and phenotypic stability indices; and (4) identify the relationship between genomic breeding values and trait stability indices. Different measures of phenotypic values of trait stability were calculated, including ASV, ASI, FW regression coefficient, RS, and YSI.

2. Materials and Methods

2.1. 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].

2.2. 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:

Y_ij = µ + B_i + ID + G + C + ε_ij

(1)

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:

Y_ijk = µ + ID + G + C + E_i + E_i × ID + G × E_i + C × E_i + B_k(E_i) + ε_ijk

(2)

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 E_i are the effects of the check identifier, genotype, check, and i^th environment, respectively; check ID was used to distinguish one check over the other checks (e.g., check 1 vs. check 2 vs. check 3) used in the analyses, as well as the number of checks present on each block; E_i × ID is the effect of environment by check identifier interaction; G × E_i is the effect of genotype by environment interaction; C × E_i is the effect of check by environment interaction; B_k(E_i) is the effect of the k^th block nested within i^th environment, and ε_ijk is the residual. Effects were considered fixed when calculating BLUE and as random when estimating BLUP values.

Broad-sense heritability (H²) for combined analysis across environments was calculated using the formula $H^2=\frac{{\sigma }_G^2}{{\sigma }_G^2+\raisebox{1ex}{${\sigma }_G^2$}\!\left/ \!\raisebox{-1ex}{$n$}\right.+{\sigma }_e^2},$ where ${\sigma }_G^2$ is the variance due to genotype, n is the number of environments; and ${\sigma }_e^2$ is the residual variance. H² values for trait stability was calculated from variance components derived from an AMMI approach using a formula previously described by Wang et al. [41].

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

$${\text{Y}}_{ge}=\mathsf{\mu }+{\mathsf{\alpha }}_g+{\mathsf{\beta }}_e+{\displaystyle \sum }_{n=1}^N{\mathsf{\lambda }}_n{\mathsf{{\rm Y}}}_g{\mathsf{\eta }}_e+{\mathsf{\theta }}_{ge}$$

(3)

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; ${\mathsf{{\rm Y}}}_g$ and ${\mathsf{\eta }}_e$ are the genotype and environment PC scores for the PC axis, n and ${\mathsf{\theta }}_{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]: D_i = $\sqrt{{{\displaystyle \sum }}_{s=l}^c{Y}_{is}^2}$, where Y_is² is 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]:

$$\mathrm{ASV}=\sqrt{ {\left[\frac{IPC1sum of square}{IPC2sum of square}\left(IPC1score\right)\right]}^2+ {\left(IPC2 score\right)}^2 }$$

(4)

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:

y_ijk = µ + g_i + (b_i + 1)h_j + ε_ijk

(5)

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 of rank (SDR), where SDR = $\sqrt{\frac{{{\displaystyle \sum }}_{j = i}^m\left(R_{ij} - R_i\right)2}{l - 1}}$, 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].

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

2.5. 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:

y = Xβ + Zu + ε

(6)

u ~ N (0, Kσ²_u)

(7)

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.

3. Results

3.1. 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. Variance components for grain yield, heading date, and plant height estimated using an Additive Main Effect and Multiplicative Interaction (AMMI) approach.

Trait	Genotype (G)	Environment (E)	GE	Residuals	No. of Sig. IPC	% Variation for IPC1	Ratio (GEsignal/G)
Grain yield	3.90	196.61	15.65	12.78	7	29.30	2.55
Heading date	44.31	1753.47	126.23	112.09	5	33.70	1.73
Plant height	474.09	8725.56	948.72	909.79	5	32.60	1.15

All genotype and genotype by environment interactions are significant at p < 0.0001.

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

Heritability values for the evaluated traits were moderate to high, where HD was the most heritable trait (H² = 0.83), followed by PH (H² = 0.79) and GY (H² = 0.56) (

Table 2. Heritability and mean of trait values and stability indices for grain yield, heading date, and plant height.

Trait 1	Heritability	Heritability, Trait Stability 2	Mean (Adj. Values)	Stability Index (Mean)
				ASI 3	ASV	FW	RS	YSI
GY (t ha−1)	0.56	0.55	6.10	2.10	1.95	0.02	350.85	457
HD (Julian Days)	0.83	0.63	151.74	5.71	7.11	0.07	-	-
PH (cm)	0.79	0.70	87.1	15.43	15.59	0.15	-	-

¹ GY, grain yield; HD, heading date; PH, plant height. ² Heritability estimates for trait stability based on using ANOVA derived from AMMI analysis (Table 1). ³ ASI, AMMI stability index; ASV, AMMI stability value; FW, Finlay-Wilkinson regression; RS, rank sum; YSI, yield stability index.

). Similarly, H² 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. Spearman’s Rank correlation coefficient (ρ) between trait values and stability indices. Trait values were represented as adjusted and genomic estimated breeding values (GEBV).

Stability Index	Grain Yield		Heading Date		Plant Height
	Adjusted Values	GEBV	Adjusted Values	GEBV	Adjusted Values	GEBV
ASI	−0.11 *	−0.11 *	−0.12 *	−0.12 *	0.05	0.01
ASV	−0.05	−0.08	−0.10 *	−0.10 *	−0.02	−0.04
FW	0.41 ***	0.35 ***	0.74 ***	0.75 ***	0.84 ***	0.77 ***
RS	−0.92 ***	−0.81 ***	-	-	-	-
YSI	−0.66 ***	−0.61 ***	-	-	-	-

* Significant at p < 0.05, ** significant at p < 0.001, *** significant at p < 0.0001. ASI, AMMI stability index; ASV, AMMI stability value; FW, Finlay-Wilkinson regression; RS, rank sum; YSI, yield stability index.

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

3.2. 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. Multi-trait loci controlling trait and stability for grain yield, heading date, and plant height in Pacific Northwest winter wheat.

SNP ID	SNP Name	Chr.	Position (cM)	Allele	Minor Allele Frequency	Positive Allele 1	Trait 2	Model	Effect	FDR-Adjusted p-Value	R2
IWB9588	BS00065601_51	1A	84.63	[A/G]	0.17	A	PH_ALL	K	1.18	2.2E-04	0.02
						A	PH_ALL	K-PC	1.15	4.08E-05	0.03
						A	PH_FW	K	0.03	4.0E-04	0.04
IWB72787	Tdurum_contig61410_542	1B	53.49	[A/G]	0.28	G	GY_ALL	K-PC	0.06	0.003	0.03
						A	GY_YSI	K	−61.66	7.8E-04	0.08
						A	GY_RS	K-PC	−44.14	0.003	0.06
						A	GY_YSI	K-PC	−79.59	2.98E-06	0.02
IWB7202	BS00022666_51	2A	115.08	[A/G]	0.27	G	HD_PUL	K	0.34	1.14E-03	0.03
						G	HD_FW	K	0.022	0.007	0.01
IWB72096	Tdurum_contig50589_867	2B	114.09	[T/C]	0.27	T	HD_PUL	K	−0.29	1.14E-05	0.02
						T	HD_PUL	K-PC	−0.22	0.015	0.02
						T	HD_FW	K-PC	−0.02	0.010	0.01
IWB7001	BS00022276_51	2D	18.22	[T/G]	0.20	G	HD_ALL	K-PC	0.27	4.21E-04	0.07
						G	HD_PUL	K-PC	0.37	5.08E-04	0.06
						G	HD_FW	K	0.03	0.006	0.05
						G	HD_FW	K-PC	0.03	3.29E-04	0.11
IWB75368	wsnp_BE405275A_Ta_1_1	4A	29.86	[T/G]	0.09	T	GY_ALL	K	−0.08	3.05E-05	0.02
						T	GY_ALL	K-PC	−0.07	0.002	0.02
						G	GY_RS	K-PC	37.42	0.037	0.02
IWB12343	BS00108019_51	5B	51.16	[T/C]	0.47	T	PH_ALL	K	−0.63	0.003	0.02
						T	PH_ALL	K-PC	−0.56	0.029	0.05
						T	PH_FW	K	−0.02	0.003	0.03
IWB27416	Excalibur_c53772_302	5B	68.36	[T/C]	0.05	C	GY_ALL	K-PC	0.06	0.032	0.01
						T	HD_LND	K-PC	−0.19	0.012	0.01
						T	GY_RS	K-PC	−55.62	0.001	0.03
IWB69770	Tdurum_contig29357_338	6A	118.07	[A/C]	0.22	A	GY_ASI	K	−0.20	0.005	0.02
						A	PH_FW	K	−0.02	0.036	0.06
IWB77743	wsnp_Ex_c31955_40681185	7A	216.28	[T/C]	0.25	T	GY_ALL	K	−0.04	0.010	0.02
						C	GY_YSI	K	38.92	0.004	0.01
IWB80605	wsnp_Ku_c5693_10079278	7A	208.71	[A/C]	0.28	A	PH_ALL	K	−0.59	0.028	0.01
						A	PH_FW	K	−0.02	0.001	0.04
IWB7147	BS00022542_51	7B	76.31	[T/C]	0.49	T	GY_ALL	K-PC	−0.04	0.011	0.03
						C	GY_RS	K	28.88	0.004	0.10
						C	GY_YSI	K	42.28	0.003	0.05
						C	GY_RS	K-PC	32.22	0.001	0.06

¹ Positive allele is the allele that increases trait value for trait and phenotypic stability. ² GY, grain yield; HD, heading date; PH, plant height; ALL, BLUP values of trait across all environments; ASI, AMMI stability index; FW, Finlay-Wilkinson regression; RS, rank sum; YSI, yield stability index.

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

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

4. 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² 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.

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