GWAS in a Collection of Bulgarian Old and Modern Bread Wheat Accessions Uncovers Novel Genomic Loci for Grain Protein Content and Thousand Kernel Weight

1. Introduction

Wheat is a unique cereal crop due to the baking qualities of its flour and occupies a central strategic role in broad-scale food security. Increasing yield and protein content in wheat grain is essential for safeguarding the human rights to sufficient and nutritious food. The quality and nutritional value of food products made from wheat flour largely depend upon the type and concentration of grain proteins [1,2]. As a consequence, total grain protein content (GPC) is one of the main determinants of both baking quality and the international market price of wheat [3]. The range of protein in wheat grains is usually between 8 and 20%, accounting for less than 8-15% of the grain dry weight. Thousand kernel weight (TKW) is a measure of grain weight and, along with spike number per unit area and grain number per spike, is a main component of wheat total yield. Grain weight and, respectively, grain yield, is formed mainly by the starch accumulation in the developing grain and, therefore, any starch gain in the endosperm if not accompanied by an adequate rise in nitrogen (N)/protein accumulation has a dilution effect that influences the concentration of grain protein and micronutrients [4]. Due to this important interaction, breeders frequently encounter conundrums when aiming for genotypes that combine high yield and high protein content in the grain [5,6,7,8].

Comparative studies on changes that have occurred in wheat varieties released or introduced since the middle of the 19th century showed increased grain yield and decreased protein over time [4,9]. This suggests that old germplasm (landraces and traditional varieties) can be screened for new genetic diversity and targeted for the construction of new varieties. Old germplasm grown in Bulgaria has some desirable traits, such as consistent yield, resistance to drought, high protein content, or good quality for making bread [10]. The semi-dwarf high-yielding varieties that emerged in the 1970s replaced the landraces and the tall varieties that were previously grown [11]. Most of this germplasm is now extinct, but some seed samples are preserved and reproduced in the major European seed gene banks [12].

Both GPC and TKW are quantitatively inherited and are controlled by multiple genes or quantitative trait loci (QTL) [13,14]. The impact of the environmental conditions on gene expression and the genotype by environment interactions further complicates the precise evaluation of these traits [15].

Genome-wide association studies (GWAS) detect genetic effects based on linkage disequilibrium (LD) in natural germplasm collections and have become effective tools for modern plant breeding [16]. Following the rapid development of DNA marker technologies, and in particular the advent of single nucleotide polymorphism (SNP) chips, association mapping has been used progressively to establish a strong connection between a genome-wide SNP and a trait of interest. Identifying marker-trait associations (MTAs) or quantitative trait nucleotides (QTNs) can make it easier for breeders to choose the best genotypes, reduce the breeding cycle and achieve higher genetic gains. A number of recent GWAS studies reported genomic regions associated with GPC and yield components [14,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31].

The association panels employed in these studies varied in terms of diversity level, genetic relatedness, and the nature of accessions. For instance, Kartseva et al. [27] used a diverse population of 255 accessions from 27 countries on 5 continents, and revealed novel stable genomic regions harboring GPC-associated markers on chromosomes 3A and 3B. QTL hotspots containing 165 significant MTAs for quality and agronomic traits were mapped on almost all chromosomes in an association panel of 170 diverse landraces from the Mediterranean region [31]. Using a set of 93 spring common wheat varieties and breeding lines adapted for cultivation in the Siberian region of the Russian Federation, another study reported eleven genomic regions associated with GPC, of which nine were physically mapped on chromosome 6A harboring the NAM-A1 gene, homoeologous to the Gpc-B1 (NAM-B1) gene [23]. Another recent investigation by Tyrka et al. [30] screened 168 Polish breeding lines of common winter wheat for a number of agronomic traits, and identified trait-associated markers for heading time, lodging resistance, plant height, and TKW. Therefore, it is essential to evaluate new large mapping populations with different underlying sub-structure, extents of genetic relatedness among individuals, and LD decays, to identify novel QTNs or QTL for yield and quality traits. Additionally, these traits are influenced by environmental factors, therefore, detecting MTAs in multi-environment or multi-year studies is crucial for their application in marker-assisted breeding.

For this study, we assembled an association panel of advanced and historic varieties, and local accessions from Bulgaria. Our fundamental focus was to explore the natural genetic variation for protein content in wheat grain and thousand kernel weight, and to reveal the associated genetic determinants. We found a considerable amount of phenotypic variation for GPC and TKW with no apparent differences between modern high-productive varieties and old germplasm, nor did the traits significantly correlate with one another. Based on GWAS findings, we identified promising genomic regions for wheat improvement and uncovered possible candidate genes. These results provide information about the genetic resources available to breeders to improve grain yield and nutritional properties of wheat products, as well as an opportunity to develop closely associated markers to aid molecular breeding of new varieties.

2. Results

2.1. Phenotypic Variation

Descriptive statistics, frequency distribution, and boxplots showed wide phenotypic variation for both GPC and TKW under three growing seasons (harvests 2014, 2017, 2021) within the set of 179 wheat accessions with total average over the crop seasons of 13.4 % ± 1.34 for GPC, and 45.2 g ± 4.28 for TKW (Table 1, Tables S1 and S2, Figure 1). The range of coefficients of variation (CV, %) across the individual years was similar for the two traits. Estimates of high broad-sense heritability (h²) for both traits showed moderate to high values in the individual growing seasons (from 0.64 to 0.78), and high values over the environments (0.82 for GPC and 0.81 for TKW). To exclude the impact of the growing season, we computed BLUE values for each accession, treating genotype as fixed and growing season as random effects. BLUEs varied across the years from 11.6 to 14.7 %, on average 13.4 % for GPC, and from 39.5 to 51.8 g, on average 45.4 g for TKW (Table 1, Tables S1 and S2, Figure 1).

The ANOVA results explained the presence of broad phenotypic variation among genotypes for GPC and TKW, revealing highly significant effects of genotype, environment (growing season), and their interaction (Table 2). For both traits, no significant differences were noted only between 2014 and 2017 (Table S3).

To assess the trait consistency across the environments and to explain the relationships between GPC and TKW, Pearson’s correlation coefficient approach was used. Low to high positive Pearson’s correlation coefficients (r) over the years were computed for GPC (spanning from 0.50 to 0.92) and TKW (ranging from 0.30 to 0.79) (Table 3). In general, consistency was noted across the growing seasons, with one exception – GPC in 2017 was not correlated with that in 2021. From the perspective of the data desirability for GWAS, correlation analysis was performed also with the BLUE mean values. The Pearson’s correlation coefficients (r) computed based on the BLUEs were positively significant (p<0.05) for both GPC and TKW. The two grain characteristics were not correlated across the growing seasons, as well as based on the BLUEs (Table 3).

According to [32], the studied population has a distinct structure, encompassing three sub-populations (SPs), of which SP1 (109 accessions) and SP3 (16 accessions) included predominantly modern releases, while the old germplasm formed a separate cluster, SP2 (49 accessions), and 5 genotypes were considered as admixed. Taking into account the distinct population structure, we tested the hypothesis of whether the old accessions differ significantly and consistently from the modern releases. Therefore, we compared the three SPs across the growing seasons. While the two sub-groups containing modern varieties did not differ significantly across the years and with the BLUEs concerning both GPC and TKW, the sub-group of old accessions (SP2) displayed variability across the environments. Sub-population SP2 had significantly lower GPC compared to the SP3-varieties in 2021, but did not differ from both groups of modern releases in 2014, 2017 and with the BLUEs (Figure 2a). The group of old germplasm showed a higher TKW mean value in 2017 but a lower value in 2021 when compared to SP1, whereas the TKW mean BLUE of SP2 was significantly lower than that of SP3 (Figure 2b).

To distinguish genotypes that could be used for improving grain protein content in breeding programs, the accessions were classified according to [2] into five groups with protein ≥13% (Group 1), ≥12% (Group 2), ≥11% (Group 3), and >10% (Group 4), <10 (Group 5). The prevailing part of the accessions (90 and 94% of SP1 and SP3 modern releases, 84% of SP2 old accessions, or 89% of the entire population) fall into Groups 1 and 2 (Figure 3).

In order to evaluate the phenotypic variability of the population throughout the crop seasons, we initially determined the deviations from the respective yearly average GPC and TKW. Subsequently, these discrepancies were averaged for each accession throughout the three seasons. The influence of genotypes fluctuated per the environmental gradient, with the variance extending from 0.012 (2014) to 0.025 (2017) for GPC (Table S4), and from 0.011 (2017) to 0.022 (2014) for TKW (Table S5). The variance within the entire population over the years spanned from 0.0 to 0.108 for GPC, and from 0.0 to 0.101 for TKW.

2.2. Linkage Disequilibrium (LD) Estimation, Significant Quantitative Trait Nucleotides (QTNs) and Quantitative Trait Loci (QTL)

From the perspective of GWAS, LD was considered. This is essential to define the interval of highly associated SNPs and to identify the most significant loci [16]. Here, we calculated LD decay using the whole association panel irrespective of the genotype status (old vs modern). The LD decay values varied from 1.5 to 3.0 Mbp on the individual chromosomes, with the highest value of LD decay in the D-genome (2.54 Mbp) and homoeologous group 3 (2.27 Mbp), on average 1.98 Mbp (Table S6).

For GPC, no significant QTNs were detected across the growing seasons (data not shown). However, based on BLUEs, a total of 293 significant QTNs were identified at -log₁₀ (p) > 5.53, distributed on 18 chromosomes with the highest number (84) on chromosome 3A (Figure 4a).

We combined the adjacent significantly associated markers into 25 LD blocks as trait-associated QTL. These loci were positioned on 13 chromosomes (1A, 1B, 1D, 2A, 2B, 2D, 3A, 5A, 5B, 5D, 6A, 6B, and 7A), with the highest number of QTLs (5) on chromosome 3A (Table 4 and Table S7, Figure 5). If we assume that the adjacent significant SNPs are independent and there is no interaction between them, and therefore, their effect is additive, the total effect of the detected QTL represents a sum of the individual SNP effects [33]. The results underscored the presence of 13 QTL with positive additive effects on protein content ranging from 5.28 to 49.26, while the rest loci had negative summary effects ranging from -0.91 to -40.52, thus reducing the protein content. The proportion of phenotypic variance explained by the significant SNPs within QTL (R² in %) ranged from 13 to 23% (Table 4). Locus QGpc.ippg-1A.2 with the highest number of individual SNPs with positive effect (11 out of 15), 14 to 20% individual markers contribution to the variation, and a total effect of 49.26 is promising for increasing the protein content.

For the yield-related trait TKW, a total of 483 significant QTNs (-log₁₀ (p) > 5.53) were identified on all chromosomes (Figure 6a). Twenty-nine of them were detected in two growing seasons (2017 and 2021), and the rest 454 were based on BLUEs. Clustering of adjacent significant trait-associated SNPs/QTNs within LD blocks formed 35 QTL on 13 chromosomes (1A, 1B, 2B, 3A, 3B, 4B, 5A, 5B, 6A, 6B, 6D, 7A, 7B) with the highest number on chromosomes 2B (6), 3B (5) and 6A (5) (Table 5 and Table S8, Figure 5). Based on the additive assumption, most of the QTL associated with TKW displayed neutral, slightly negative, or slightly positive effects on the trait. Our findings emphasized one QTL (on 6A) with а strong negative effect, and two loci with a high positive impact (on 6D and 7A) on the trait (Table 5). Noteworthy, 14 SNPs distributed across chromosomes 1A, 1B, 2B, 6A, and 6B were deemed environmentally stable because of their strong association (-log₁₀ (p) > 5.53) with TKW BLUE values and with phenotypic data gathered throughout at least one growing season (Table 5, Figure 5). Eight of them (on chromosomes 1A, 2B, 6A and 6B) were not included in LD blocks (Table 5). Both types of TKW-associated loci (QTL and independent stable QTNs) explained between 8 and 22% of the trait variance. The highest total additive effect (26.55) and the highest additive positive effect (50.01) on the trait was calculated for locus QTkw.ippg-6D having 28 positive SNPs out of 40, with the individual amount of variation 16 to 21% indicating that this region may be a relatively major QTL.

Collectively for the two traits, the majority of loci were mapped on A (26) and B (30) sub-genomes, while the D sub-genome carried a minor portion of QTL (4) (Table 4 and Table 5). The Q-Q plots supported the detected strong association signals for GPC and TKW (Figure 4b and Figure 6b).

Two genomic regions were shared between GPC and TKW (Figure 5). A short overlapping region on chromosome 2B (length of 55,849 bp) includes the interval from 646.9 to 647.0 Mbp and represent an intersection of QGpc.ippg-2B.1 and QTkw.ippg-2B.4, located in the support region of the flanking SNPs. The second one is a larger region on chromosome 6A extending from 607,9 to 613.0 Mbp (length of 5,087,923 bp), and includes the entire loci QGpc.ippg-6A.2 and QTkw.ippg-6A.5 (Table 4 and Table 5). The overlapping QTL in the chromosome 2B common region showed negative effects on both traits, while the loci in the shared region on chromosome 6A had high negative effect on protein content and slightly positive effect on TKW.

2.3. Potential Candidate Genes

To hypothesize potential candidate genes underlying the analyzed traits, we explored the 25 and 35 QTL detected for GPC and TKW, respectively, as well as the eight environmentally stable QTNs for TKW that were outside LD blocks. The search for candidate genes was performed based on the physical position of the flanking QTNs on T. aestivum L. cv. Chinese Spring reference genome extended by LD interval estimated for each chromosome based on LD decay. Within the defined QTL, we found a total of 1,460 high-confidence genes associated with GPC (Table 4 and Table S9). For TKW, a total of 2,477 high-confidence annotated genes within the 35 QTL, and 215 genes in the vicinity of the stable independent QTNs were evaluated as potential candidate genes (Table 5, Table S10).

In Table 6 and Table 7, we assembled shorter lists of selected potential candidate genes based on gene annotations and a literature survey. Relevant candidate genes associated with GPC and TKW encoded proteins involved in senescence-associated proteolysis, in nutrient remobilization and allocation from leaves to grains; storage proteins; a variety of transporter proteins; proteins involved in starch metabolism; the components of the ubiquitin-protease proteolytic system, and a number of transcription factors, such as NAC-domain proteins. Certain genes associated with early seed germination, such as subtilisin-like serine protease, alpha-amylase, thioredoxin and germin-like proteins, were shown to be specific to GPC. Regarding TKW, genes were found that code for proteins controlling grain size and weight, including expansins, cell wall invertase, and serine carboxypeptidase.

3. Discussion

3.1. Phenotypic Variation

To meet the growing demand for sufficient amounts and quality of food and overcome the challenges posed by environmental changes, it is essential to effectively utilize the available genetic resources of bread wheat [34]. Here, we used a collection of old and modern bread wheat accessions, and explored the genetic variation of TKW and GPC, the two most important characters determining yield and end-use quality, and eventually, the economic value of bread wheat. The study revealed large phenotypic variability for the target traits with highly significant contributions of genotype, environment, and genotype by environment interaction. The observed moderate to high values of broad-sense heritability and the consistency of trait records evidenced by the significant correlations between the crop seasons suggest that a considerable part of the variation is due to inherent genetic differences among the accessions. These results agree with similar findings in bread and durum wheat for GPC [19,25,27,35], and TKW [13,36]. Given that the panel contains historical and contemporary varieties released or collected between 1925 and 2010, it is possible that the recorded genetic variance—especially for TKW—is inflated. This is because throughout the course of nearly nine decades, genetic factors have contributed significantly to increases in yield and yield-related attributes. The landmark of these genetic gains for the Bulgarian wheat collection is the introduction of semi-dwarfing genes in the 1970s [11]. The pleiotropic effects that these genes exert on plant responsiveness to N applications, on photosynthetic rates, and on the accumulation of carbohydrates in the grain ultimately have improved grain yield.

This study shows that a large proportion (89%) of the accessions have protein above 12%, and are suitable for preparing leavened breads [2]. The average grain protein values (13.4%) are similar to or higher than the reported information on protein contents in other bread wheat collections [19,27,28,31,37]. Additional research revealed that some accession sets had greater protein values (up to 20%), demonstrating the impact of genotype by environment interactions [22,23,38].

Wheat yield and quality are affected by climatic factors (temperature, precipitation, drought type) [39]. In contrast to the abundant rainfall that was recorded from April to July in 2014, stretching from heading time to maturity, the amount of rain that fell in May and June 2017 was consistently less than the average for the climate (Figure S1), which suggests a sustained moderate drought during anthesis and grain filling. Prolonged water insufficiency early in grain development reduces the number of amyloplasts and endosperm cells, which lowers the capacity of starch accumulation, and ultimately lowers grain weight [40,41]. In 2021, a modest drought during anthesis and a more severe final drought were seen (Figure S1). Minor water stress following anthesis can speed up grain filling in wheat by boosting the activity of key catalytic enzymes that convert sucrose to starch, and by remobilizing nonstructural carbohydrates from the vegetative tissues to the grain [41]. The effects of these patterns of drought (timing and severity) explain well the significantly higher mean TKW value obtained in 2021 than that in 2017. The considerably lower mean TKW recorded in 2014 is consistent with research showing that prolonged soil wetness after anthesis restricts the amount of assimilates available to growing grains, hence diminishing the development of grain yield [42]. The observed variance in the protein content in the collection under study may potentially be explained by the precipitation oscillations among the three crop seasons. In our study, we found a significant difference (p < 0.001, Table 1) in the average protein values for 2014 (13.2%) and 2021 (14.2%), when prolonged waterlogging or terminal drought periods, respectively, were suggested during late vegetation. A decline in GPC has been associated with waterlogging [43]. Conversely, during a drought an increase in the seed proteins has been shown [44], possibly connected to modified carbon (C) partitioning and, hence, to a shift in the C/N ratio, that favors greater N-assimilation [45].

Significant genotype by environment interaction was seen (Table 2), and the genotypic effects varied along with the environmental gradient (Tables S4 and S5). It is interesting to note, that throughout the growing seasons varieties Mustang, Bozhana and Levent displayed consistently high positive deviation from the average GPC and high values of TKW. These accessions are potential sources for concurrent improvement of the two traits. Ancestral history indicates that the Mustang variety is descended from the old accession Yubilejna-3, which also exhibits consistently high protein levels in the grain. Additionally, its pedigree involves a hybrid derived from Agropyron sp., a wheatgrass that has been shown to contribute to seed storage proteins in wheat-Agropyron introgressions [46,47].

The lack of association between TKW and GPC in our data using Pearson’s approach shows that there is little to no decrease in seed weight in the presence of the identified loci for GPC. Similarly, no significant correlations were found between GPC and TKW across environments in studies of wheat lines derived from wild emmer [48,49]. This observation is consistent with findings by Oury and Godin [6] that genotype by environment interactions for grain yield and GPC may obscure the strong genetic background of the yield-protein interrelationship. Thus, Lindeque et al. [50] detected limited significant correlations between grain yield and protein content in a study of wheat accessions of various yielding capacities grown in environments of different precipitation trends. This result holds promise for the simultaneous genetic enhancement of the two traits.

Historical germplasm, such as landraces and traditional varieties, are an important source for bringing new genes into contemporary crops [31,51,52]. Therefore, it is strategically important to characterize these genetic resources in order to properly utilize them in pre-breeding. Our comparative analysis of variation for GPC showed that although the old germplasm (SP2) was characterized in general by high values of grain protein concentrations (Figure 3, Table S1), these values were lower than those of modern releases or did not differ notably from them (Figure 2). Moreover, the percentage of high-protein genotypes (≥ 12%) was higher in the two modern sub-populations, SP1 (90%) and SP3 (94%), compared to the old accessions of SP3 (84%). The idea that the older germplasm had superior yield-related metrics than the more recent varieties was also not supported by the observed variability in mean TKW values of SP2-accessions across the growing seasons. One likely reason is that a lot of contemporary Bulgarian varieties are either descended from old accessions, or have highly productive and high-protein Russian and Serbian ancestors in their pedigrees as Table S1 and the research by Kartseva et al. [27] show. These results suggest that wheat breeders have managed to construct improved varieties in terms of the studied traits.

3.2. Genomic Regions Associated with Grain Protein Content (GPC) and Thousand Kernel Weight (TKW)

Significant verified SNPs or genomic regions linked to GPC and TKW have been detected on all wheat chromosomes, according to recent research [14,53,54,55]. We compared the strongly associated SNPs found in this study to previously published loci or markers using the IWGSC RefSeq v.1.0 map as a reference. Our analysis identified four loci linked to GPC (QGpc.ippg-1A.2, -1A.3, -3A.2, and -3A.3) that include known QTNs [27,37,56,57], while QGpc.ippg-1A.2 and additional three loci were situated in close proximity (0.3 to 6 Mb distance) to other previously published genomic regions [19,28,57,58] (Table 4). Five SNP clusters and seven stable SNPs associated with TKW overlapped with previously reported QTL [21,53,58,59] (Table 5). Meanwhile, 14 TKW-linked SNP clusters were close (0.5 to 10 Mb distance) to reported QTL in both durum and bread wheat collections [30,36,53,58,59,60,61,62,63] (Table 5).

To our knowledge, the remaining QTL and/or stable QTNs associated with GPC and TKW detected in this study have not been previously linked to the studied traits suggesting that these loci are novel. Hence, we can say with a fair degree of confidence that some of the noteworthy loci associations with GPC on chromosomes 1A, 1B, 1D, 2A, 2B, 2D, 3A, 5A, 5B, 5D, 6A, 6B, and 7A found in the current research are unprecedented. Certain SNP/QTL related to TKW on chromosomes 1A, 1B, 3B, 5B, 6A, 6B, 6D, and all six QTL on 2B are also considered novel. The loci QGpc.ippg-1A.2, -2A, -3A.1, and -6B.2, and QTkw.ippg-6D with the highest additive positive effects and a high proportion of phenotypic variance explained are potentially distinct or close to known QTNs (QGpc.ippg-1A.2) trait-associated genomic regions that hold promise for use in breeding programs aimed at improving GPC and TKW, respectively.

3.3. Putative Candidate Genes Related to Grain Protein Content (GPC) and Thousand Kernel Weight (TKW)

We searched for potential candidate genes not only at suggestively novel loci, but also in genomic regions that validated previously established associations, given the wider genome coverage of the detected regions. The most well-known putative candidate genes that encode for various protein types related to GPC and TKW are summarized in Table 6 and Table 7, and the complete lists of all high-confidence genes are given in Tables S9 and S10. Grain production and nutrient content of cereal crops are governed by the buildup of nutrients in the grain near the end of the plant lifespan [64]. The quantity and quality of these nutrients is greatly affected by the remobilization of C and N from vegetative tissues to developing grain during senescence, as well as by the subsequent accumulation of resources for seed germination. All these processes are accompanied by expressional changes in a vast number of genes that ultimately impact both grain output and protein content [65]. Therefore, it is not surprising that we detected common putative candidate genes in the genomic regions that seemed to affect GPC and TKW in the studied association panel. The next section covers potential genes that may impact both traits and those that may affect only one.

3.3.1. Genes Related to Senescence-Associated Proteolysis, Nutrient Remobilization and Allocation from Source to Sink

Senescence-associated proteolysis in leaves results in hydrolysed transportable products which is essential for the source-to-sink allocation of nutrients. Within the GPC-related loci, we identified genes encoding cysteine-, aspartyl-, and serine-proteases, metalloproteases, and senescence-associated family proteins (DUF581). Several genes encoding the main constituents of the ubiquitin/26S proteasome system (UPS), particularly the proteasome subunit alpha type, 26S proteasome regulatory subunit family proteins, ubiquitin-activating enzyme E1, ubiquitin-conjugating enzyme E2, and E3 ubiquitin-protein ligase, were also found in the GPC-associated regions (Table 6). As for TKW, we found 23 putative genes for E3 ubiquitin-protein ligase (Table 7). TaGW2, a wheat gene that codes for E3 RING-type ubiquitin ligase, has two SNP polymorphisms that indicate its relationship to TKW, suggesting it as a potential candidate for grain development [66]. Within the TKW-associated loci, we retrieved six putative candidate genes translated into NRT1/PTR (Nitrate transporter 1/peptide transporter) proteins, and 13 genes annotated as high-affinity nitrate transporter proteins (Table 7). Products of the NRT1/PTR genes, together with the high-affinity transporter genes NRT2 (Nitrate transporter 2) are known for the uptake and translocation of nitrates and small peptides, thereby being involved in the N remobilization from source leaves to reservoir organs [67]. In the GPC-coding regions, we found genes encoding amino acid and peptide transporter proteins that are known to facilitate the transfer of amino acids across membranes and from source to sink. One gene, encoding glutamine synthetase, a key enzyme in the N assimilation, was detected on chromosome 6B (Table 6). The wheat glutamine synthetase genes were suggested to have functions for protein accumulation in grains [68], and for grain size traits, including TKW [69].

3.3.2. Genes Coding for Storage Proteins

The output of our search for GPC-related candidate genes includes three genes on chromosome 1A encoding for 11S globulin seed storage proteins. In wheat, the 11-12S globulins, called ‘triticins’, account for about 5% of the total seed protein [70], but are considered to be nutritionally rich due to the presence of a unique lysine-rich repeat motif [71]. Interestingly, six putative candidate genes with the annotated function of gamma gliadins and low-molecular-weight glutenins were discovered in the TKW-associated loci QTkw.ippg-1A.1 and -1A.2. Similarly, Giancaspro et al. [72] and Schierenbeck et al. [73] identified candidate genes for the key prolamins in TKW-coding regions. This finding is not surprising given that gliadins and glutenins make up the majority of endosperm reserve proteins in wheat grain. Two putative candidate genes for GPC and four genes for TKW were annotated as bifunctional inhibitors/lipid transfer protein/seed storage helical domain. Others also reported on such function of suggestive candidate genes underlying marker associations with grain architecture traits and TKW [59,73]. Although these genes might represent breeding interest for boosting yield components, it should be noted that the gene products, including lipid transfer proteins, inhibitors and seed reserve proteins, especially 2S albumins, are recognized as food allergens [74].

3.3.3. Genes Related to Sugar Transport and Starch Metabolism

Starchy endosperm makes up over 80% of cereal grains, hence increasing starch synthesis can improve TKW and yield. Starch synthesis in wheat grain requires the cooperation of numerous proteins with enzymatic and transporter functions [75,76]. In GPC-associated chromosome regions we detected putative candidate genes encoding sucrose transporter protein and PROTEIN TARGETING TO STARCH (PTST), while in TKW-related loci we found genes coding for SWEET (Sugars will eventually be exported transporters) proteins, sucrose synthase, and ADP-glucose pyrophosphorylase (AGPase) (Table 6 and Table 7). The gene TraesCS3A02G289800 translated into PTST (Table 6) has been recognized as putative candidate for grain protein in a recent study in a highly diverse population [27]. Candidate genes coding for SWEET sugar transporter [77] and AGPase [78] that may impact grain weight per spike and TKW, respectively, have also been proposed. Interestingly, we discovered a TKW candidate gene on chromosome 6B, which translated into cytokinin oxidase/dehydrogenase (CKX), an enzyme that controls cytokinin levels. An intriguing hypothesis suggests that a grain weight-associated locus (TaGW2-6A), which codes for E3 ubiquitin ligase, modulates the cytokinins’ expression levels through CKXs and the AGPase gene expression, which indirectly alters the rate of starch accumulation, size of endosperm cells, and ultimately the size and weight of grains [79].

3.3.4. Regulatory Genes

Regulatory networks during grain filling stage in cereals include differentially expressed transcription factors mainly enriched in NAC, WRKY, MYB, bZIP and others [64,65,80]. The output of our candidate gene search for GPC involved 11 genes encoding NAC-domain-containing proteins, seven genes for proteins with WRKY domain, and five genes for bZIP transcription factor family proteins (Table 6). As for TKW, we discovered 17 putative candidate genes encoding NAC-domain-containing proteins, and two genes for a SQUAMOSA promoter binding protein-like transcription factor (Table 7). NAC genes have been shown to regulate the transcriptional changes during early leaf senescence in wheat, including genes associated with the N and nutrient remobilization from source to sink [81,82], and with the grain accumulation of glutenin and starch [83,84], thereby presenting candidates for improved yield and quality breeding. Previous studies in rice [85] and wheat [72] reported on putative candidate genes for grain size and TKW, encoding for SQUAMOSA promoter binding protein-like, a star player for plant growth and development [86].

3.3.5. Genes Related to Early Seed Germination

In this work, we found prospective candidate genes related with GPC translated into subtilisin-like serine protease, alpha-amylase, thioredoxin, and germin-like protein (Table 6). Subtilisin-like proteases contribute to the degradation of reserve proteins during germination [87]. Genes for subtilisin-like proteases have been proposed as possible candidates for GPC [19] and seed longevity [88] in wheat. Alpha-amylase is one of the primary enzymes responsible for starch degradation to fuel early germinating grain [89]. Thioredoxins are ubiquitous tiny proteins (about 12 kDa) that can regulate the redox environment in the cell and play a key role in preharvest sprouting resistance and seed germination [90]. Germins and germin-like proteins are a broad and diverse class of plant proteins that were first identified in germinating wheat grains [91].

3.3.6. Genes Related to the Regulation of Grain Size and Weight

In TKW-associated regions, we detected two putative candidate genes coding for cell wall invertase (CWI), four genes coding for serine carboxypeptidase, and three expansin-coding genes. CWI is an enzyme that is essential for sink tissue formation and C partitioning, and has a strong correlation with kernel weight [92,93]. Orthologous genes encoding serine carboxypeptidase in rice (OsGS5) and wheat (TaGS5) have been associated with changes in grain size and weight [94,95]. Expansins affect the structure of the cell wall, cell growth, and grain tissue extension, and so play a crucial role in determining grain size [96].

4. Materials and Methods

4.1. Plant Material

The plant material consisted of 179 winter wheat (Triticum aestivum L.) accessions that originated from Bulgaria, of which 129 were modern varieties, and 50 were old accessions (historic varieties of tall stature and landraces). Seeds from the old accessions were procured from the seed gene banks at the Leibniz Institute for Plant Genetics and Crop Research (IPK), Gatersleben, Germany, and the Crop Research Institute, Prague, Czech Republic. According to the available information, the period of early breeding and research expeditions to aggregate seeds from landraces and historic varieties for gene bank collections expands from 1925 to 1970. Seeds from the contemporary varieties released until 2010 were made available from the two major breeding centers in Bulgaria (Dobrudzha Agricultural Institute, General Toshevo and the Institute of Plant Genetic Resources, Sadovo), as well as by breeders. Information about botanical variety, status (modern vs old), year of release and known genealogy is given in Table S1.

4.2. Phenotyping

The seed material was collected from field experiments conducted in Sofia, Bulgaria (42^o41´N, 23^o19´E) during three crop seasons (2013/14, 2016/17 and 2020/21), denoted by the year of harvest (2014, 2017 and 2021). For each growing season, the accessions were sown in a random design in double 1 m-long rows and two replications. The soil type at the experimental field is leached vertisol, pH 6.1, with 3.1 % humus, 1,420 mg total N per kg soil, of which 18 mg inorganic N. Plants received 120 kg N/ha as ammonium nitrate in two split doses, 40 kg/ha two weeks after sowing and 80 kg/ha before stem extension growth stage. For pest control, standard agricultural practices were applied. The average monthly temperature ranged between -5.8°C and 23.2°C, during the vegetation cycles (October-July). The monthly precipitation ranged from 4.7 mm to 150.5 mm, with average annual precipitation of 682 mm, which is ca. 20 % higher than the climate norm for Sofia. Information on the monthly weather statistics for the region of the experimental field is presented in Figure S1.

Plant material was hand-harvested and hand-threshed. For each genotype and crop season, the thousand kernel weight (TKW) was determined based on grain number and grain mass of 10 main spikes per replication. The grain protein content (GPC) was measured in three pooled whole-grain powder samples (1.0 g each) for each accession and year. Before analyses, seeds were dried to a constant weight and ground to a fine powder with IKA Tube Mill Control (IKA Werke GmbH & Co, Staufen, Germany). The total N concentration was determined using a UDK 159 Automatic Kjeldahl Nitrogen Protein Analyzer (Velp Scientifica, Italy), and the corresponding GPC values (in %) were automatically inferred from the N values, according to [97].

4.3. Statistical Analyses

The significant differences in GPC and TKW among accessions, growing seasons, and interaction effects between genotype and environment were tested by analysis of variance (ANOVA). The Best Linear Unbiased Estimators (BLUEs) for each accession across the growing seasons were obtained to eliminate the environmental impact by assuming the genotype as a fixed effect and the growing season as a random effect. The relationships for GPC and TKW among the growing seasons and with the BLUEs of these traits were tested by Pearson correlation coefficients (r).

Broad sense heritability h² for GPC and TKW was calculated using the formula:

$$h^2=\frac{{\sigma }_G^2}{\left({\sigma }_G^2+\frac{{\sigma }_E^2}{nE}\right)}$$

where ${\sigma }_G^2$ is the genotype variance, ${\sigma }_E^2$ is the variance of the residual, and nE is the number of environments (growing seasons). Heritability in each environment was estimated using the same formula, where nE in the denominator is the number of replications in a given environment.

All phenotypic data analyses were accomplished using STATISTICA 14 [98].

4.4. Association Mapping and Candidate Gene Search

Before performing marker-trait association analysis, the population stratification, the genetic relatedness among population entries and the LD were considered. The genotypic data for the association panel of 179 accessions were already available from the 25K Infinium iSelect array (SGS Institut Fresenius GmbH TraitGenetics Section, Gatersleben, Germany) and described in [32]. Population structure was modeled using the Bayesian clustering algorithm in STRUCTURE 2.3.4, and the membership coefficients (Q-values) were determined [32]. Kinship was estimated as a similarity matrix (K) from 17,083 SNPs in TASSEL v. 5 (https://tassel.bitbucket.io/, accessed on 10 December 2023). The values of LD decay (in Mbp) were determined for each chromosome as described in [32]. Considering the distinctive population structure and high genetic relatedness among population individuals, GWAS was undertaken with a mixed linear model (MLM), Q+K model, to control pseudo associations [99,100,101]. Hence, the filtered set of 17,083 SNPs, phenotypic data from the three crop seasons, and the calculated BLUE mean values for the traits, along with population structure (Q values) and kinship similarity matrix (K) as covariates were used for association mapping analysis in TASSEL v.5. Additive effects and R² (percent phenotypic variation) of QTNs were estimated in TASSEL v.5 and Excel (Microsoft). The percentage of phenotypic variation explained by each QTN (R²) was calculated as the difference of R² with and without the strongest associated SNP.

The threshold of statistically significant QTNs was set by the Bonferroni algorithm. Bonferroni corrected p-value was calculated with the formula: p-value =$\frac{\alpha }{n}$, here α=0.05 and n is the number of SNPs. The threshold above which the significant markers on the Manhattan plots are located was calculated by -log₁₀ (corrected p-values). According to our data, for the 17,083 SNP markers, the estimated Bonferroni threshold is 5.53. On quantile-quantile (Q-Q) plots, distributions of important p-values (expected vs. observed -log₁₀ (p-values)) for GPC and TKW were shown. To draw the Manhattan plots and Q-Q plots, R package qqman was employed [102]. Clusters of three or more SNPs were considered as suggestive QTL, when the flanking markers met the above criteria for significance, and the adjacent SNPs were located within a support interval based on the LD decay values estimated for each chromosome. Trait-associated SNPs obtained with BLUEs that were shared with at least one individual growing season, and at least one of them was significant (-log₁₀ (p)>5.53) were considered as environmentally stable ones.

Candidate gene regions were defined as LD blocks harboring significant flanking trait-associated markers (-log₁₀(p)>5.53) and their extension to the left and right within the LD-estimated interval. For significant trait-associated SNPs outside of LD blocks (referred to as independent QTNs), flanking regions within the LD interval on either side of the marker were searched for candidate genes. Within these regions, high-confidence putative candidate genes were predicted by blasting against the cv. Chinese Spring reference genome IWGSC RefSeq v1.0. [103]. Gene annotations, gene ontologies (GOs), InterPros and details for the potential candidate genes were obtained using EnsemblPlants, and Persephone web-based platforms (http://plants.ensembl.org/Triticum_aestivum/Info/Index; https://web.persephonesoft.com/?data) (accessed on 10 January 2024).

5. Conclusions

The GWAS analysis revealed a number of QTL underlying both studied traits. Additionally, 14 environmentally stable QTNs on chromosomes 1A, 1B, 2B, 6A, and 6B were found for TKW. While some of the loci confirmed already published ones, others might be new. Among these, four loci for GPC (QGpc.ippg-1A.2, -2A, 3A.1 and -6B.2) and one locus for TKW (QTkw.ippg-6D) had notably positive effect on the traits, thereby representing promising genomic regions for use in breeding projects aiming at improving grain protein and yield, respectively. Besides, to our knowledge all of the targeted regions except for QGpc.ippg-1A.2, have not been previously linked to the studied traits. Based on their gene ontology and annotation, a few of the many high-confidence genes that were retrieved from the trait-associated genomic areas were selected as intriguing candidates. Annotations to senescence-associated proteolysis, nutrient remobilization and source-to-sink allocation, synthesis of storage proteins, starch metabolism, seed germination, and regulatory transcription factors are present in these genes.