GWAS identifies an ortholog of the rice D11 gene as a candidate gene for grain size in an international collection of hexaploid wheat

Phenotypic characterization of grain yield components

To explore components of grain yield in wheat, we measured four phenotypes: grain length (Gle), grain width (Gwi), 1000-grain weight (Gwe) and grain yield (Gyi) over two years at two sites. Those phenotypes are referring only to the international panel of wheat and do not include the Canadian accessions. As shown in Table 1, means (± standard deviation) observed for these traits corresponded to: 3.28 mm (± 1.42) for grain length, 1.77 mm (± 0.88) for grain width, 36.17 g (± 21.77) for 1000-grain weight and 2.30 t/ha (± 1.44) for grain yield. The broad-sense heritability estimates were 90.6% for grain length, 97.9% for grain width, 61.6% for 1000-grain weight and 56.0% for grain yield. An analysis of variance revealed significant differences due to genotypes (G) for all traits and, for two traits (Gwe and Gyi), the interaction between genotype and environment (GxE) proved significant. A correlation analysis showed a high significant positive correlation between grain yield and grain weight (r = 0.94; p < 0.01) and also between grain length and grain width (r = 0.84; p < 0.01). Also, significant positive correlations were identified between grain yield and grain length (r = 0.50; p < 0.01) and between grain yield and grain width (r = 0.43; p < 0.01). Interestingly, a bimodal distribution was observed for grain length and width (Fig. 1). Together, these results suggest that a major gene controls two important characters related to grain size with a high heritability within this collection.

In examining the relationship between 1000-grain weight and grain length/width using bagplots on the collection of 159 accessions, no outliers were found when considering the relationship between grain weight and width. In contrast, two accessions (Attila3, Babax8) were indeed detected as outliers when comparing grain weight and length (Supplementary Fig. S1). In the later steps (analysis of population structure and GWAS) we excluded these two accessions considered to be outliers.

Genome-wide SNP marker discovery and validation

To genetically characterize our wheat collection and study the genetic determinants of grain size, we used a double digestion (PstI/MspI) GBS approach to genotype this collection. Overall, 77,124 and 73,784 SNPs were discovered for the set of 71 Canadian wheat accessions and 157 exotics wheat accessions, respectively.

To assess the reproducibility and accuracy of genotypes called via the GBS approach, we genotyped 12 different plants of CS (i.e. biological replicates), which were added to the set of 288 wheat samples for SNP calling and bioinformatics analysis. Sequence reads of the full set of 300 wheat samples obtained from GBS were analyzed following the standard steps of SNP calling and bioinformatics analysis described below. This yielded a total of 129,940 loci that were used for the assessment of accuracy and reproducibility of SNP calls. For each individual plant of CS, the GBS calls were compared between replicates and with the Chinese Spring reference genome (at the corresponding positions).

On the non-imputed data, we detected a very high level of concordance (99.9%) between the genotypes of each CS individual and the reference alleles for the 1,196,184 called genotypes ([130 K SNPs × 12 samples]—missing data; Supplementary Fig. S2). Among those 12 biological replicates of CS, we found a very high reproducibility of genotype calls, as the pairwise identity of genetic distance calls varied from 1.56E−04 to 5.08E−04, with an average of 2.86E−04. In order to ensure about identity of each CS plant, we have found that this value between the individual w56_Guelph (Canadian wheat variety) and each of the CS plant is greater than 0.1.

After imputation of the missing genotype calls, we observed a mean concordance of 93.8% between the CS individuals and the CS reference genome. Furthermore, 76.7% of genotypes were called initially and 23.3% of genotypes were imputed. It should be noted that the accuracy rate for imputing missing data is 73.4%. More details of SNP data set are provided in supplementary Table S1.

As a further examination of data quality, we compared the genotypes called using both GBS and a SNP array on a subset of 71 Canadian wheat accessions that had been previously genotyped using the 90 K SNP array. A total of 77,124 GBS-derived and 51,649 array-derived SNPs were discovered in these 71 accessions (Supplementary Table S2). Of these, only 135 SNP loci were common to both platforms and among these potential 9,585 datapoints (135 loci × 77 lines), only 8,647 genotypes could be compared because the remaining 938 genotypes were missing in the array-derived data. As shown in Fig. 2, a high level of concordance (95.1%) was seen between genotypes called by both genotyping approaches. To better understand the origin of discordant genotypes (4.9%), we inspected the set of 429 discordant SNP calls and observed that: (1) 3.5% of discordant calls corresponded to homozygous calls of the opposite allele by the two technologies; and (2) 1.4% of discordant calls were genotyped as heterozygous by GBS while they were scored as homozygous using the 90 K SNP array. More details are provided in Supplementary Table S3. From these comparisons, we conclude that GBS is a highly reproducible and accurate approach for genotyping in wheat and can yield a greater number of informative markers than the 90 K array.

Genome coverage and population structure

For the full set of accessions, a total of 129,940 SNPs was distributed over the entire hexaploid wheat genome. The majority of SNPs were located in the B (61,844) and A (50,106) sub-genomes compared to the D (only 17,990 SNPs) sub-genome (Table 2). Although the number of SNPs varied two to threefold from one chromosome to another within a sub-genome, a similar proportion of SNPs was observed for the same chromosome across sub-genomes. Typically, around half of the markers were contributed by the B sub-genome (47.59%), 38.56% by the A sub-genome and only 13.84% by the D sub-genome.

The analysis of population structure for the accessions of the association panel showed that K = 6 best captured population structure within this set of accessions and these clusters largely reflected the country of origin (Fig. 3). The number of wheat accessions in each of the six subpopulations ranged from 6 to 43. The largest number of accessions was found in northwestern Baja California (Mexico) represented here by Mexico 1 (43) and the smallest was observed in East and Central Africa (6).

GWAS analysis for marker-trait associations for grain size

To identify genomic loci contributing to grain size in wheat, we performed a GWAS analysis on 157 accessions (excluding the two accessions considered to be outliers) and 73,784 SNPs. As seen in Fig. 3, both Q–Q plots suggest that the confounding effects of population structure and relatedness were well controlled. For both traits, the greatest marker-trait associations were detected at the end of chromosome 2D, while another weaker association was shared at the beginning of chromosome 1D. For grain width only, a marker-trait association was detected on chromosome 4A. In total, seven SNPs were found to be associated with one or both traits, with respectively one, five and one significant SNPs being located on chromosomes 1D, 2D and 4A. Except for two SNPs (chr2D:442798939 and chr4A:713365388), all other SNPs were significant for both grain length and grain width. The SNP at 4A:713365388 was significant only for grain width while the SNP at 2D:442798939 was significant only for grain length.

The most significant association was observed on chromosome 2D and contributed to both grain length and grain width (Table 3, Fig. 3). For this QTL, a total of four SNPs was observed and the SNP most significantly associated to both traits was located at position 2D:452812899. A fifth SNP located at 2D:442798939 was significantly associated to grain length only, but was just below the significance threshold (p-value = 4.34E−05) for grain width.

A high degree of LD was detected among some of the seven SNPs from chromosome 2D displaying association with grain traits. These formed one discontinuous linkage block as the LD between markers belonging to this block was higher (mean of r² = 0.90). For this reason, we considered these to define one quantitative trait locus (QTL) on chromosome 2D (Supplementary Fig. S3). This QTL included 5 SNP markers (chr2D:403935865, chr2D:442798939, chr2D:444560418, chr2D:452644656 and chr2D:452812899) and the peak SNP (chr2D:452812899) explained between 7 and 13% of the phenotypic variation for grain length and width. The minor allele frequency (MAF) at this locus was 0.31 and exerted an allelic effect from − 0.81 to − 0.35 mm (Table 3).

On chromosome 1D, the SNP marker chr1D:166874041 defined a QTL for both grain length and width. The percentage of phenotypic variation explained by this marker for grain length and width was 11% and 6% respectively, with a MAF of 0.30 and allelic effects of 0.76 and 0.33 mm for grain length and width, respectively. Furthermore, a high degree of interchromosomal LD was observed among the peak SNPs between chromosomes 1D and 2D (r² = 0.94) displaying association with grain traits. In addition, almost all accessions which have the major allele on chromosome 1D are the same which have the major allele on chromosome 2D. Thus, the combined impact of these two loci could explain the observed bimodal distribution.

On chromosome 4A, the SNP marker chr4A:713365388 defined a QTL for grain width only and it explained 6% of the variation, had a MAF of 0.14 and exerted an allelic effect of 0.36 mm. However, we reported a very weak LD between this peak SNP marker and the two others on chromosomes 1D and 2D.

In summary, a total of three QTLs significantly associated with grain length and/or width were identified on chromosomes 1D, 2D and 4A.

Candidate gene detection for grain size

To identify candidate genes contributing to grain size within the studied wheat collection, we investigated the genes residing in the same linkage block as the peak SNP for each QTL. On chromosome 2D, the QTL with the largest number of associated SNPs (chr2D:403935865 to chr2D:452811303) included a total of 315 high-confidence genes of which 66 genes are expressed during embryogenesis and grain development in wheat. On chromosomes 1D and 4A, the linkage blocks harboring SNP markers chr1D:166874041 and chr4A:713365388, each defining a QTL, did not include high-confidence genes. Upon examination of the annotations and gene expression profile for the candidate genes, the most promising appears to be the TraesCS2D01G331100 gene in the QTL on chromosome 2D, which is most highly expressed in the developing embryo during embryogenesis and grain development in wheat (Fig. 4). As well, it is expressed in both endosperm and pericarp, and was found to encode a cytochrome P450 (CYP724B1), which showed homology to enzymes involved in brassinosteroid biosynthesis, indicating the mechanism by which seed size may be regulated in wheat. It is an ortholog of the rice CYP724B1 gene, commonly known as the D11 gene. The D11 gene was previously reported as being involved in the regulation of internode elongation and seed development due to its role in the synthesis of brassinosteroids, key regulators of plant growth promoting the expansion and elongation of cells. More details are provided in Supplementary Table S4.

Haplotypes at the wheat orthologue of the rice D11 gene and their phenotypic effects

To provide a useful breeding tool for the main QTL identified in this research, we defined SNP haplotypes around our candidate gene. Using HaplotypeMiner, we identified two SNPs (chr2D:423365752 and chr2D:425474599, Supplementary Fig. S4) that best captured the SNP landscape in the vicinity of the candidate gene. These markers reside in the same haplotype block as the SNP markers, but were not individually found to be significantly associated with grain width and length. These SNP markers define three haplotypes (AT, CT or CC) around the candidate gene, with 99, 18 and 40 individuals carrying these haplotypes, respectively. To investigate the phenotypes associated with these haplotypes, we analyzed the trait value for each haplotype. Interestingly, we observed that for all traits, the mean values of accessions with haplotype AT were significantly larger (p < 0.001) than those obtained for the other haplotypes. As shown in Fig. 5, accessions carrying haplotype AT showed mean values of 3.76 mm for grain length, 2.02 mm for grain width, 40.87 g for grain weight and 2.55 t/ha for grain yield, compared to 2.16 mm, 1.05 mm, 26.87 g and 1.75 t/ha (respectively for grain length, width, weight and yield) for accessions carrying haplotype CC and 1.65 mm, 0.78 mm, 26.89 g and 1.69 t/ha (respectively for grain length, width, weight and yield) for accessions carrying haplotype CT. Furthermore, the relation between the 3 haplotypes and the 6 groups found in the population analysis showed that the haplotype AT predominates in the populations of Mexico 1 and North Africa (Supplementary Fig. S5, Supplementary Table S5). To conclude, we suggest that SNP markers corresponding to haplotype AT will provide a useful tool in marker-assisted breeding programs to improve wheat productivity. Therefore, we point out that the relationship between yield and haplotypes around the D11 gene would allow the selection of high-yielding wheat lines in a breeding program.