Deciphering Genomic Regions and Putative Candidate Genes for Grain Size and Shape Traits in Durum Wheat through GWAS

1. Introduction

Wheat is one of the most important cereal crops consumed as a staple food by humans on almost every continent in the world. It was the first domesticated crop about 10,000-12,000 years ago around Fertile Crescent and played an important role in initiating the agricultural revolution [1,2,3,4]. The wild form of durum wheat, wild emmer (Triticum turgidum ssp. dicoccoides, AABB), is an allopolyploid in the tetraploid group and probably evolved from natural hybridization between two diploid species Triticum urartu (AA) and Aegilops speltoides (BB). After the domestication of wild emmer, the emmer wheat, T. turgidum ssp. dicoccum, formed and this led to the evolution of modern durum wheat, T. turgidum ssp. durum [5,6,7,8]. Durum wheat is an important cereal used for making pasta, flatbread, couscous, bulgur, etc., especially in the Mediterranean region [9]. About 36 million tons of durum wheat are produced on around 13 million hectares worldwide [10]. Turkey and Canada are the top producers of durum wheat, with 2 million hectares per year in each country [11,12]. Durum wheat comprises approximately 5% of total wheat production and is an economically important cereal due to its unique properties [13].

In general, the market price of durum wheat is 20-40% higher than that of bread wheat, sorghum, and corn [14]. Therefore, any approaches to increasing durum wheat yield, such as investigating grain size and shape traits, are important targets for breeders because they have a direct association with wheat yield and milling quality [15]. The grain size and shape are determined by the weight, area, length, width, perimeter, sphericity, and horizontal axis proportion of the grain [15,16]. Additionally, grain size affects seedling vigour [17,18], which is an important factor for improved stand establishment and yield. Therefore, revealing the genetic basis of grain size and shape may provide significant information to enhance wheat yield.

Linkage mapping or quantitative trait loci (QTLs) mapping is a practical method for dissecting the genetic mechanism of target traits, including yield and yield components. To date, several QTL mapping studies have been reported for grain size and shape-related traits in durum [19,20,21], bread [16,22,23,24,25], and einkorn wheat [26,27]. Despite the success of linkage mapping, it has some fundamental limitations. In linkage mapping, only the allelic polymorphism between parental lines of F2, RIL (Recombinant Inbred Lines), or back-cross populations can be evaluated to identify QTLs [28]. Therefore, linkage mapping has a lower power to identify QTLs with minor effects and may span 15-20 cM large genomic regions [29,30]. These large distances limit the resolution of mapping, especially in species like wheat with large and complex genomes. However, the genome-wide association study (GWAS) method is implemented in a variety of genotypes from different ancestors; hence, it has greater allelic variation and higher map resolution, due to long-term recombination events and larger population sizes [29]. GWAS uses single nucleotide polymorphism markers (SNPs) that are dispersed all over the genome and identify associations with agronomically important traits [31].

Previous GWAS studies reported several MTAs for grain size and shape-related traits in bread wheat [32,33,34,35,36,37,38,39,40,41], but there are limited reports for durum wheat [31,42] and diploid wheat species, such as einkorn (Triticum monococcum) [43], Triticum urartu [44], and Aegilops tauschii [45,46]. In durum wheat, Wang, et al. [42] reported five MTAs for grain length on chromosomes 2A, 3A, 3B, 6A, and 7A, three MTAs for grain area on chromosomes 3A and 7A (2), and three MTAs for grain width on chromosomes 3A (2) and 4A. In another durum wheat study, Alemu, et al. [31] identified five MTAs for grain length on chromosomes 2B, 4B, 5A, 6A, and 7B, whereas four MTAs were reported for grain width on chromosomes 2A (2), 5A and 7B. The grain size and shape-related traits in durum wheat still need to be investigated, and their genetic basis needs to be uncovered by newly developed genome-wide scan sequencing methods. SNP genotyping technologies, such as the DArTseq genotyping system (Diversity Array Technology) [47] provide a genome-wide detailed scan with thousands of molecular markers simultaneously. In addition to genotyping systems, phenotyping technologies also increase the accuracy of GWAS analyses. So, digital image analysis software has a higher superiority in terms of precise measurements in contrast to manual phenotyping. These software are increasingly used in plant phenotyping studies [48] and provide detailed imaging of grain characteristics, such as grain size and shape-related traits [49].

In this study, a panel of durum wheat was genotyped using the DArTseq genotyping system, which comprises highly polymorphic SNP markers with known chromosomal locations and sequences. The seed samples of these genotypes were screened by digital image analysis software to evaluate grain size and shape-related traits, and a GWAS was performed to identify marker-trait associations (MTAs) for these traits.

2. Materials and Methods

2.1. Plant materials

A durum wheat panel consisted of 146 durum wheat (Triticum turgidum L.) advanced lines were used as plant materials. The pedigree of the panel is submitted in Table S1.

Experiment designs and measurements

The panel was assessed in randomized complete block design (RCBD) with three replications over two years (2017-18 and 2018-19) in the experimental area of Department of Field Crops in Faculty of Agriculture, Cukurova University at Sarıçam/Adana, Turkey. The trial area is located at 29 m altitude and features an alluvial, medium pH, sandy-clay soil type with a deep and well-drained structure (37°00’45’’ N and 35°21’20’’ E). Each season was accepted as a different environment and named as follows: E1 (2017-2018) and E2 (2018-2019).

The panel was sown in the plots, which had two 2-m rows with a row spacing of 20 cm and a genotype spacing of 10 cm. In each row, there were 20 seeds. During the growth season, common agricultural practices as fertilization, irrigation, disease management, and pest and weed control were used. After physiological maturity, completely mature spikes (Zadoks Scale, GS93) [50] were harvested for every individual plant. Ten randomly selected spikes were manually threshed, and the seeds were bulked for each replication to calculate thousand-grain weight (TGW) and measure grain size and shape characteristics. TGW was determined by manually counting 250 grains for each replicate and then weighing them in an ultra-analytical balance to convert them to thousand-grain weight. Thirty seeds were used per replication for grain size and shape analyses. All seed samples were first photographed, and then the images were transferred to a computer environment to measure the grain size and shape traits. The following grain size and shape traits were measured using Smart Grain 1.2 software [49]; area size (AS), perimeter length (PL), grain length (L), grain width (W), length-width ratio (LWR), and circularity (CS).

2.2. Genotyping-by-sequencing (GBS) analysis

A total of 76.265 SNP markers were obtained from the DArTseq genotyping system. The durum panel was screened by the markers that were produced from the Wheat Chinese Spring IWGSC RefSeq v1.0 genome assembly. Therefore, the D genome markers were eliminated from the raw data, and then the A and B genome markers were updated in the Durum Wheat Genome (cv. Svevo) V1 (https://urgi.versailles.inra.fr/blast/). Each marker sequence query was set to “Best score” to update the marker in the Svevo genome. The chromosome and chromosomal location of each marker were reassigned according to the alignment similarity, which exceeded 95%. The remaining updated markers were then subjected to additional filtering procedures, such as the removal of markers with more than 20% missing data and fewer than 10% minor allele frequency, in order to generate acceptable polymorphic markers. After data refinement, a total of 3,251 high-quality and polymorphic SNP markers were obtained for use in GWAS analysis. Heterozygous alleles, by the way, were disregarded and marked as missing alleles for each marker.

2.3. Basic statistical analysis

The ANOVA was performed by the “metan” R package [51], which was developed for multi-environment variance analysis. The same package was used to calculate the broad-sense heritability and coefficient of variation (CV) for all traits. JASP software Version 0.11.1 was used to generate distribution plots and Pearson’s correlation coefficients between traits as well as between environments [52]. Best linear unbiased predictor (BLUP) values were calculated for each trait across two environments in JMP Genomics 9.0 software [53]. BLUP values will be considered a third environment for use in GWAS analyses.

2.4. Population structure and linkage disequilibrium (LD) analyses

The Eugene vector principal component analysis was used to describe the population structure and was plotted in the GAPIT package of R software [54] to determine the number of principal components for use in GWAS analysis. LD between pairwise comparisons of SNP markers with a sliding window size of 50 markers was estimated using the squared correlation coefficient (r²) in TASSEL 5.2.86 [55]. The LD decay plot was created by the LD results (pairwise r² values) from TASSEL 5.0 for the whole genome in R Studio 2022.07.1. LD blocks were created by HAPLOVIEW v4.1 software (https://www.broadinstitute.org/haploview/haploview) [56].

2.5. Genome-wide association analysis

The mean values of grain size and shape-related traits and their BLUP values were used in GWAS analysis to identify MTAs using the GAPIT package of R studio [54]. The HapMap data format was used for genotype files. GWAS was performed by the FarmCPU (fixed and random model circulating probability unification) GWAS method [57]. FarmCPU uses a multi-locus model for testing markers across the genome. This method uses the mixed linear model and stepwise regression model iteratively to eliminate the disadvantages of the general linear model (GLM) and mixed linear model (MLM). In the GWAS analysis, the markers that exceeded the threshold of an FDR (false discovery rate)-adjusted p value at the 0.01 level (−log₁₀ P-values ≥ FDR) were considered to have significant associations with the related traits. Manhattan and Q–Q plots were also created and illustrated in the GAPIT package. DArTseq markers were distributed on the durum chromosome by the rMVP package in R Studio [58].

2.6. Candidate Gene Identification

The flanking sequences that cover 1 Mb upstream and downstream of significant environmentally stable markers (MTAs) obtained from GWAS results were screened against the Triticum turgidum genome (Svevo.v1) using the Ensembl database platform’s BioMart tool (https://plants.ensembl.org/biomart) to find any possible candidate genes. Any match within the 2 Mb flanking region of the MTAs was listed with the gene stable ID, gene starts and ends, and gene description. Additionally, we have benefited from the published literature to know the detailed roles of plausible genes in plants, especially in wheat.

3. Results

A durum wheat panel consisting of 146 genotypes was evaluated for grain size and shape-related traits to identify significant MTAs through GWAS. The preliminary phenotypic results revealed significant genetic diversity for the studied traits among the genotypes. It was appropriate to use this population variation in the following GWAS processes (Table 1).

3.1. Phenotypic evaluation of grain traits

Grain size and shape traits were evaluated using the image analysis software Smart Grain 1.2 [48]. Almost completely normal distributions were observed for all measured traits in the durum population in the BLUP data and the individual environments (Figure S1). The descriptive statistics for all traits over the two environments are shown in Table 2. Broad-sense heritability values were ranged between 0.96 and 0.98. The coefficient of variation percentages were 3.37 (AS), 1.84 (PL), 1.88 (L), 2.16 (W), 2.04 (LWR), 1.23 (CS), and 3.77 (TGW) (Table 2). The variance analyses demonstrated significant differences (P<0.001) among the genotypes for all calculated traits (Table 1). Significant environmental effects were observed in AS, PL, W, LWR, and TGW. A significant genotype-by-environment interaction was observed for all traits (Table 1). With rare exceptions, significant positive and negative correlations were observed between the traits in the two environments and the BLUP data (Table S2). By the way, significant positive correlations were calculated between environments and BLUP data for every single trait, except for E1 vs. BLUP for TGW (Table 3).

3.2. Structure of durum population and SNP density on the genomes

As a result of the strict refinement process, 3,251 high-quality SNP markers were obtained and distributed in the A and B genomes of durum wheat. Marker coverage of the A (2,112 markers) genome was 64.96% of all markers, whereas it was 35.04% in the B (1,139 markers) genome. The A genome chromosomes contained 260 (1A), 454 (2A), 307 (3A), 143 (4A), 282 (5A), 276 (6A), and 389 (7A) markers, whereas the B chromosomes contained 186 (1B), 199 (2B), 173 (3B), 96 (4B), 164 (5B), 155 (6B), and 167 (7B) markers. The number of markers on the homologous chromosome in group 2 was highest (20.08 percent), while those on the homologous chromosome in group 4 were lowest (7.35 percent).

The population was divided into three subpopulations based on the eigenvector principal component analysis. There was a sharp drop in the second principal component, but a significant drop was still observed in the third principal component. Accordingly, we chose the third principal component to cluster the populations (Figure 1A and 1B). The heatmap also verified the clusters using a dendrogram (Figure 2).

3.3. Linkage disequilibrium analysis

LD was calculated using 3,251 SNP markers. Of the 161,276 marker pairs, 106,939 showed a significant linkage disequilibrium at p<0.01 level, which corresponds to 66% of marker pairs, whereas 54,337 marker pairs had r²>0.1. LD decay was estimated based on the r² values for the whole durum wheat genome. The LD between the marker pairs decayed at r²=0.2 value. The drop point of the LD decay was 3,601,053 bp in the whole genome (Figure 3).

3.4. Genome-wide association analysis

Genome-wide association analysis detected 41 MTAs for all evaluated traits except grain width (W). These MTAs were distributed on chromosomes 1A, 2A (21), 3A (2), 4A, 5A (2), 6A (3), 7A (2), 1B (2), 5B, 6B, and 7B (5) of durum wheat (Table S3, Figure 4). Of the 41 MTAs, only 11 were stable across environments (Table 4). These results are almost completely in agreement with the environmental correlations (Table 3). Therefore, we are going to focus on these MTAs in the following sections.

For grain area size (AS), nine MTAs were identified on chromosomes 2A (4), 4A, 6A, and 7B (3). Here, only QAS.su.2A1 was stable in two environments. For grain perimeter length (PL), nine MTAs were detected on chromosomes 1B, 2A (4), 5A, 6A, 7A, and 7B. QPL.su.1B1 was the only stable MTA in three environments. For grain length (L), eight MTAs were found in chromosomes 1B, 2A (4), 5A, 6A, and 7B. Only one MTA, QL.su.1B1, was stable in all three environments for this trait. For the grain length-width ratio (LWR), three MTAs were identified on chromosomes 2A (2) and 3A. Two MTAs on chromosome 2A, QLWR.su.2A1and QLWR.su.2A2, were stable in two environments. Regarding grain circularity (CS), five MTAs were discovered on chromosomes 2A (3), 3A, and 7A. Three MTAs, QCS.su.2A1, QCS.su.2A2, and QCS.su.7A1 were stable in two environments. Seven MTAs were found on chromosomes 1A, 2A (4), 5B, and 7B for the thousand-grain weight (TGW), the trait that is influenced by all the traits mentioned above. For this trait, the QTGW.su.2A2, QTGW.su.2A3, and QTGW.su.7B1 were stable in two, three, and three environments, respectively. Although there were 41 different MTAs, some were identified on the same SNP markers, and these SNPs seem to be associated with multiple grain traits in durum wheat. For example, SNP-1095449 was found to be associated with AS, PL, and L (Table S3). Other markers were SNP-991434 (AS, PL, L, and TGW), SNP-1127014 (AS, PL, and L), SNP-1006957 (AS, PL, and L), SNP-1091721 (AS and PL), SNP-10983760 (PL and L), SNP-991737 (L, LWR, and CS), SNP-1150369 (LWR and CS), and SNP-1127543 (LWR and CS) (Table S3).

3.5. Putative candidate genes underlying grain size and shape-related traits in durum wheat

The LD blocks were created for the stable MTAs to decide the BLAST search border on the Svevo genome (Figure S2). However, we had to use only the flanking sequences that span 1 Mb upstream and downstream (a total of 2 Mb) of the markers because the LD block intervals were very large (up to 55,532 kb) and corresponded to a large number of candidate genes that are far away from the peak markers (Table 5). However, we are aware that LD blocks contain markers with substantial linkage disequilibrium and high r², which is always significant to take into account. In this context, a BLAST search against the Svevo genome detected 118 high-confidence putative candidate genes for all grain traits (Table S4).

The identified genes in the 2 Mb region encode different proteins that have many functions in plants, such as growth and development, stress responses, cell elongation, and seed germination and development. Eighteen of these were found to be associated with seed traits in different plants, especially in wheat (Table 6). Therefore, only the genes associated with seed size and shape will be focused on here. In this context, for AS, QAS.su.2A1 was found in the genes TRITD2Av1G047210 and TRITD2Av1G047390, which encode UDP-glycosyltransferase and glycosyltransferase, respectively. For PL and L, QPL.su.1B1 and QL.su.1B1 were found on the same marker, SNP-100083695, and both MTAs were found in the same genes: TRITD1Bv1G011760, TRITD1Bv1G012100, TRITD1Bv1G012160, TRITD1Bv1G012200, and TRITD1Bv1G012290. The first two genes encode protoheme IX farnesyltransferase and ubiquitin carboxyl-terminal hydrolase, respectively, whereas the last three genes encode histone deacetylase. For LWR, QLWR.su.2A1 was found in TRITD2Av1G065030, which encodes BRI1-EMS suppressor 1 (BES1)/brassinazole-resistant 1 (BZR1) family (BES1/BZR1 homolog 1). The other MTA, QLWR.su.2A2 was found in TRITD2Av1G180930 and TRITD2Av1G181430, which encode transcription factors and digalactosyldiacylglycerol synthases, respectively. For CS, QCS.su.2A1 was detected together with QLWR.su.2A1 on the same marker, SNP-1150369; hence, it coincided with the same gene, TRITD2Av1G065030. Additionally, QCS.su.2A2 was found together with QLWR.su.2A2, whose coincident genes were already mentioned above. Other stable MTA, QCS.su.7A1 found in TRITD7Av1G256220, which encodes B3 domain-containing protein. For TGW, QTGW.su.2A2 was found in TRITD2Av1G048230, TRITD2Av1G048320, and TRITD2Av1G048480, which encode cytochrome P450, patatin, B3 domain-containing protein, respectively. The other MTA, QTGW.su.2A3 was found in TRITD2Av1G191770 and TRITD2Av1G191850, which encode phospholipase C and pentatricopeptide repeat-containing protein. The last stable MTA, QTGW.su.7B1 was found in TRITD7Bv1G159220 and TRITD7Bv1G159310, which encode elongation factor like protein and ABC transporter B family protein, respectively.

4. Discussion

This study was designed to evaluate the grain size and shape traits of a durum wheat panel with 146 genotypes, and to conduct a GWAS to identify MTAs related to grain size and shape-related traits.

4.1. Phenotypic evaluation

After data curation and basic statistical analyses, high genetic diversity was identified among the genotypes (Table 1), and all traits showed a high heritability in the population (Tables 2). The majority of trait comparisons revealed a positive correlation, indicating that all traits had an increasing impact on TGW and consequently increased grain weight (Table S2). All traits showed normal distribution, indicating that multiple genes may be responsible for controlling the grain size traits (Figure S1).

4.2. MTAs identified for grain size and shape traits

Several GWAS studies have been conducted on grain size traits in durum [31,42] and bread wheat [33,34,35,36,37,38,39,40,59,60,61]. Additionally, a few diploid wheat species such as Triticum monococcum [43], Triticum urartu [44], and Aegilops tauschii [45,46]. In the present study, 41 MTAs were identified; however, only 11 MTAs were stable across the environments (Table 4). Therefore, only stable MTAs on chromosomes 2A, 1B, and 7A were compared with previously reported MTAs.

For AS, Gao, et al. [33] reported a significant MTA in bread wheat on chromosome 2A at location 742,132,445 bp. The stable MTA for the current study, QAS.su.2A1, was found to be at a distance of 104,655,222 bp, which makes both of them distinct from one another. In another study, Schierenbeck, et al. [61] reported a significant MTA in bread wheat on chromosome 2A at the position of 82,350,302 bp. This MTA is very close to QAS.su.2A1 (104,655,222 bp). Rabieyan, et al. [59] reported a significant MTA on 2A in bread wheat. However, the position of this MTA was reported in cM length; therefore, its short sequence was searched against the bread wheat genome and was found to be 17,954,870 bp, which is relatively close to QAS.su.2A1. In addition to these, Yu, et al. [26] reported a crucial MTA on 2A in einkorn wheat. Nevertheless, this MTA was also identified as a cM length; therefore, it could not be compared to QAS.su.2A1. In a recent einkorn study, Sesiz, et al. [27] reported two QTLs at around the tip of two arms of chromosome 2A, which were 34,773,385–53,795,616 bp (1) and 581,712,653–600,943,973 bp (2), respectively. Here, QAS.su.2A1 is incontrovertibly close to the first QTL of einkorn wheat. These findings collectively suggested that QAS.su.2A1 may target a genetic area linked to durum wheat’s grain characteristics.

For PL, no MTAs were encountered on 1B in previous reports. However, in the present study, QPL.su.1B1 was identified on 1B in all three environments at 29,327,461 bp. This makes QPL.su.1B1 unique and a new region for grain perimeter length in wheat. Interestingly, for L, the only stable MTA, QL.su.1B1, was detected on the same marker (SNP-100083695) as the QPL.su.1B1. This is not surprising because these traits are known to be highly associated and demonstrated a high positive correlation in this study (Table S2).

For L, Li, et al. [35] reported two important MTAs on chromosome 1B in bread wheat. One MTA was found at 642.6–642.7 Mb, which is far from QL.su.1B1 (29,327,461 bp), whereas the other MTA was found at 26.9–30.8 Mb, which covers QL.su.1B1 completely. This provides strong evidence that QL.su.1B1 may associated with grain length in wheat. In a different study, Muhammad, et al. [37] reported a significant MTA in bread wheat on 1B at 637.0 Mb, which overlapped with the MTA reported by, Li, et al. [35]. Based on previous reports and present findings, some of the grain-length-associated regions appear to be in the distal region of the short and long arms of chromosome 2A in wheat.

For LWR, two MTAs were identified in the present study, namely QLWR.su.2 and QLWR.su.2A2 at 148,130,749 bp and 505,958,255 bp on 2A, respectively. In a previous study, Gao, et al. [34] reported an important MTA on 2A at 724,513,384 bp. This MTA is far away from our MTAs, which means that our MTAs may be the new genomic region associated with LWR in wheat. No more MTAs were encountered in 2A for this trait in the literature.

For CS, three stable MTAs were identified on 2A (2) and 7A, namely QCS.su.2A1 (148,130,749 bp) QCS.su.2A2 (505,958,255 bp) and QCS.su.7A1 (673,131,697bp). For this traits, Gao, et al. [33] reported one MTA on 2A at 742,132,445 bp, which is relatively near to QCS.su.7A1. ,In another, Sesiz, et al. [27] reported a QTL on 2A at start position 106,445,919 bp in einkorn wheat. Our MTA, QCS.su.2A1 is comparatively close to this einkorn QTL. This might be a clue to the position of CS trait in the wheat genome. No additional MTAs were reported in 2A or 7A for the CS trait.

The most significant grain trait, TGW, is a major determinant of grain yield and is mainly affected by a combination of other grain architecture traits. To date, several GWAS studies have reported many MTAs for TGW. However, since MTAs were detected only on chromosomes 2A and 7B in the present study, we focused only on these chromosomes. In this context, Rasheed, et al. [38] identified two MTAs on 2A and 7B at 10.5 cM and 222.0 cM, respectively. However, they could not be compared to our MTAs because they were reported in cM length. In another study, Li, et al. [35] reported one MTA on 2A at 760.6–760.7 Mb, which is distant from QTGW.su.2A2 (106,204,569 bp) and QTGW.su.2A3 (531,237,720 bp). Besides, Schierenbeck, et al. [61] identified an important MTA on 2A at 82,350,302 bp. This MTA is comparatively positioned near QTGW.su.2A2 (106,204,569 bp). Additionally, two important QTLs were reported on 2A at 33,062,393–59,383,738 bp (1) and 609,490,374–676,558,749 bp (2) in einkorn wheat [27]. Here, QTGW.su.2A2 is close to the first einkorn QTL, whereas QTGW.su.2A3 is close to the second einkorn QTL. In light of this knowledge, when considering the relationship between einkorn and bread wheat in terms of the ancestral A genome. These findings may be pointing out the TGW location on the A genome in wheat species.

4.3. Candidate gene prediction

Seed size- and shape-related traits are highly correlated with each other and are important agronomic traits that determine grain yield in wheat. Therefore, all the identified MTAs are individually or collectively important in terms of finding new genetic pathways to improve grain yield in wheat.

The BLAST results provide 85 high-confidence putative candidate genes related to stable MTAs (Table S4). Here, only the putative genes that play a proven role in plants’ grain, especially wheat, were selected to examine our results (Table 6).

For AS, QAS.su.2A1 is corresponded to UDP-glycosyltransferase. Dong, et al. [62] reported a QTL (GSA1) regulating grain size and abiotic stress tolerance by modulating cell proliferation and expansion, which encodes a UDP-glucosyltransferase in rice. They also announced that overexpression of GSA1 resulted in larger grains. In the present study, QAS.su.2A1 may show the genomic region responsible for the formation of large grains by regulating the grain area size.

For PL, QPL.su.1B1 is overlapped with a genomic region that encodes Protoheme IX farnesyltransferase, Ubiquitin carboxyl-terminal hydrolase, and histone deacetylase. Vergès, et al. [63] reported the function of protein farnesylation in the seed development of Arabidopsis. Ubiquitin plays a role in regulating seed and organ size in plants [64,65,66]. Wang, et al. [67] reported that histone deacetylase interacts with plant steroid hormones, brassinosteroids (BRs), which play a role in many plant characteristics, including seed size. Based on these reports, QPL.su.1B1 and its candidates may regulate the grain size in durum wheat. QPL.su.1B1 were discovered at the same marker with QL.su.1B1 on chromosome 1B. Thus, all putative candidate genes identified for QL.su.1B1 are also valid for this MTA, which provides more clues regarding the role of this MTA in grain size in durum wheat.

For LWR and CS, QLWR.su.2A1 and QCS.su.2A1, and QLWR.su.2A2 and QCS.su.2A2, were detected on the same markers on 2A, separately. Therefore, putative candidate genes for these MTAs were evaluated together. In this context, QLWR.su.2A1, and QCS.su.2A1 coincided with BES1/BZR1 homolog 1. Jiang, et al. [68] reported that brassinosteroid plays an important role in determining the size, mass, and shape of Arabidopsis seeds. The other MTAs, QLWR.su.2A2 and QCS.su.2A2 overlapped with transcription factor, and late embryogenesis abundant (LEA) hydroxyproline-rich glycoprotein. Some transcription factors have been reported to be a group of proteins that regulate grain size [69]. For example, Huang, et al. [70] report that the WIDE AND THICK GRAIN 1 (WTG1) gene functions as a significant factor in determining grain size and shape in rice. Importantly, the orthologous of this gene was described in wheat as TaWTG1 on the short arms of group 7 chromosomes in bread wheat [71]. The other protein, LEA is formed during the late period of seed development, and LEA proteins have been detected in seeds of different crops, fruits, and vegetables to date [72]. The other CS-related MTA QCS.su.7A1 coincided with the B3 domain-containing protein. B3 TFs are plant-specific proteins and were first described and cloned in maize (Zea mays) [73]. In addition, Yang, et al. [74] demonstrated that a B3 TF, namely ZmABI19, plays a role as a grain-filling induction regulator. The reported information supports the potential role of the identified MTAs in regulating grain traits in durum wheat.

The increase in grain size leads to an increase in TGW and, thereby, an increase in grain yield in wheat. In this study, almost all the grain size and shape-related characteristics demonstrated a significant positive correlation with each other, and they played roles in increasing TGW. In the present study, three important and stable MTAs, QTGW.su.2A2, QTGW.su.2A3, and QTGW.su.7B1, were identified for TGW and corresponded to some protein products that were directly associated with grain characteristics in different plants. For example, QTGW.su.2A2 coincided with cytochrome P450 (CYP), patatin, and B3 domain-containing proteins. It is reported that CYP family members regulate seed size in Arabidopsis [75], tomato [76], sweet cherry [77], and soybean [78]. In wheat, Ma, et al. [79] reported a gene, TaCYP78A3, that encodes wheat cytochrome P450 CYP78A3, which is expressed in wheat reproductive organs. Their results show that TaCYP78A3 increases size in wheat. Huang, et al. [80] described the role of patatin in seed size in Arabidopsis. In addition, Liu, et al. [81] reported a patatin-related protein, OspPLAIIIα, and found its role in seed size in rice. The B3 domain-containing protein was found for QCS.su.7A1 and discussed for this MTA above. Despite QCS.su.7A1 and QTGW.su.2A2 located on different chromosomes, both of which coincided with the same protein products. These results support a plausible role of the B3 domain in grain architecture in wheat.

The other MTA, QTGW.su.2A3 overlapped with phospholipase C and pentatricopeptide repeat-containing protein. The role of phospholipase C was studied by Yu, et al. [82] who reported that phospholipase C1 modulates grain size in rice. Yang, et al. [83] reported the function of pentatricopeptide repeat-containing protein EMP9 on maize seed development. In addition, Liu, et al. [84] also reported the role of this protein on maize seeds.

The third MTA, QTGW.su.7B1, coincided with the genes encoding elongation factor like protein and ABC transporter B family protein. Transcript elongation factors (TEFs) play a significant role in the regulation, proliferation, and differentiation of cells, and control different stages of growth processes [85]. A member of TEFs, TaTEF-7A, was reported on chromosome 7A in wheat by Zheng, et al. [85]. This protein showed the highest expression in young spikes and developing seeds, and it was reported that it regulates the grain number per spike in wheat. In another study on durum wheat, Giancaspro, et al. [86] identified a perfect candidate on 5B involved in the determination of grain weight and encoding the protein elongation factor. In the present study, QTGW.su.7B1, was found on chromosome 7B, which is orthologous with 7A. Although these MTAs were reported on different wheat chromosomes, they have consistently coincided with TEFs in wheat genomes. The other overlapped protein, ABC transporter B family proteins are essential for plant development, and they have many functions in seed development [87].

The stable MTAs identified in the present study coincide with the genomic regions that encode protein products that have important roles in regulating plant seed traits. As seen and understood from the literature, some of these MTAs overlap with previously reported genomic regions, whereas some are new regions for grain traits in wheat. In both cases, the identified MTAs demonstrated perfect aspects as candidate genes for related traits. We know that the combination of grain size and shape-related traits serves to increase or decrease grain yield. Increasing TGW without decreasing the seed number per spike is the key to increasing wheat yield. However, TGW is determined by many genes with minor effects. Therefore, not only TGW-related MTAs but also other grain trait MTAs have an important effect on regulating grain yield in wheat.

5. Conclusions

In the current study, a durum wheat panel was assessed for traits related to grain size and shape. Using SNP markers and phenotypic data, a GWAS approach was used to discover marker-trait associations (MTAs) for associated traits. A total of 41 MTAs were detected for grain size- and shape-related traits. Of these, only 11 MTAs were stable across the environments. The positions of these stable MTAs were BLAST searched against the Svevo genome, and 118 high-confidence putative candidate genes were identified. These genes encode different protein products that play important roles in plant growth and development, stress response, cell elongation, and seed germination and development. However, 18 of these were found to be associated with seed traits in different plants, particularly in wheat grains. In summary, based on the information reported and the annotation results, the identified MTAs may possibly regulate grain architecture traits in wheat. Ultimately, the results need to be supported by extended approaches, such as converting SNP markers to Kompetitive allele-specific PCR (KASP) markers and using them in different populations or should be deeply investigated in gene expression studies.