A total of 272 wheat landraces, obtained from 10 major wheat-growing zones in China, were utilized in this study (Supplementary Table S1). All landraces were planted in six environments: Ya’an (103°370 E, 22°420 N) in 2012 (E1), Wenjiang (103°410 E, 30°360 N) in 2013, 2014, and 2015 (E2, E3, and E4), and Chongzhou (103°390 E, 30°330 N) in 2014 and 2015 (E5 and E6). Each row of material was 30 cm apart and 1.5 m long, and contained 15 seeds. Field management referred to criteria commonly practiced in wheat production. Four yield-related traits were evaluated in at least four environments: SL – the average length of the main spikes from five plants; SN – the average total number of spikelets from five main spikes of the plant; TN – the average number of tillers in five plants; TKW – the average weight of five samples of 1,000 kernels randomly selected from a given genotype. Two bi-parental populations [Huimai × Datianquxiaomai (HD) and Huimai × Heshangmai (HH)] were used to validate the results. The parents of bi-parental populations were selected from 272 wheat landraces. The TKW of these two RIL populations was evaluated in Chongzhou in 2019, and Huimai resulted significantly more productive than Datianquxiaomai and Heshangmai (higher TKW).

To eliminate environmental effects, the best linear unbiased prediction (BLUP) value for each trait was calculated across environments and used for statistical analysis. The BLUP was calculated using the methodology of Piepho et al. (2008) as previously described (Liu et al., 2017). The broad-sense heritability (H²) value was calculated using SAS v8.1 (SAS Institute Inc., Cary, NC, United States) and is defined as:

where V_g, V_ge, and V_e are the estimates of genetic variance, the genotype × environment interaction, and environmental variance, respectively (Smith et al., 1998). Correlation analyses were performed using SPSS 20 (IBM SPSS Statistics; IBM Corp., Armonk, NY, United States).

Genomic DNA was extracted from a single plant for each of the accessions using the cetyltrimethylammonium bromide method (Murray and Thompson, 1980). DNA samples with an A260/280 ratio of 1.8–2.0 and a concentration of >80 ng μl⁻¹ were used for genotyping. The 272 accessions had been already genotyped by wheat 660 K array, including 630,517 probes (Axiom^® Wheat660 SNP array; Affymetrix, Santa Clara, CA, United States) by Zhou et al. (2017). The software, bowtie2, was employed in this study to obtain physical locations of 172,711 SNPs based on IWGSC RefSeq v1.0.¹

The population structure was performed by Structure2.3.4 based on the Bayesian clustering algorithm (Pritchard et al., 2000). Ten runs of STRUCTURE were performed with a K between 1 and 10 using the admixture model with 100,000 replicates for burn-in and 100,000 replicates for MCMC. The R package “mrMLM” in R Project for Statistical Computing was used to examine the association between markers and yield-related traits. Marker-trait associations were performed using five multi-locus models, including FASTmrEMMA (Wen et al., 2017), ISIS EM-BLASSO (Tamba et al., 2017), mrMLM (Wang et al., 2016), pKWmEB (Ren et al., 2018), and pLARmEB (Zhang et al., 2017). All five models were adjusted by both population structure and family relationship. A logarithm of odds (LOD) value ≥3 was used as the screening criterion (Guan et al., 2019). According to previous studies, the linkage disequilibrium decay of wheat ranges from 3.5 to 23 Mb (Kidane et al., 2019; Li et al., 2019; Luján Basile et al., 2019). Significant SNPs were therefore selected with a physical distance ≤10 Mb and referred to as a conservative QTL.

To further validate the significant QTL identified for TKW, the peak SNP with this locus was converted into a kompetitive allele-specific polymerase chain reaction (KASP) marker based on the probe sequence. KASP primers were designed and analyzed as described in previous studies (Lin et al., 2020b, 2021) and produced by Sangon Biotech (Shanghai, China; Supplementary Table S2). The KASP marker detected different alleles at this locus in the two bi-parental populations. Eighty-two lines were selected for each of the alleles from both populations. These lines were used to evaluate differences in TKW between the two allele groups using a Student’s t-test in SPSS 20 (IBM SPSS Statistics; IBM Corp., Armonk, NY, United States).

Based on IWGSC RefSeq v1.0, predicted genes within ~10 Mb of the physical location of the QTL were selected. We undertook two different methods to predict the possible existence of candidate genes. The first method was expression analysis. Based on data from WheatExp,² genes expressed highly at stages Z71 and Z75 were the most important, due to the key stages in kernel development (Zadoks et al., 1974). Fragments per kilobase of exon model (FPKM) represented gene expression level. As in Lin et al. (2020b), genes expressing less than two per million mapped reads at any stage were removed. FPKM is fragments read per thousand bases per million mappings and represented gene expression level.

The second method to predict the possible existence of candidate genes was annotation. All genes were also used to perform homologous annotations in rice and Arabidopsis using the KEGG Orthology Based Annotation System 3.0 (KOBAS 3.0; Xie et al., 2011).³ Functional annotations were performed via UniProt⁴ and EnsemblPlants.⁵

The four traits among the 272 wheat landraces varied considerably (Table 1). Based on the BLUP values, the SL ranged from 6.33 to 14.63 (cm). SN ranged from 19.43 to 27.36 (count). TN ranged from 8.11 to 18.10 (count), and TKW ranged from 17.90 to 40.47 (g). The heritability of these four traits ranged from 0.64 to 0.93 (Table 1). TN and TKW showed medium heritability, whereas SL and SN showed high heritability. Correlation analysis showed that all correlations were significant (p < 0.05), except the correlation between SL and SN (Table 2), indicating that the changes in these two traits are independent.

A total of 172,711 polymorphic SNPs were obtained [after filtering for missing values ≤10% and minor allele frequency (MAF) ≥0.05] from the Wheat 660 K SNP arrays. This subset was used to perform GWAS. Based on LOD values ≥3, a total of 308 significant SNP markers were identified using the five multi-locus models. Detailed significant SNP markers information and the number of significant SNP markers detected by different models are shown in Supplementary Table S3 and Figure 1, respectively. In FASTmrEMMA, the least number of significant SNP markers was detected, with only 35. A total of 10, 8, 6, and 11 significant SNP markers were detected for SL, SN, TN, and TKW, respectively, with phenotypic variation explained (PVE) up to 6.55%. In the model of ISIS EM-BLASSO, 73 significant SNP markers were revealed. A total of 20, 23, 15, and 15 were detected for SL, SN, TN, and TKW, respectively, with PVE up to 9.27%. The mrMLM model could detect 67 significant markers. A total of 18, 20, 14, and 15 were detected for SL, SN, TN, and TKW, respectively, with PVE up to 13.09%. In the model of pKWmEB, the most number of significant SNP markers was detected (81), and a total of 23, 25, 19, and 14 were detected for SL, SN, TN, and TKW, respectively, with PVE up to 18.78%. The pLARmEB model detected 52 significant SNP markers, and 12, 15, 11, and 14 were detected for SL, SN, TN, and TKW, respectively, with PVE up to 11.09%. Significant markers existed in the above three models, and no more than 10 Mb was considered as a robust QTL.

Twenty-seven robust QTL were identified by more than three different multi-locus models and were considered as robust QTL (Table 3).

Seven QTL associated with SL were identified on chromosomes 1A, 1D, 2A, 2B, 3A, 4A, and 7B, with the PVE up to 11.44%, and the LOD values ranging from 3.45 to 14.68. QSl.sicau-3A, located at 650.86–651.84 Mb, was identified in four models, with LOD values up to 14.68. Nine QTL associated with SN were identified on chromosomes 1B, 2A, 2B, 3B, 5A, and 7A, with PVE up to 7.61%, and LOD values up to 15.85. QSn.sicau-7A, located at 671.48 Mb, was detected by four models. The highest LOD (15.85) and PVE (7.61%) values on this QTL were detected by the pKWmEB model. Four QTL for TN were identified on chromosomes 1B, 2B, 5B, and 6D, with a PVE up to 9.58%, and the LOD values up to 8.42. QTn.sicau-6D, located at 469.60–472.74 Mb, had the highest PVE value which was 9.58%. Seven QTL associated with TKW were identified on chromosomes 2B, 3A, 4B, 6B, 7A, and 7B, with PVE values up 18.78%, and the LOD values up to 12.04. QTkw.sicau-4B, identified by the pKWmEB model, explained the highest PVE (18.78%). Interestingly, among all robust QTL detected in the present study, QTkw.sicua-4B resulted that with the highest PVE value. It can be regarded as a major QTL, contributing to breed more productive wheat cultivars, so it was validated in the first instance.

To validate the genetic effect of the peak SNP for QTkw.sicau-4B, a KASP marker (KASP-AX-108886949) was developed and the differences between TKW of landraces carrying alternatively A/G allele were calculated, both in natural population and in two RIL populations. Among the 272 landraces tested, the average TKW of landraces carrying genotype allele A was heavier than those carrying allele G (p < 0.01; Table 4; Figure 2A). The difference was 22.87%. In the HD population, genotypes were divided into two groups: 42 with allele A and 40 with allele G. The TKW of genotypes with allele A ranged from 28.75 to 44.04 g, whereas those with allele G ranged from 22.28 to 39.27 g. The average TKW with allele A was heavier than that with allele G (p < 0.01; Table 4; Figure 2B). In the RIL population HH, there were 39 with allele A and 43 with allele G. The TKW of genotypes with allele A ranged from 22.52 to 44.56 (g), whereas those with allele G ranged from 19.03 to 37.47 g. The average TKW with allele A was heavier than with allele G (p < 0.01; Table 4; Figure 2C). The difference in TKW between genotypes ranged from 15.33 to 16.81%, with an average value of 16.07% among the two RIL populations. Thus, the developed KASP marker may be useful for breeding cultivars carrying a high-TKW allele.

We observed superior allele A frequency in 10 major wheat-growing zones (Figure 3). The frequencies show that the superior allele A showed an 8% frequency in Chinese wheat landraces, and more than half the materials from the Northern Spring Wheat Zone and the Northwestern Spring Wheat Zone carried allele A. The Northern Winter Wheat Zone, Yellow and Huai River Valleys Facultative Wheat Zone, Middle and Low Yangtze Valleys Autumn-Sown Spring Wheat Zone, Southwestern Autumn-Sown Spring Wheat Zone, and Qinghai-Tibetan Plateau Spring–Winter also had some materials with allele A. Materials from the Southern Autumn-Sown Spring Wheat Zone, Northeastern Spring Wheat Zone, and Xinjiang Winter–Spring Wheat Zone did not carry allele A.

A total of 46 high-confidence genes which were identified in IWGSC RefSeq v1.0 flanked QTkw.sicau-4B by no more 5 Mb. The search for candidate genes was undertaken following two pipelines. The expression analysis of the candidates was investigated in the two most interesting stages for kernel development: Z71 and Z75. Nine genes were expressed at both stages, namely, TraesCS4B01G272300, TraesCS4B01G272500, TraesCS4B01G273400, TraesCS4B01G273500, TraesCS4B01G274200, TraesCS4B01G274500, TraesCS4B01G275400, TraesCS4B01G275500, and TraesCS4B01G276300. Based on these sequences, we found homologous genes using KOBAS 3.0 (see footnote 3). Three genes were identified as TKW-regulating candidate genes (Table 5). TraesCS4B01G272300, TraesCS4B01G275400, and TraesCS4B01G276200 were homologous to OsMCA1 (Liu et al., 2015c), H2B (Ma et al., 2019a), and EDR1 (Prasad et al., 2009) in rice, respectively.