Multi-model genome-wide association and genomic prediction analysis of 16 agronomic, physiological and quality related traits in ICARDA spring wheat

Identification and exploration of the genetic architecture of traits related to yield, quality, and drought and heat tolerance is important for yield and quality improvement of wheat through marker-assisted selection. One hundred and ninety-two spring wheat genotypes were tested at two heat-stress locations in Sudan (Wad Medani and Dongula), a drought stress site in Morocco (Marchouch) and a site with high yield potential in Egypt (Sids) in replicated trials during the 2015–2016 and 2016–2017 cropping seasons. A total of 10,577 single nucleotide polymorphism markers identified from the 15 K wheat SNP assay were used in a genome-wide association (GWA) study and genomic prediction for 16 phenotypic traits related to yield, quality and drought and heat tolerance. Significant marker-trait associations were detected across GWAS models for all traits. Most detected marker-trait associations (MTAs) were environment-specific, signifying the presence of high quantitative trait loci-by-environment (QTL x E) interaction. Chromosome arm 5AL had significant multi-model MTAs for grain yield and yield-related traits at the heat-stress locations. Highly significant QTLs were detected on chromosome 2D for waxiness. Homoeologous group 2 and 6 chromosomes were with significant MTAs for grain protein content, gluten content, alveograph strength and Zeleny sedimentation test while chromosome arm 3BL was significant for both Z and W traits. Genomic prediction analysis with ridge regression-best linear unbiased prediction model estimated the breeding values of the studied traits with prediction accuracies ranging from 0.16 for leaf rolling to 0.72 for peduncle length. The identified QTLs could be targeted for marker-assisted selection or further studies aimed at fine mapping and cloning the causative genes and detecting favorable haplotypes with positive effects for agronomic, physiological or quality-related traits.


Introduction
Bread wheat (Triticum aestivum L., 2n = 6x = 42), as one of the major food crops, providing approximately 19% of calory and 21% of protein needs of human populations (Braun et al. 2010). Meeting the demand of wheat production for the ever-increasing world population is a challenging task.
The challenge of meeting the increasing demand for wheat in developing countries is increasing due to rapid population growth coupled with the effects of climate change (Lobell et al. 2011). Hence, breeders actively develop varieties with higher yield and yield stability, resistance to biotic and abiotic stresses, and quality parameters that fulfill the requirements for specific end-products.
Grain yield is an outcome from an aggregation of several agronomic and physiological traits (Chen et al. 2012;Sukumaran et al. 2018;Tshikunde et al. 2019;Jamil et al. 2019). Agronomic traits for attention include grain number per spike, thousand kernel weight, above-ground biomass, plant height, whereas important phenological traits include heading and maturity dates as the grain filling period under specific seasonal conditions has an enormous effect on final grain yield (Reif et al. 2011b;Ihsan et al. 2016;Sun et al. 2017;Wang et al. 2017;Ma et al. 2018;Liu et al. 2018;Li et al. 2020). Drought and heat stresses are the two most important environmental constraints that curtail wheat production globally. These stresses cause various physiological changes in plants leading to accelerated growth and premature senescence which in turn reduce yield potential. Various physiological traits, including chlorophyll content, canopy temperature, stomatal conductance, waxiness (glaucousness), and leaf rolling have been targeted as adaptive traits in abiotically stressed environments and it is important to understand their roles in grain yield performance in heat/drought affected environments (O'Toole et al. 1979;Reynolds et al. 2007; Reynolds and Tuberosa 2008;Saint Pierre et al. 2010;Paliwal et al. 2012;Cossani and Reynolds 2012;Adamski et al. 2013;Guo et al. 2016). In addition to yield and yield-related agronomic and physiologic traits, grain quality traits are also crucial since they determine the unique characteristics and value of a wide range of end-use products such as pan bread, cookies, cakes, and pastries. Several physical and chemical properties contribute to the quality of bread wheat. Quality parameters determined by grain protein content, gluten content, dough rheological properties, and Zeleny sedimentation test are targeted for wheat quality (Zanetti et al. 2001;Groos et al. 2004;Reif et al. 2011a;Carter et al. 2012;Würschum et al. 2016;Battenfield et al. 2018;Guo et al. 2020).
With the surge in high throughput sequencing technologies, GWAS has become a widely used approach to identify genes and dissect QTLs in plants for almost all traits (Brachi et al. 2011). GWAS overcomes the two common limitations (i.e. restricted allelic diversity and limited genomic resolution) of the bi-parental QTL mapping approach (Brachi et al. 2011). The main challenge for GWAS is control of false positives caused by population structure and family relatedness (Kaler et al. 2020). This problem has been addressed by incorporating these two confounding factors as covariates using the mixed linear model (MLM) (Price et al. 2006). However, due to overfitting, this model frequently leads to falsenegatives that might exclude key loci. Multi-locus GWAS analysis methodologies such as multiple locus mixed linear model (MLMM), fixed and random model circulating probability unification (FarmCPU) and Bayesian-information and linkage-disequilibrium iteratively nested keyway (BLINK) have emerged to overcome false-negatives (Zhang et al. 2019). In MLMM, marker association tests are used to select associated markers that are fitted as cofactors for marker tests. The cofactors are then adjusted through forward and backward stepwise regression in the mixed model (Segura et al. 2012). FarmCPU is a more recently developed multi-locus GWAS model with robust tradeoffs for controlling false-positives without compromising false negatives (Liu et al. 2016). FarmCPU uses the MLMM model and incorporates multiple markers simultaneously as covariates in a stepwise MLM to partially remove the confounding effect between tested markers and kinship. This model uses both the fixed-effect model (FEM) and random effect model (REM) iteratively to eliminate confounding factors. The FEM encompasses testing markers and multiple associated markers as covariates to control false positives whereas REM estimates kinship using associated markers that prevent model overfitting (Liu et al. 2016). However, the REM model in FarmCPU demands high computing effort for large numbers of individuals. A new multi-locus model known as BLINK was developed to overcome this limitation (Huang et al. 2019). This model approximates maximum likelihood using Bayesian information criterion (BIC) in a FEM model. BLINK also removes the assumption that causal genes are evenly distributed across the genome as required by the SUPER and FarmCPU methods, making the model superior in statistical power with discovery of less false-positives (Huang et al. 2019). Genomic prediction is effectively used in several crop plants to predict the genotypic estimated breeding values (GEBVs) of individuals solely based on their overall genome-wide marker effects without the need of phenotyping (Meuwissen et al. 2001;Crossa et al. 2017). This method overcomes the limitation of GWAS or other linkage-based mapping methods by including QTLs with minor effects (i.e., QTLs below a specified significance threshold) and has been used in several plant breeding programs (Crossa et al. 2017).
The objective of this study was to identify QTLs for 16 agronomic, physiologic and quality traits and to predict the genomic breeding values of these traits. To accomplish this, 192 spring bread wheat genotypes was assembled from International Center for Agricultural Research in the Dry Areas (ICARDA) materials and grown at four different locations in Egypt, Morocco and Sudan with contrasting environmental conditions for two consecutive years.

Plant materials
A panel of 192 spring bread wheat genotypes from ICARDA were assembled for this study. The pedigrees and selection histories of these genotypes are provided in Table S1. This panel comprised a set of synthetic derivatives, cultivars from Central and West Asia and North Africa (CWANA) region, and elite breeding lines from ICARDA's shuttle breeding program (Tadesse et al. 2016).

Phenotyping
Field evaluations were carried out at four locations including two sites in Sudan (Dongula and Wad Medani), Sids in Egypt and Marchouch in Morocco for two years (2015-16 and 2016-17). Wad Medani station is located at 14°24 0 04 00 North and 33°31 0 11 00 East at 410 masl and is a global field trial platform for heat tolerance studies managed by a collaboration of the Agricultural Research Corporation of Sudan, International Maize and Wheat Improvement Center (CIMMYT) and ICARDA. The annual temperature at this site is in between 22 -40°C with a mean of 37°C during March which is the flowering and grain-filling stage of wheat at this location. The wheat cropping cycle is short (December -March). The soil is a calcareous Vertisol with pH 8.5 and is deficient in N (300 mg L-1) and P (4 -6 mg L-1) (Tadesse et al. 2019). The Dongula station is located at 19°10' North and 30°29' East at 226 masl. Both Wad Medani and Dongula stations are irrigated. Sids station in Egypt is located at 29°3 0 58.06 00 North and 31°5 0 57.79 00 East at 322 masl. The soil is a highly fertile clay loam with pH 7.8 and is a high-yielding site for wheat production. The temperature ranges from 20 to 35°C with an average 30°C during the cropping cycle (from December to April). The Marchouch station in Morocco, located at 33.6°N and 6.7°W at 410 masl, is rainfed with annual average rainfall of 350 mm and temperature range of 10 to 23°C with an average of 16.6°C. The soil is a Chromic Luvisol with a surface pH of 6.
Genotypes were planted at both stations in Sudan and in Egypt at the first week of December in both years under irrigated conditions. Trials were irrigated by flood irrigation applied at regular intervals, with seven events at Dongula and Wad Medani and nine at Sids. Planting at Marchouch was in mid-December under rainfed conditions. The experiment was conducted using an alphalattice design with two replications and 20 incomplete blocks each comprising 10 genotypes with a plot size of 2.5 m length and six rows with 0.2 m spacing.
Seeding rate was 100 kg/ha. Local management practices were applied at each location (Tadesse et al. 2019). The phenological observations of days to heading (DH) and days to maturity (DM) were recorded when 50% and 90% of plants in each plot reached heading and maturity, respectively. The grain filling period (GF) was recorded as days taken from DH to DM. Above-ground biomass (BM), number of grains per spike (GRSP), plant height (PH), and peduncle length (PL) were determined in a 1 m 2 area from the central part of each plot at the time of harvesting. Plant height (PH) was recorded as the distance from the ground to the tip of the spike excluding awns. Grain yield (GY) was recorded for each plot and thousand-kernel weight (TKW) was measured by weighing 1000 grains per plot.
Canopy temperature (°C) was evaluated using a handheld infrared IR 1000 thermometer (Fluke 566 series, Everett, Washington, USA). The readings included average leaf temperature of one linear meter in each plot at the vegetative (CTD1) and at grain filling (CTD2) stages. Leaf waxiness (WX) and leaf rolling (L) were based on observations carried out at optimal growth stages with appropriate scales according to guidelines listed by Torres and Pietragalla (2012).
Grain protein content (PR) was estimated by a nearinfrared transmittance spectrophotometer according to the approved method of the American Association of Cereal Chemists (AACC) (ICC,Vienna,Austria). Grain gluten content (GL) was measured with an automatic gluten measuring apparatus (GlutomaticÒ) according to AACC approved method 38-12.02. Dough rheological behavior was measured by alveograph (Chopin S.A., Villeneuve la Garenne, France) strength (W) using a white flour sample (ICC standard method No. 122) and sedimentation value (Z) was determined following the Zeleny method (ICC standard 116/1 (ICC 2008)).

Data analysis
A multi-environment trial analysis with the R (META-R) v6.04 (Alvarado et al. 2020) software package was used to conduct analysis of variance (ANOVA) and to estimate best linear unbiased estimations (BLUEs) of phenotypic traits including years, replications, blocks and genotypes for all test sites. The repeatability of traits for each environment was calculated according to the following formula: Broad-sense heritability for the combined analysis was calculated using the formula: where r 2 G and r 2 E are the genotype and error variance components, respectively, nRep is the number of replicates, r 2 GE is the G 9 E interaction variance component and nloc is the number of environments. Frequency distributions of the phenotypic data were analyzed using Minitab 18 (Minitab Ltd., Coventry, UK) and Pearson's correlations between phenotypic traits were computed using the cor function in R environment (R Core Team 2020).

Genotyping
Leaf samples from 5 two-week-old seedlings from the same seed source as used for field trials were pooled for genomic DNA extraction. The pooled leaf tissue samples were frozen in liquid nitrogen and stored at -80°C before DNA extraction according to Ogbonnaya et al. (2001). The extracted genomic DNA was sent to TraitGenetics (Gatersleben) for genotyping with the 15 K SNP HD Custom Genotyping BeadChip, which is an optimized and reduced version of the 90 K wheat SNP assay described by Wang et al. (2014b). SNP markers with \ 0.05 of minor allele frequency (MAF) and [ 0.1 missing values per genotype were excluded. After quality assessment, 10,577 SNPs comprising 3,403 from the A-genome, 4,533 from the B-genome, 954 SNPs from the D-genome, and 1,582 SNPs with unassigned chromosomes were included for GWAS study. Chromosome positions of SNPs were assigned based on the consensus map of wheat generated by Wang et al. (2014b).

Genome-wide association analysis
Genome-wide association analysis was performed using six models encompassing single-and multi-locus analysis methods using the GAPIT v.3 software package in R (Lipka et al. 2012). Population structure and linkage disequilibrium of the panel were reported in a previous study (Tadesse et al. 2019). The singlelocus methods including the general linear (GLM) and mixed linear (MLM) models were implemented in which the first included only population structure (Q) modeled using principal component analysis according to Price et al. (2006). MLM on the other hand included both the kinship similarity matrix (K) and Q as random and fixed effects, respectively, (Yu et al. 2006) making the model more stringent and able to control false-positive MTAs. Settlement of MLM under a progressively exclusive relationship (SUPER) (Wang et al. 2014a) was the other single-locus analysis method used in this study. Multi-locus based GWAS models including MLMM (Segura et al. 2012), FarmCPU (Liu et al. 2016) and BLINK (Huang et al. 2019) were applied for MTA analysis. Days to heading was incorporated as a covariate in all models to avoid confounding effects from this phenological factor. Hence, GWAS was executed using the data set generated from the BLUEs mean values of phenotypic traits and 10,577 SNP markers with six GWAS analysis models. Population structure and the kinship similarity matrix were also incorporated according to the model type. The exploratory threshold with P B 0.001 (Àlog 10 PC3) was used to report significant marker-trait associations.

Genomic prediction and cross validation
The ridge regression-best linear unbiased prediction (RR-BLUP) mixed model was utilized to estimate the GEBVs of 192 genotypes for the 16 studied traits using the R package ''rrBLUP'' (Endelman 2011; Endelman and Jannink 2012) with the formula: where X and Z were the design matrices and b and l represented the vectors for fixed and random effects, respectively; y was a vector of phenotypic values and e was the residual error. To evaluate the prediction accuracy, cross-validation was estimated by allocating genotypes randomly as training and testing sets with a size of 80 and 20%, respectively, and repeated in 500 iterations. Predictive ability was estimated from the correlation coefficient between the observed phenotypic values and predicted values of the test set. Prediction accuracy was then calculated from prediction ability divided by the square root of repeatability of the trait (Legarra et al. 2008;Chen et al. 2011).

Agronomic traits
The BLUEs for the 16 traits tested at four locations in two consecutive years are given in Table S2 and the phenotypic variation within the four locations and combined years is summarized in Tables S3, S4 and S5. The analysis of variance detected highly significant variation among genotypes for all traits, except CTD1, across locations and seasons. In the year 2015-16, grain yield ranged from 1.8 to 4.4, 0.7 to 4.1, 7 to 13 and 2.5 to 8.9 t/ha with mean values of 2.8, 2.3, 10 and 6 t/ha at Wad Medani, Dongula, Sids and Marchouch, respectively. The corresponding values in 2016-17 were 1.6 to 4.3, 0.8 to 5.9, 4.5 to 13 and 2.2 to 5.9 t/ha with average values of 3, 3.5, 9.7 and 3.6. A wide range of repeatability values was recorded for the tested traits. For instance, there was high repeatability for grain quality traits (Z (0.9), W (0.81), PR (0.81), and GL (0.82)); yield and yield-related traits (DH (0.91), DM (0.88), GY (0.8), PH (0.74), TKW (0.76), and GFP (0.73)) and drought-related traits (WX (0.8), and L (0.8)) for the trials at Wad Medani in 2016 (Table S4). Frequency distributions for all recorded traits across the four locations and two years are given in Fig. S1. Normal distributions were obtained for all traits across environments and years, except for leaf rolling (L) at Sids in both years where the distribution was bimodal. The box plot distribution clearly shows that the Sids location outperformed others in yield and yield-related traits whereas Wad Medani was a better performing for grain quality traits than Dongula and Marchouch (Fig. S2).

Correlations among traits
Pearson's correlations among agronomic, physiological and quality traits tested across locations and years are presented in Figs. S3A, S3B, S3C and S3D. Grain yield was significantly correlated to BM, DH, PH, DH, GRSP, and WX but negatively correlated to GFP, PR, GL, W, Z, CTD1 and CTD2 at Wad Medani ( Fig. S3A). A positive correlation was recorded between grain quality traits (PR, GL, W and Z), but these traits were negatively correlated to GRSP, BM, GY and WX at Dongula (Fig. S3B). Days to heading (DH) was negatively correlated to PR, GL and W at Marchouch (Fig. S3C). Grain filling period (GFP) showed a highly significant negative correlation with DH but a positive correlation with TKW at Sids (Fig. S3D).

Genome-wide association analysis
Significant MTAs were detected for the 16 traits. Six different GWAS models comprising both single-locus (GLM, MLM and SUPER) and multi-locus (MLMM, FarmCPU and BLINK) methods were implemented to detect reliable MTAs. All identified MTAs are listed in Tables S6-S11 and the Manhattan and quantilequantile (Q-Q) plots are given in Figs. S6-S11. The Q-Q plots show that models MLM, MLMM, Farm-CPU and SUPER controlled false positives quite well (Figs. S7-S10).

GWAS for yield and yield-related traits
Chromosome regions were identified with significant MTAs for eight yield and yield-related traits using the six GWAS models. Highly significant multi-model MTAs were detected on the long arm of chromosome 5A for GY at the two heat stressed locations in Sudan. SNP marker Ku_c19858_2078 on chromosome arm 5AL revealed highly significant MTA with GY at Wad Medani in 2016 -17 across models (Table 1). Several SNP markers within a 3 cM interval uncovered significant multiple MTAs for this trait (Fig. 1). In both years, all models detected significant multiple MTAs on arm 5AL for GY tested at Wad Medani and Dongula. SNP markers RAC875_c38693_319 and RAC875_c34939_467 on chromosome arm 7BL showed an environment-specific, highly significant MTA with GY at the drought-prone Marchouch site. Five models, except SUPER, identified marker BS00057445_51 on the short arm of chromosome 3A as a significant MTA with GY at the high yielding Sids location. Regions in chromosomes 1A, 2B, 6B, and 7B showed highly significant environment-specific MTAs for GY across the six models (Figs. S6-S11). SNP marker wsnp_Ex_c17303_25979191 on chromosome arm 3BL revealed a significant multi-model MTA for TKW at Sids (Table 1). Excal-ibuR_c6255_1119 on chromosome arm 1AL identified a significant MTA for TKW with all six models at Dongula. Markers on other chromosomes, including 5B (BS00024993_51), 6A (KukRi_Rep_c104648_106) and 6B (ExcalibuR_c99745_169) also revealed environment-specific significant MTAs for TKW across models. Various environments and model-specific MTAs for TKW were also detected (Tables S6-S11). SNP marker wsnp_RFL_Contig2403_1927045 on the long arm of chromosome 1B displayed a significant MTA with GRPS at Dongula. Environment-specific MTAs were detected for GRPS on chromosomes 1A, 3A, 3D, 5A and 6A across the different test models. In addition to GY, the SNP marker Ku_c19858_2078 on chromosome arm 5AL had a concurrent effect on DM, BM, GFP and PH across environments.
Marker Ku_c19858_2078 that is very close to BS00022071_51 had a highly significant association for DM at Wad Medani in 2016 -17. This SNP also showed significant MTA for the same trait at Dongula. Multi-model significant MTAs were also detected for DM on chromosomes 2A, 6A, 6D, 7A and 7D across environments. SNP marker TduRuM_contig14130_315 on chromosome arm 5BS showed a significant MTA with GFP at Dongula. RAC875_c29314_291 on chromosome arm 7DL and wsnp_Ex_c10084_16572374 on arm 2AS were other SNP markers that exhibited robust MTAs across the six models for GFP tested at Marchouch. Various MTAs were identified for aboveground biomass (BM) such as CAP8_c710_140, BS00022071_51, RAC875_c17347_312, and KukRi_R-ep_c70199_506 on chromosomes 4A, 5A, 6B and 7A, respectively, across the six models.
Highly significant MTAs for PH across models and environments were detected on the long arm of chromosome 5B. The SUPER model, in particular, identified 195 significant multiple MTAs for PH within a 4 cM interval on chromosome arm 5BL across the four environments (Table S10). Chromosomes 1A, 2D, 4A, 5A and 7B revealed significant MTAs for PH across environments and models. The BLINK model identified a highly significant MTA with marker TA001068-0306-w on chromosome arm 3AL for PL with LOD of 13.0 at Wad Medani. SNP  in 2016-2017 (Table S10). All six models detected environment-specific significant multi-MTAs for CTD2 on chromosome arm 1BL; e.g., ExcalibuR_c23992_436.
All six models identified significant MTAs on chromosomes 2B, 2D, 5B, 6B, and 7B for WX. Marker RFL_Contig4849_702 on chromosome arm 2BL was significantly associated with WX at Wad Medani across the six models (Fig. 2), whereas Excal-ibuR_Rep_c109101_94 on chromosome arm 2DS was associated with WX at Sids across the six models. Both environment and model-specific significant MTAs for WX were identified (Tables S6-S11).
Leaf rolling is a well-established drought/heatrelated trait. Various significant MTAs were detected for this trait across environments and within years using the six GWAS models. Multi-model MTAs with leaf rolling were identified on chromosomes 5A, 1B, 4B, 6A and 6B (Table 3). For instance, marker Ku_c47168_563 on the long arm of chromosome 5AL displayed a consistently significant MTA across all six models.

GWAS for quality-related traits
Highly significant and stable MTAs were detected for the four bread wheat quality traits. Significant multi-MTAs were identified on homoeologous group 2 and 6 chromosomes for PR and GL (Tables S6-S11). For instance, all models identified wsnp_Ex_c1988_ 3742022 on chromosome arm 6DS with significant MTA for PR at Wad Medani (Table 4). Markers IAAV8501 and GENE-3709_393 on unassigned chromosomes were consistently associated with PR at Wad Medani and Marchouch in 2016-2017, respectively. Chromosomes 2D, 5A, 6A and 6D had stable associations with GL across models. Chromosomally unassigned marker IAAV8501 was associated with GL at Marchouch.
Marker RAC875_c24504_119 on chromosome arm 3BL was significantly associated with alveograph strength at Dongula (Fig. 3). Significant multi-MTAs were detected for W on chromosomes 2A (RAC875_c26214_505, KukRi_c3882_2021 and CAP12_c259_307) and 2D (BS00093111_51, BS00022532_51 and TA009010-0422) across all models at Marchouch (Tables S6-S11). Marker BS00065960_51 on chromosome arm 6BS also showed a significant MTA at the same location. IAAV8501 with an unassigned chromosome position also had a stable significant MTA with W across the six GWAS models at Marchouch.
In addition to alveograph strength, the SNP Marker RAC875_c24504_119 was associated with the Zeleny sedimentation test at Dongula across all six models (Table 4). RAC875_Rep_c118667_79 on chromosome arm 2BS and CAP8_c3568_256 on chromosome arm 3AL were associated with Z evaluated in Dongula across all models (Fig. 3). Several significant markertrait associations were detected for Z with the six models across the four environments tested for two years (Tables S6-S11).

Traits co-localized on chromosome regions
SNP markers with multiple trait effects were identified in numerous chromosomes regions. Several SNP markers within a 4 cM interval on chromosome arm 5AL exhibited MTAs with possibly pleiotropic effects on yield and yield-related traits (Fig. 4). For example, Ku_c19858_2078 on chromosome arm 5AL was one such marker associated with multiple phenotypic effects involving GY, PH, GFP, DM and BM detected in all four environments with all six models (Fig. S4). Except for BLINK and FarmCPU, four models identified an environment-specific significant MTA between Ku_c19858_2078 and CTD2 at Wad Medani. This marker also had an environment-specific significant MTA with quality-related traits, protein and gluten content at Wad Medani. Except for FarmCPU, the other models detected significant multiple MTAs on chromosome arm 1AL for both TKW and PH across environments and years. Chromosome arm 5BL had multiple SNPs with MTAs for PH and PL across environments and years. The model SUPER detected more than 140 highly significant MTAs on chromosome 5BL (within 4 cM), mainly for both PH and PL but also for GFP and L (Table S10). In both years, chromosome arm 6AL showed significant multi-MTAs for phenological traits (GFP and DM) at Sids, Wad Medani and Dongula. BS00093111_51 on the short arm of chromosome 2D exhibited concurrent effects with W, GL and PR across environments. The MLM model detected significant MTAs on 2AS with concurrent effects on quality traits (PR, GL, W and Z) at the same three environments  Table 2 and footnote to Table 1 for abbreviations (Table S7). All models identified marker-trait associations on chromosome arm 5AS with concurrent effects for various quality traits across environments. Several chromosome regions with SNPs having concurrent effects on quality, drought, heat, yield and yield-related traits were detected using the six models across the four environments (Tables S6-S11).

Genomic prediction
Genomic breeding values of the 192 spring bread wheat genotypes were determined for the 16 traits using the RR-BLUP model. Prediction accuracy was measured through the cross-validation method using 500 randomly generated datasets from 80 and 20% training and prediction sets, respectively. A wide range of prediction accuracies (PA) was found for these traits at each location ( Fig. S5A-S5D). The highest PA for grain yield was 0.62 from the 2016-2017 trial at Marchouch (Fig. 5C). High prediction accuracies across locations were obtained for peduncle length with 0.72, 0.69 and 0.64 from the trials at Wad Medani, Dongula and Sids, respectively. The prediction accuracy for Zeleny test was relatively high in all trials with 0.63 at Wad Medani and Marchouch and 0.61 at Dongula.

Phenotypic variation in agronomic, physiologic and quality traits
The aims of the present study were to investigate the phenotypic variation in 16 agronomic, physiologic and quality traits, locate the causative genes, and evaluate the genomic prediction accuracy of these traits in a panel of wheat lines grown in four very different environments. Highly significant variation among genotypes and the high heritability for most traits confirmed the suitability of GWAS analysis. Grain yield under drought/heat stress conditions is highly variable and the genetic mechanism is quite different from that under well-watered conditions (Pinto et al. 2010). Hence, adaptive morphological and physiological traits are widely targeted as coping mechanisms to give higher yield under these stress conditions (Reynolds et al. 2007). Sids station in Egypt is a well-known site with favorable conditions for wheat Fig. 4 Co-localised QTL detected on the long arm of chromosome 5A for yield and yield-related traits at Wad Medani in 2016. See Table 2 and footnote of  Table 1 for abbreviations production and varieties produced 3.6, 4.4 and 2.0 times more grain than at heat-stressed Wad Medani and Dongula and drought-stressed Marchouch, respectively, in 2015-2016. Grain yield was positively correlated with DH, DM, BM, PH, GRPS and GFP, in agreement with previous reports (Dodig et al. 2012;Sukumaran et al. 2015) and with WX (Guo et al. 2016), but was negatively correlated with canopy temperature at flowering across environments (Saint Pierre et al. 2010). Furthermore, grain yield was negatively correlated with grain quality parameters across environments, a well-known relationship (Kaya and Akcura 2014;Thorwarth et al. 2018).

GWAS analysis
The six implemented GWAS models frequently identified significant MTAs in similar chromosome regions for a particular trait. However, some models or only a single model also detected QTLs for different traits. For instance, all six models detected a major QTL for grain yield on chromosome arm 5AL at the heat stressed sites, whereas only the three single-locus models identified a QTL for grain gluten content on chromosome arm 2DS and the SUPER model detected highly significant multi-MTAs (195 MTAs within 4 cM) on chromosome arm 5BL for plant height. This result clearly illustrates one of the advantages of multimodel GWAS analysis for thorough QTL analysis. Most identified MTAs were environment-specific, signifying the presence of QTL x E interaction. However, several across-environment MTAs were detected at the Wad Medani and Dongula sites due to similar environmental conditions. These results provide the basis for further marker-assisted selection research to enhance yield, yield stability and end-use quality of spring wheat in heat stressed environments.

GWAS and GP for yield-related traits
Multi-model highly significant MTAs were detected on chromosome arm 5AL for grain yield and yieldrelated traits in the heat stressed environments, in agreement with previous studies (Reif et al. 2011b;Sukumaran et al. 2015;Tadesse et al. 2019;Li et al. 2019). These MTAs were located close to the vernalization locus Vrn-A1 (Yan et al. 2004). However, analysis using days-to-heading as a covariate has indicated that these markers are independent of flowering dates of the varieties (Tadesse et al. 2019 Gahlaut et al. (2019) and another on chromosome 6A (wsnp_Ex_c29368_38408543, 77 cM) close to BobWhite_c30930_192 (71 cM) discovered by Sukumaran et al. (2015). Regions in chromosomes 5B, 7D and 2A were revealed as multimodel significant MTAs for GFP across the heat and drought environments and some of the regions were detected and validated previously (Paliwal et al. 2012;Tiwari et al. 2013). TG0019 in an unassigned chromosome that revealed a MTA in all six GWAS models at Wad Medani could mark a newly discovered QTL for GFP in heat stressed environments.
Maintaining higher GRPS in drought and/or heat stress environments is an important breeding target. This trait positively correlates with grain yield under stress conditions (Shokat et al. 2020). The multimodel MTA for GRPS on the long arm of chromosome 3A detected at Sids was reported previously (Li et al. 2019;Shokat et al. 2020). The GLM, MLM and MLMM models detected significant multi-MTAs for GRPS on chromosome arm 2BL as reported in previous studies (Liu et al. 2018;Tadesse et al. 2019). Li et al. (2019) identified an environmentally stable MTA for kernels per spike on chromosome 1AL, close to SNP marker Ku_c10813_1122, revealed a multi-model MTA for GRPS at Dongula. The multi-model significant MTA for PH on chromosome arm 1AL was previously reported by Sukumaran et al. (2015). KukRi_c205_223 at 77 cM on chromosome arm 2DS, according to the wheat consensus map (Wang et al. 2014b), exhibited a multimodel significant MTA for PH at Dongula. Likewise, RAC875_c48703_148 at 9.2 cM on the same chromosome was reported by Sun et al. (2017). This QTL could be linked to the Rht8 gene (Korzun et al. 1998). Ma et al. (2018) also detected a PH QTL with pleiotropic effects on other yield-related traits in the same chromosome. GWAS model SUPER detected significant multi-MTAs for PH on chromosome arm 5BL, including marker wsnp_Ku_R-ep_c70220_69775367 (69 cM) that unveiled a significant MTA across all six GWAS models. Wang et al. (2017) detected Excalibur_c1925_2569 (132 cM) with significant MTA for PH on the same chromosome. Ku_c15539_433 (157 cM) on chromosome arm 7BL had a multi-model significant MTA for PH at Marchouch. This marker was quite distant from markers wsnp_Ex_rep_c68762_67626384 and Excal-ibur_c50612_409 (28 cM) detected by Wang et al. (2017) with a multi-environment significant MTA. SNP markers AX-111600193 (Li et al. 2020) and Kukri_c77040_87 (Ain et al. 2015) on chromosome arm 4AS with effects on PH were 50 cM from ExcalibuR_c13276_1322 with a significant MTA for PH at Sids. Considering the distance between the markers, ExcalibuR_c13276_1322 could mark a novel QTL with a profound effect on yield under favorable conditions.
Peduncle plays a considerable role in wheat productivity, such as facilitating transporting assimilates to the filling grain, lessening the risk of leafborne pathogen infections, and maintaining significantly higher water potential under drought stress or high temperature conditions (Kong et al. 2010). Multimodel significant MTAs were detected for PL on chromosome arms 3AL (two), 4AS, 5BL and 7AL and some of them were reported previously. For example, SNP marker BS00039498_51 detected by Sukumaran et al. (2015) at 122 cM was adjacent to the currently detected marker wsnp_Ex_Rep_c66357_64540369 (124 cM) on chromosome arm 3AL with significant MTA for PL at Sids. Since they were not reported in any previous study, the multi-model significant MTAs detected on chromosome arms 4AS, 5BL and 7AL might be new.
This study estimated the genomic breeding values of genotypes for yield and yield-related traits and a wide range of genomic prediction accuracies was determined. The prediction accuracies for yield and yield-related traits were generally promising, ranging from 0.21 for GY at Sids to 0.72 for PL at Wad Medani. Lozada et al. (2019) reported a genomic prediction accuracy of up to 0.56 for yield and yieldrelated traits in winter wheat. Juliana et al. (2020) reported genomic prediction accuracies in the range 0.34 -0.59 for grain yield in the CIMMYT wheat breeding program after analysing a large number of trials conducted over several years and locations.

GWAS and GP for physiological traits
Multi-model significant MTAs for canopy temperature were discovered on multiple chromosomes. Except BLINK and FarmCPU, four models detected wsnp_CAP12_c1101_569783 with a significant MTA with CTD1 on chromosome arm 4BS. Tahmasebi et al. (2017) also found two QTLs on chromosome 4B with high phenotypic effects for canopy temperature at the vegetative stage. However, due to the difference in fingerprinting platforms, we could not determine whether those QTL overlapped with the currently detected QTL. Our QTL was possibly the semi-dwarf gene Rht-B1 (Pearce et al. 2011) as Rebetzke et al. (2013 reported that the reduced height was associated with higher canopy temperature. Pinto et al. (2010) detected environmentally stable QTL on chromosome 7A in agreement with the present marker TduRuM_-contig12722_779 associated with CTD1 at Wad Medani. In agreement with this study, Pinto et al. (2010) detected QTL for CTD2 on chromosomes 1B, 3B and 5A. These could be targets for marker-assisted selection aimed at wheat improvement in drought/heat stress environments. SNP markers wsnp_CAP12_c1101_569783 on chromosome arm 4BL and IACX5821 (unassigned chromosome location) exhibited multi-model significant MTAs for CTD1 and CTD2, respectively, could represent new QTL for these traits.
Waxiness or glaucousness is a polygenic trait involved in multiple biosynthetic pathways and its expression is affected by environment (Broun et al. 2004). RFL_Contig4849_702 on the chromosome arm 2BS exhibited an environment-specific MTA across the six models at Wad Medani (Fig. 2B). The six models also detected ExcalibuR_Rep_c109101_94 on chromosome arm 2DS having a significant association with waxiness at Sids. Two wax inhibitors and two wax production genes underlining variation in glaucousness; Iw1 on chromosome arm 2BS (Wu et al. 2013;Adamski et al. 2013), Iw2 on chromosome arm 2DS (Nelson et al. 1995;Wu et al. 2013), W1 on chromosome arm 2BS (Lu et al. 2015), and W2 on chromosome arm 2DS (Tsunewaki and Ebana 1999), respectively, were identified previously. An additional gene loci, W3 was also detected previously on the chromosome arm 2BS and its mutation led to a reduction of 99% in b-diketones accounting for 63.3% of the total wax load of the wild-type (Zhang et al. 2015). Genomic regions in chromosome arms 5BL, 6BL and 7BL possibly represent new QTL for waxiness.
Leaf rolling was identified as a physiological indicator of drought-and heat-avoidance (O'Toole et al. 1979). The genetic bases of leaf rolling in drought/heat stress environments has not been extensively studied in wheat, but limited studies have identified QTLs with variable effects (Peleg et al. 2009;Tahmasebi et al. 2017). We detected multimodel MTAs for leaf rolling on chromosomes 1B, 4B, 5A, 6A and 6B and some of these could have a role in marker-assisted selection, fine mapping and/or gene cloning aimed at utilizing this trait to improve wheat production in environments exposed to heat and drought stress.
The genomic prediction accuracies of genomic breeding values of genotypes for the four physiological traits were highly variable, ranging from 0.16 for leaf rolling to 0.57 for canopy temperature at grain filling at Wad Medani. A higher average PA (0.75) obtained by Würschum et al. (2020) for epicuticular waxiness in winter wheat was much higher than the current study (0.47), but the number of worldwide winter wheat genotypes used in the former was fivefold more than the present study.

Marker-trait associations and GP for qualityrelated traits
Chromosome regions revealing multi-model significant MTAs for quality traits PR, GL, W and Z were detected in this study. Both PR and GL exhibited environment-specific QTL indicating significant effects of genotype x environment (G x E) interaction. Homoeologous group 6 chromosomes were involved in multi-model significant MTAs for all quality traits. These MTAs were likely related to storage proteins gliadins and low molecular weight (LMW) glutenins Nelson et al. 2006). The multimodel significant MTAs for PR, GL and W detected on the short arms of homoeologous group 6 chromosomes (Table 4) were likely determined by variation at the Gli-2 gliadin loci (Payne 1987). Homoeologous group 2 chromosomes also contributed multi-model significant MTAs for all quality traits; these were likely associated with the photoperiod-sensitivity Ppd1 genes (Scarth and Law 1984;Royo et al. 2016). The alleles coding for insensitivity to short day length enhance grain yield in low latitude and short-season conditions by allowing earlier heading. This leads to earlier grain filling before more extreme heat and/or drought stress periods that come with the onset of summer (Worland et al. 1988;Nelson et al. 2006). Although a thorough analysis of the relationship between Ppd1 genes and quality parameters of wheat has not been documented, Nelson et al. (2006) reported a strong correlation between early heading and protein concentration among 114 recombinant inbreed lines generated from a synthetic-hexaploid (W7985) 9 Opata 85 cross. Marker RAC875_c24504_119 on chromosome arm 3BL was another source of multi-model highly significant MTAs with concurrent effects on alveograph strength and Zeleny sedimentation across three environments. Carter et al. (2012) extensively studied chromosome 3B as a source of genes affecting several milling and baking quality traits in soft white spring wheat but did not report this particular QTL. Genomic PA for quality traits was in a range of 0.38-0.63 recorded for PR at Wad Medani and Zeleny at the Wad Medani and Marchouch sites, respectively. Lado et al. (2018) reported PA for these baking quality traits in a similar range to the current study.

Conclusions
A comprehensive GWAS and GP analysis was conducted for 16 agronomic, physiologic and quality related traits using 192 spring bread wheat genotypes grown in four environments for two years. Several multi-model significant MTAs were detected. Most MTAs were environmentally specific indicating the presence of high QTL x E interaction. Chromosome arm 5AL proved to be the source of QTL for yield and yield-related traits at the heat stress sites. The long arm of chromosome 3B was with multi-model significant MTA for TKW detected at Sids. A stable MTA for PH was detected on chromosome arm 5BL across environments and models. Chromosome arm 4BS had a QTL for CDT1 and all models detected significant MTAs for WX on chromosome arm 2DS. This last chromosome region and homoeologous regions were previously associated with both wax production and wax inhibition. QTL for quality traits were associated with loci controlling storage proteins gliadins and LMW glutenins. All six GWAS models identified highly significant MTAs for Zeleny sedimentation and alveograph strength on chromosome arm 3BL. The identified QTL can be targeted for marker-assisted selection for developing superior haplotypes with positive effects for targeted traits or for further studies on fine mapping and cloning of the underlying causal genes. Genomic prediction was applied to various traits to estimate the genomic breeding values of individual genotypes in specific environments. The genomic prediction accuracies were generally encouraging for most traits, indicating that implementing the wheat-breeding program at ICARDA could be worthwhile.
Acknowledgements The Grains Research and Development Corporation (GRDC) of Australia is acknowledged for supporting this research.
Author contributions WT, SS, ST, AH, A-A, AAA conceived and designed the study. ST, WT, AH, SS collected the phenotypic data. AA conducted data analysis and prepared the manuscript. WT, SS, AAA, A-A reviewed and edited the manuscript. All authors approved the final version of this manuscript.
Funding Open access funding provided by Swedish University of Agricultural Sciences. Grains Research Development Corporation (GRDC), Australia (grant number-ICA00013), supported this study.

Declarations
Conflicts of interest The authors declare that no competing financial interests or personal relationships influenced the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.