A total of 768 bread wheat accessions were used to identify QTLs of grain size, including 683 Chinese resources (560 improved varieties and 123 landraces) and 85 introduced accessions. Field experiments were performed at three locations, i: the Shandong Agricultural University Agronomy Experimental Station in Tai’an from 2016 to 2019, ii: Weifang Academy of Agricultural Sciences in Weifang in 2019, and iii: Jining Academy of Agricultural Sciences in Jining in 2019. Each accession was planted in five-row plots with 5 cm distance between plants and 25 cm distance between rows. The interval between adjacent plots was 50 cm. At the mature stage, we harvest 10 spikes without any mechanical damage, disease or insect infestation. After threshing, we measured thousand grain weight (TGW), grain length (GL), grain width (GW), grain area (GA), grain perimeter (GP), grain roundness (GR), grain diameter (GD), length-to-width ratio (LWR), grain protein content (GPC) and grain starch content (GSC) for each accession using a Crop Grain Appearance Quality Scanning Machine (SC-E, Wanshen Technology Company, Hangzhou, China).

Genomic DNA was extracted from the seedling leaves of all 768 wheat accessions, followed by further genotyping via an Illumina 55K assay. Finally, a total of 47,743 of 53,063 SNPs were identified in the wheat panel. We estimated the whole-genome distribution and minor allele frequency (MAF) of these SNPs using an in-house Python script. Additionally, we performed quality control of SNPs to exclude those with high missing rate (> 50%) and low MAF (< 5%) for further analysis.

We first extracted 45,298 SNPs with miss rate ≤ 0.5 and MAF ≥ 0.05 from 53,063 SNPs using an in-house Python script. Using PLINK (window size 50, step size 50, r ² ≥ 0.3), a total of 4,360 independent SNPs were further screened out based on r ² of LD ≤ 0.3 (Purcell et al., 2007). The software STRUCTURE was used to calculate varying levels of K (K = 1-20) (Pritchard et al., 2000). We also performed principal component analysis (PCA) and kinship analysis using these independent SNPs and GAPIT software (Lipka et al., 2012; Tang et al., 2016). The phylogenetic analysis of qGL-2D was performed by generating a neighbor-joining tree using Mega 7 (Kumar et al., 2016).

Only 45,298 un-imputed SNPs with miss rate ≤ 0.5 and MAF ≥ 0.05 were used to conduct GWAS for GL, GW, GA and TGW, respectively. The first three PCs were used to construct the PC matrix. We performed GWAS with a Compressed Mixed Linear Model (CMLM) via PCA and kinship analysis using default settings of GAPIT (Lipka et al., 2012; Tang et al., 2016). Additionally, the threshold to determine significant association was set at 1.0 × 10^-5 after Bonferroni-adjusted correction (Pang et al., 2020).

Gene expression data from different wheat cultivars were used to analyze the gene expression profiles of the candidate region. Expression data were download from wheat-URGI website (https://wheat-urgi.versailles.inra.fr/Seq-Repository/Expression). Then the transcriptomic information of candidate genes were exacted by a custom python script. Epidermal tissues were peeled off using tweezers under a stereomicroscope. Then, the cell layers were stained with safranin and mounted on glass slides (Matsunami Glass Ind., Japan). The tissue specimens were subjected to observation with a light microscope (BX50F Olympus Optical Co., Ltd, Japan).

In order to identify candidate genes for qGL-2D, LD heatmaps surrounding peaks were constructed using the R package “LD heatmap” (Shin et al., 2006). Using pairwise LD correlation (r ² > 0.6), we mined the candidate regions of qGL-2D (Yano et al., 2016). We further investigated the expression of these candidate genes in bread wheat grain using typical materials belonging to different haplotypes.

To identify genetic loci associated with grain weight, a panel of 768 bread wheat accessions were constructed, including 560 improved varieties, 123 landraces and 85 introduced accessions. Using a 55K SNP assay, we obtained 47,743 SNPs of the panel. Subsets of these data were further filtered and used in additional analyses ( Figure S1 ). A reasonable assessment of population structure facilitates the identification of real marker-trait associations (Crowell et al., 2016; Juliana et al., 2019). Therefore, we calculated varying levels of K means using un-imputed SNPs and STRUCTURE software (Golbeck, 1987). Landrace, improved and introduced varieties appeared clearly at K = 3 ( Figure 1A ). Further PCA indicated that top three PCs accounted for 17.09%, 6.15% and 3.38% of genetic variation within the bread wheat panel ( Figure 1B ). The results suggested obvious genetic differentiation between landrace and improved varieties of bread wheat.

A total of 10 traits were identified in three different plots for two years, including eight grain shape components (TGW, GL, GW, GA, GP, GR, GD, and LWR) and two grain quality components (GPC and GSC). All traits showed high heritability from 89.30% (GSC) to 95.27% (TGW) ( Table S1 ). After obtaining the best linear unbiased prediction (BLUP) of each accession with respect to each trait across all traits, the coefficient of variation (CV) of all traits ranged from 1.44% GSC to 15.48% TGW ( Table S1 ). GPC was proved to be negatively correlated with the eight grain size components, suggesting that larger, heavier and longer bread wheat grains usually had lower GPC ( Figure S2 ). During the domestication of landrace to improved variety, bread wheat grains increased in size, weight, and length, but their GPC decreased ( Figures 1C, D , 2C ).

Focusing on four key grain shape traits (GL, GW, GA and TGW), GWAS was performed to identify QTLs based on their respective multi-year and multi-location data and BLUP. A total of 23 QTLs were detected on 12 chromosomes, including 5, 6, 6 and 6 QTLs for GL, GW, GA and TGW, respectively ( Table 1 and Figures S3 , S4 ). Seven of 23 QTLs showed common association within different locations or years, including qGW-2B, qGL-2D, qGW-2D.1, qTGW-4A, qTGW-5A.1, qGA-6D, qTGW-6D and qTGW-7D. Consistent with the positive correlations between GL, GW, GA and TGW ( Figure S2 ), close linkage, and overlapping or one-factor-to-many-effects (pleiotropy) were detected on chromosome 2D (for qGA-2D and qGL-2D), chromosome 5A (for qGA-5A, qTGW-5A.1 and qGL-5A.1), chromosome 6D (for qGA-6D and qTGW-6D), and chromosome 7D (for qGA-7D, qGW-7D and qTGW-7D ( Table 1 ).

To validate the results of GWAS, we compared the localization of the QTLs identified in this study with previously detected QTLs associated with bi-parental mapping population. Twelve of 23 QTLs in this study were co-localized with previously reported QTLs, including 1, 6, 3 and 6 QTLs for GL, GW, GA, and TGW, respectively ( Table 1 ). The qGA-6D and qTGW-6D were detected most frequently (five times), followed by qGA-5A, qTGW-5A.1, qGL-5A.1, qGL-5A.2, qGA-7D, qGW-7D and qTGW-7D (twice), whereas qGW-2A, qGW-2B, qGW-2D.1, qGW-2D.2, qGW-3D, qTGW-4A, qGL-5A.2, qTGW-5A.2 and qTGW-5B were detected rarely (once). Additionally, we also identified six new grain size QTLs, including qGA-1D.1, qGA-1D.2, qGA-2D, qGL-2D and qGL-4B.

The qGL-2D was a key locus for GL, as it was detected using the data for each location every year and BLUP ( Figures 3A, B and Figure S3 ). Using BLUP of GL yielded five significant SNPs (-log(p) > 5) representing qGL-2D. Thus, the five SNPs were identified via qGL-2D haplotype analysis. A total of seven haplotypes were detected, including two high-frequency haplotypes (HAP1 and HAP4, 36.6% and 56.4%), two low-frequency haplotypes (HAP2 and HAP3, 3.6% and 2.8%) and three rare haplotypes (HAP5-7, < 1%) ( Figure 3C ). Among them, GL was the shortest in HAP1 (6.56 mm), followed by HAP2 (6.57 mm) and HAP3 (6.70 mm), whereas HAP4 had the longest GL ( Figure 3C ). For other five traits were related to grain shape (GA, GW, GD and HGW) and grain quality (GPC). The HAP4 exhibited the greatest GA, GW, GD, and HGW, and the lowest GPC ( Figure 3D ). The results suggested that qGL-2D was widely involved in grain shape and grain quality.

To determine the evolutionary features of qGL-2D, we conducted a phylogenetic analysis of the seven haplotypes. Two major clades were formed ( Figure 3E ). One clade contained a widely divergent group, including HAP4, HAP3, HAP2 and HAP7, the most prevalent haplotypes associated with improved varieties of bread wheat. Another major haplotype in bread wheat landrace, HAP1, was clustered in the other clade ( Figure 3E ). In summary, the qGL-2D allele associated with improved varieties of bread wheat showed substantial genetic differences compared with bread wheat landrace, which could be attributed to selective effects on large grain during the process of modern bread wheat improvement.

To analyze the candidate gene within qGL-2D, we defined the QTL region based on local LD. As indicated in the LD heatmap, an interval from 522,544,495 to 533,987,666 bp on chromosome 2D was an LD block with r ² > 0.6 ( Figure 2A ). The qGL-2D contains 125 annotated genes. To further reduce the candidate number, we performed transcriptome analysis using one short-grain accession (Chinese Spring (HAP1)), two long-grain bread wheat accessions (Aikang 58 (HAP4), 04chu122 (HAP5)) and 5 BC₂ near isogenic lines (NILs) carrying qGL-2D 04chu122 or aikang58 segment ( Figure 2B ). A total of 29 expressed genes were identified in eight accessions mentioned above ( Table S2 ), and only TraesCS2D02G414800 showed higher expression within two long-grain and eight NILs than in one short-grain accession ( Figures 2D , E ). Homology analysis showed that TraesCS2D02G414800 encodes oleosin, which is involved in seed maturation and germination. Taken together, the results provide possible key candidates for further investigation of the molecular mechanism underlying GL within bread wheat.