Genome-wide association study of yield components in spring wheat collection harvested under two water regimes in Northern Kazakhstan

Phenotyping of the collection under irrigated and rainfed conditions

In this study, we considered a collection of 179 spring bread wheat accessions including 92 commercial and prospective cultivars of Kazakhstan and Russia, 86 breeding lines from Alexandr Barayev Scientific-Production Center for Grain Farming (SPCGF, Shortandy, Akmola region), and a check cultivar for the Akmola region, Tselinnaya yubileinaya (TY, Table S1). Field experiments were conducted on the experimental plots of SPCGF (51°36′09″N, 71°02′24″E, 391 m above sea level) in 2018, 2019, and 2020, both under irrigated and non-irrigated (rainfed) conditions. Each accession was grown in 1 m² complete randomized block plots composed of seven rows with 50 seeds per row in two repetitions. Field management was consistent with local practices for wheat production. Irrigation (45 mm) was applied at two critical stages: tillering and booting. The raw meteorological data registered for experiments are provided in (Raw_meteo_data.xlsx).

The nine agronomic traits associated with drought resistance and grain productivity and used for phenology and phenotyping included the heading date (HD, days), the seed maturity date (SMD, days), plant height (PH, cm), peduncle length (PL, cm), number of productive spikes (NPS, pcs), spike length (SL, cm), number of kernels per spike (NKS, pcs), thousand kernel weight (TKW, g), kernels yield per m² (YM2, g/m²).

Genotyping of the collection

Genomic DNA was extracted from a single seedling of each individual accession using the cetyltrimethylammonium bromide (CTAB) method (Doyle & Doyle, 1990). The DNA concentration for each sample was adjusted to 30 ng/ μL. All samples of the 179 wheat accessions were genotyped using a 20,000 (20K) Illumina iSelect SNP array at the TraitGenetics Company (TraitGenetics GmbH, Gatersleben, Germany). A total of 8,662 polymorphic SNP markers were selected for GWAS using previously published criteria (Miyagawa et al., 2008). According to these criteria, markers with call rate ≥ 90%, Hardy–Weinberg equilibrium fit at P ≥ 0.001, a confidence score of 0.5, and minor allele frequency (MAF) ≥5% were considered to meet the requirements. Accessions with greater than 15% missing data were excluded from further analysis.

Analysis of linkage disequilibrium, kinship, population structure, and statistics

The statistics for yield trials were assessed using GraphPad Prism Version 9.0 (GraphPad Prism, 2021). GEI was analyzed using the genotype and genotype-environment (GGE) biplot method and the Finlay and Wilkinson (FW) regression analysis in GenStat software Version 19.1 (VSN International, 2020). The correlation analysis was calculated using Rstudio software (RStudio Team, 2015)).

A model-based clustering method (admixture models with correlated allele frequencies) in STRUCTURE v.2.3.4 software (Pritchard, Stephens & Donnelly, 2000) was used to study the population structure of the entire collection. Five runs were conducted for each K ranging from 2 to 10 with a 100,000 burn-in length and 100,000 Markov chain Monte Carlo (MCMC) iterations. The optimal number of clusters (K) was chosen based on the ΔK as described by (Evanno, Regnaut & Goudet, 2005). The obtained values were then transformed into a population structure (Q) matrix.

The linkage disequilibrium (LD) in the studied collection was separately calculated for each hexaploid common wheat genome (genome A, genome B, and genome D), as well as the average LD for three genomes using Java-based TASSEL v.5.2.53 software (Bradbury et al., 2007). The R statistical platform was used to build a plot between the pairwise R² and the genetic distance (LD decay plot) (RStudio Team, 2015). TASSEL was also applied to calculate a population kinship matrix (Kin) based on the scaled identity using state (IBS) method (Stevens et al., 2011).

Marker-trait association analysis

The TASSEL and mixed linear model (MLM) method with the application of K and Q matrices was used for the identification of QTLs for agronomic traits both under irrigated and rainfed conditions. The analysis was based on QTL phenotypic data for the nine traits obtained from field trials in 2018, 2019, and 2020 and their average values over three years. P < 1 × 10⁻³ was used as a significant threshold for identified MTAs. The positions and sequences of the SNP markers were obtained from the 90K Array Consensus map set of the common wheat genome (Wang et al., 2014). For confirmation of the correction due to K and Q matrices, the distribution lines in each quantile-quantile (QQ) plot were analyzed. In the case of several significant MTAs positioned closely to each other, the SNP with the lowest P-value was chosen. MapChart v.2.32 (Voorrips, 2002) was used to draw a genetic map.

Candidate gene analysis

For the search for protein-coding genes that overlapped the identified significant MTAs, the sequence for each marker was used in the BLAST tool of Ensembl Plants (Ensembl Plants, 2020) for comparison against the reference genome of T. aestivum.

Phenotypic variation and correlation analysis

The field performance of the 179 local spring wheat accessions was analyzed at the SPCGF (Northern Kazakhstan) under irrigated and rainfed conditions during three field seasons in 2018, 2019 and 2020. The two-tailed t- test suggested that average values in all nine studied traits were significantly different between the tested irrigated and rainfed conditions. The average PH values showed the largest difference (P < 0.0001) between the two tested conditions (Table 1), 73.6 cm in irrigated in contrast to 61.6 cm in rainfed conditions.

On average, YM2 declined by 5.7% under rainfed (332.3 ± 5.68 g/m²) compared to the irrigated (352.3 ± 4.30 g/m²) conditions. In total, 51 accessions exceeded the YM2 of the local standard cultivar Tselinnaya yubileinaya (TY, 374.5 g/m²) under rainfed conditions, including nine accessions that outperformed the standard also under irrigated conditions (Fig. 1A). The Finlay and Wilkinson (FW) regression analysis (Fig. 1B) suggested that the YM2 of four wheat accessions, particularly WS10, WS32, WS82, and WS85, was stable in all three tested years (2018, 2019 and 2020) out of the nine accessions highlighted in the box in Fig. 1A, showing YM2 values of 400 g/m² and higher.

The analysis of the average YM2 using a scattered GGE biplot indicated that 52.8% of the total variance was explained by Principal Coordinate 1 (PC1), and 47.2% by PC2 (Fig. 2). PC1 effectively separated accessions that showed the highest yield performance in irrigated and rainfed conditions, while PC2 split the entire collection into groups with higher and lower YM2 for both conditions. The GGE biplot graph essentially confirmed the results in Fig. 1A and identified the accessions with high average YM2 under irrigated and rainfed trials as well as lines that showed high yield under both conditions (i.e., WS93 and WKZ19).

Pearson’s correlation assessment in both conditions indicated that average YM2 was positively correlated with NPS and TKW (Fig. 3). Interestingly, the highest correlation value of YM2 under irrigated condition was with NPS (0.39), whereas under the rainfed condition, was with NKS (0.36). We found that earlier HD was advantageous for higher TKW under rainfed conditions, while it was not a significant factor for the yield under irrigated conditions. Under rainfed conditions, a higher PH value was not a contributing factor to YM2 (Fig. 3A). Interestingly, under rainfed conditions, HD influenced both SL and NKS (Fig. 3). Expectedly, PL was highly correlated with PH (P < 0.0001), but negatively associated with NPS (Fig. 3A).

Genetic map, population structure, and LD

The DNA genotyping data for studied 179 spring wheat accessions based on the use of the 20K SNP Illumina array resulted in the identification of 8,662 polymorphic SNP markers. The total map length for those 8,662 SNPs was 3407.6 cM, with an average chromosome length of 162.2 cM. The density of markers in chromosomes varied from 0.1 SNP per cM (chromosome 4D) to 1.7 SNP per cM (chromosome 6B).

Based on the results of the population analysis performed using STRUCTURE for the 179 accessions of wheat genotypes and STRUCTURE Harvester analyses, ΔK was optimal at K = 3. LD decay occurred at 18.5 cM (genome A), 13.1 cM (genome B), and 53.8 cM (genome D) in different genomic regions with a genome-wide LD decay of 5.0 cM at 0.1 R² (Fig. S1).

Marker-trait associations under irrigated and rainfed conditions

The phenotypic data for nine agronomic traits of the 179 wheat accessions harvested under rainfed and irrigated conditions were subjected to GWAS using the 8,662 polymorphic SNP markers. Out of 237 total QTLs , 50 stable QTLs were identified for irrigated and rainfed conditions in the Akmola region, Northern Kazakhstan, for HD, SMD, PL, SL, NPS, NKS, and TKW, but no QTLs were detected for PH (Table 2, Table S2, Dataset S1). The highest number of stable QTLs was localized on the chromosomes of genome B (26), followed by genomes A (16) and D (8). In general, 25 QTLs were identified in both rainfed and irrigated conditions (Table 2, Fig. S2 ).

In total, eleven stable QTLs were identified for HD (Table 2). Three of them were detected only in rainfed conditions, while eight QTLs were found from the irrigated sites. One of the QTLs that was common for both tested conditions (QHD.ta.ipbb-3B) also affected NKS (Table 2). The largest number of QTLs identified in the rainfed trials was associated with SMD (8 of 12 total QTLs). PH and PL were one of the key traits to this study, as the two tested sites highly significantly differed in these traits (P < 0.0001). However, only two stable QTLs (QPL.ta.ipbb-3B.1 and QPL.ta.ipbb-3B.2) were identified for PL, and no stable QTL was identified for PH. The small number of QTL found under both conditions were YM2 (two QTL), NKS and SL (5 QTLs), and NPS (6 QTLs). Although all seven QTLs for TKW were identified in both tested sites, only two QTLs were detected in rainfed (QTKW.ta.ipbb-7B) and irrigated (QTKW.ta.ipbb-7A) conditions (Table 2).

Among the 46 identified QTLs, 7 pleiotropic MTAs were detected under both conditions. Those pleiotropic MTAs were mapped on the 1B, 1D, 2B, 3B, 5A, 6A, and 7A chromosomes and associated with HD (wsnp_Ex_c8240_13914674 and Excalibur_c20376_615) and spike-related traits (NPS, NKS, and TKW) (Table S3, Fig. S2).

Yield assessment in the spring wheat collection under irrigated and rainfed conditions in Northern Kazakhstan

The analysis of the collection under the two tested conditions suggested substantial differences between irrigated and rainfed fields (Table 1), indicating that limitation of water supply significantly affected all nine analyzed traits in the study. Particularly, the largest difference between two sites was found for PH, which is congruent with results reported earlier (Tsenov et al., 2015; Lehari et al., 2019; Pour-Aboughadareh et al., 2020). The analysis of average YM2 revealed that under rainfed conditions, more than 51 common wheat accessions exceeded the YM2 of the local standard cultivar, TY, in Northern Kazakhstan (Fig. 1A). Particularly, nine accessions showed outstanding yield performance in comparison to the local check cultivar TY, and four of those lines displayed grain yield stability across all three years (Fig. 1B), suggesting that they might be an excellent source for local drought -tolerance -related breeding projects. The GGE biplot scatter method indicated that the two main principal coordinates explained 52.8% (PC1) and 47.2% (PC2) of total variability, supporting the assumption of the high diversity of the tested collection (Fig. 2).The correlation analysis for traits under rainfed conditions suggested a significantly negative correlation between HD and TKW, an important yield component, hinting that early HD might be favorable for increased yield under stressed conditions (Fig. 3A).

Comparative QTL identification for agronomic traits in irrigated and rainfed conditions using GWAS

The GWAS analysis of the wheat collection grown under irrigated and rainfed conditions in the Akmola region of North Kazakhstan allowed the identification of 50 stable out of 237 total QTLs that were significant for eight of the studied traits (HD, SMD, PL, NPS, SL, NKS, TKW, and YM2) under both conditions (Table 2, Fig. S2). In total, 25 common QTLs for these traits were identified in both conditions (Fig. 4). Thirteen QTLs under rainfed conditions were revealed for HD (4), SMD (8), and TKW (1) that were not detected in trials at the irrigated conditions. Twelve QTLs were identified only under irrigated conditions for the following traits: SMD (2), NPS (6), NKS (1), TKW (1), and YM2 (2).

The assessment of QTLs for phenological traits (HD and SMD) showed that only 3 out of 23 QTLs showed effects with longer days. The remaining 20 associations were found with QTLs with shorter flowering and seed maturation time (Table 2). Two QTLs for HD, QHD.ta.ipbb-2B and QHD.ta.ipbb-3B, showed pleiotropic effects for NPS, NKS, and TKW. an SNP marker in QTL for SMD, QSMD.ta.ipbb-2A, was also significant in the QTL for TKW on chromosome 2A (Table 2). Interestingly, despite the stark differences between PH under irrigated and rainfed conditions, no QTL was detected in GWAS for this trait. Evidently, the collection’s accessions are fixed for this trait, and the low variation within the studied conditions was not enough to identify any MTA.

A total of eight QTLs were identified in irrigated and rainfed trials for the spike-related traits SL and NKS. Interestingly, seven of those eight QTLs were revealed both irrigated and rainfed conditions, and three QTLs were detected only under rainfed trials (Table 2). QNKS.ta.ipbb-3B, one of the QTLs identified only in irrigated conditions, should be particularly highlighted, as it was the factor with the highest QTL effect for the trait. Still, this QTL effect may rather be a result of the pleiotropic effect of the QTL for HD, which was identified under both conditions (Table 2). TKW is known to be one of the major yield components in wheat (Quarrie et al., 2005). In this study, eight of nine QTLs for TKW were detected using trials under rainfed conditions (Table 2). Particularly, QTKW.ta.ipbb-1B and QTKW.ta.ipbb-2B showed equally high QTL effects for TKW.

A comparison of the mapped QTLs analyzed in this study with those from other previous studies indicated that ten QTLs had known associations. For TKW, two QTLs QTKW.ta.ipbb-2B.1 and QTKW.ta.ipbb-6A were genetically mapped close to genomic regions to QTLs for this trait identified by Tura et al. (2020). Another two associations (QNPS.ta.ipbb-1B and QYM2.ta.ipbb-7B) were in the same genetic positions of QTLs identified with the analyses of six traits using GWAS based on the assessment of spring wheat in Kazakhstan (Turuspekov et al., 2017). QTKW.ta.ipbb-5A was identified in the same genetic position as in the GWAS of 285 elite spring wheat lines of wheat association mapping initiative population grown in temperate irrigated environments (Sukumaran et al., 2015). The only QTL identified for TKW that seems to be novel is QTKW.ta.ipbb-7B, which was identified under rainfed conditions and has not previously been reported. With respect to the other identified MTAs of yield components, QNPS.ta.ipbb-2B, QNPS.ta.ipbb-5A, QSL.ta.ipbb-2D, and QTKW.ta.ipbb-6A were located in close proximity to QTLs for the same traits in the study of UK reference mapping population Avalon × Cadenza in the Northern, Central, and Southern regions of Kazakhstan (Amalova et al., 2021).

The assessment of the identified only in rainfed conditions suggested that 12 of 13 QTLs are associated with phenological traits HD and SMD (Fig. 4), which underlines the importance of plant growth-related traits in avoiding water deficiency stress. Particularly, QHD.ta.ipbb-6A and QHD.ta.ipbb-6D for HD and seven QTLs for SMD (from chromosomes 3B to 6A) were found to be essential for early flowering and seed maturation time under rainfed conditions (Table 2). The locations of QTLs for HD compared to known flowering genes showed that the position of the most important SNP for QHD.ta.ipbb-5A.1 completely coincides with the physical position of Vrn1 (587,4 Mb; Table S3). QHD.ta.ipbb-5A.3 for HD and QSMD.ta.ipbb-3B.3 for SMD were also previously identified by Sukumaran et al. (2015), confirming the robustness of the identified QTLs for HD in this study. Our literature survey indicated that the remaining nine HD and eleven SMD associations seem to be novel QTLs, as none of them were reported elsewhere (Table S3).

Localization of significant SNPs in the identified QTLs for the studied agronomic traits in the wheat physical map

The alignment of the most significant SNPs in the 50 identified stable MTAs with sequences in the Wheat Ensembl database (https://plants.ensembl.org/Triticum_aestivum/Info/Index) suggested that SNPs in 43 and 7 MTAs were in genic and intergenic positions, respectively (Table S3). Interestingly, two SNPs in the QTLs for HD under rainfed conditions on homeological chromosomes 6A and 6B (QHD.ta.ipbb-6A and QHD.ta.ipbb-6B) were aligned with E3 ubiquitin-protein ligase (UPL) (Table S3). Similar to the above two QTL locations for HD, the homeological region on the distal part of the short arm of chromosome 6D also carries an MTA for HD (Table 2, Figs. 4A, 4B). However, the SNP in this MTA aligned to a different gene (Table S3), most probably because of poor representation of polymorphic markers in the genome D. The UPL, along with E1 ubiquitin-activating and E2 ubiquitin-conjugating enzymes, is known to participate in the ubiquitylation of proteins (Liu et al., 2020). Ubiquitylation is essential for the regulation of various biological processes, including growth and development, response to biotic and abiotic stress, and regulation of chromatin structure (Ramadan et al., 2015; Xu et al., 2021). Additional confirmation of the relationship between ubiquitylation and HD in this study is provided by the SNP alignment of QHD.ta.ipbb-1B and QHD.ta.ipbb-3B with the Wheat Ensembl database. Particularly, QHD.ta.ipbb-1B was aligned with a putative ubiquitin carrier protein, and QHD.ta.ipbb-3B with ubiquitin core domain-containing protein (Table S3). In other sequence alignments of the identified QTLs, the SNP in the most significant QTL for SMD (QSMD.ta.ipbb-3B.3) was aligned with the unknown function protein. The position of the SNP for QTKW.ta.ipbb-1B showed high confidence alignment with the position of a gene from Tetratricopeptide-like helical domain superfamily, which enables plants to cope with adverse environmental stresses and allows them to rapidly acclimate to new conditions (Sharma & Pandey, 2016).

Hence, the identified SNP markers for discovered 50 stable QTLs of eight agronomic traits, including eleven QTLs for HD and twelve QTLs for SMD, can be recommended for further validation tests in spring wheat projects for efficient construction of new and highly competitive cultivars in arid zones.