A total of 207 common wheat accessions were used in this study, including both elite cultivars and landraces mainly from China and seven other nations (Supplementary Table 1). These accessions were grown in Yuanyang (35°5′N, 113°97′E), Henan province, a main wheat growing area of China during the 2018–2019 and 2019–2020 cropping seasons in a randomized complete block design. Each genotype was planted in two rows of 2 m length with a 0.3 m row spacing. Field management was conducted based on standard agronomic practices.

At maturity, wheat grains were harvested, dried at 55°C for 24 h, and milled to fine powders. Then, the milled samples were dried at 55°C for another 24 h. After that, 200 mg of dried powder for each sample was digested in 8.0 ml HNO₃ in a microwave reactor with a gradient of temperature from 120°C to 180°C for 30 min. After dilution in deionized distilled water, the Mo concentration was measured by inductively coupled plasma mass spectrometry (ICP-MS, NexION 1,000, Perkin Elmer, United States).

The best linear unbiased predictor (BLUP) of Mo concentration for each accession across the 2 years was calculated using the R package lme4 (Bates et al., 2015). The broad-sense heritability (H²) was calculated using the equation H 2 = σ g 2 / ( σ g 2 + σ e 2 / e ) , where σ²_g is the variance of genotype, σ²_e is the variance of environment, and e is the number of years (e = 2 in this study). Pearson’s correlation coefficient (r) for the grain Mo trait across the 2 years was calculated in R version 4.1.1.

All wheat accessions were genotyped using the Wheat 660 K SNP array. To avoid spurious SNPs, the SNPs with minor allele frequency (MAF) < 0.05 and missing data >10% were removed. After filtering, a total of 224,706 SNPs were used for GWAS. The population structure and kinship matrix of the association panel were calculated as described in a previous study (Liu et al., 2020a).

GWAS for grain Mo concentration was performed by a mixed-model approach with the FaST-LMM (factored spectrally transformed linear mixed models) program (Lippert et al., 2011). The effective number of SNPs (n = 37,350 in this study) was calculated with the GEC software (Li et al., 2012). Accordingly, the Bonferroni value of p threshold of 2.68E-05 (p = 1/n; n is the effective number of SNPs) was used to determine significant SNPs. R package CMplot¹ was used to visualize the Manhattan and quantile-quantile (QQ) plots.

To explore the candidate genes responsible for GMoC, genes in the genomic region of 500 Kb upstream and downstream of each significant SNP consistently identified in two or more sets of data were screened based on published slow LD decay in common wheat (Sun et al., 2017; Zhou et al., 2021a). The potential candidate genes were selected based on their annotation information of IWGSC RefSeq v1.1 and functions of homologous genes in Arabidopsis and Oryza sativa. The phylogenetic tree was constructed by using the Neighbour-Joining (NJ) method in MEGA 11 with a bootstrap of 1,000 replicates (Tamura et al., 2021).

As described in a previous study (Liu et al., 2020a), the grains at 20 days after pollination of each wheat accession were collected for RNA-sequencing. The expression levels (FPKM) of the potential candidate genes were extracted for gene expression level analysis. Moreover, the expression levels of these candidate genes in different tissues were downloaded from the Wheat Expression Brower,² which was powered by expVIP platform (Borrill et al., 2016; Ramirez-Gonzalez et al., 2018). Heatmap of the expression levels of the candidate genes was visualized by R package pheatmap.

Two hundred and seven wheat accessions, including elite cultivars and landraces from eight nations representing abundant genetic diversity, were grown at Yuanyang, China in the year of 2019 and 2020 to evaluate the Mo concentration in mature grains (Supplementary Table 1). As a result, continuous and extensive variations in grain Mo concentration were observed among different accessions in the 2 years. As shown in Figure 1A, the GMoC in wheat ranged from 496.00 to 1941.58 ng/g, with an average of 1072.67 ng/g and a coefficient of variation (CV) of 24.05% in 2020 (Table 1 and Figure 1A). In the year of 2019, the GMoC ranged from 293.36 to 1324.41 ng/g, with an average of 729.66 ng/g and a CV of 24.67% (Table 1 and Figure 1A). In addition, the BLUP analysis revealed that the average GMoC was 901.16 ng/g and ranged from 528.13 to 1329.19 ng/g across the 2 years. GMoC in 2 years and the BLUP value all showed continuous variations and approximately normal distributions (Figure 1B; Supplementary Figure 1). In addition, the broad-sense heritability (H²) of GMoC across the 2 years was 0.497, and the Pearson’s correlation coefficient was 0.62 (p < 0.001; Figure 1C).

The GMoC data of 207 wheat accessions in 2019 (E1), 2020 (E2), and their BLUP values (B) were employed for GWAS with the 224,706 SNPs using the FaST-LMM program. The Bonferroni value of p threshold of 2.68E-05 (−log₁₀ (p) = 4.57) was used to determine significant SNPs. A total of 45, 51 and 71 SNPs were determined to be significantly associated with GMoC in E1, E2, and BLUP, respectively (Figure 2; Supplementary Table 2). Interestingly, 39, 47, and 66 SNPs associated with GMoC in three analyses (E1, E2, and BLUP) were located in the region of 726,761,412–728,409,806 bp on chromosome 2A, which could explain 7.23–15.94% of the phenotypic variation (PVE). The remaining significantly associated SNPs were located on chromosome 7A, 7B, and 7D, explaining 4.46–13.88% of the phenotypic variation. In BLUP dataset, AX-108792390, AX-111609105, and AX-109440948 were the most significant SNPs on chromosome 2A, 7B, and 7D, respectively. The GMoC was significantly different among wheat accessions carrying different homozygous genotypes of these SNPs (Figure 3; Supplementary Table 3).

The significant SNPs detected in at least two analyses were defined as common SNPs and used for further exploration of candidate genes. Finally, 52 common SNPs were screened, which were distributed on the chromosome of 2A (48 SNPs), 7B (two SNPs), and 7D (two SNPs) with PVE of 5.76–15.94% (Table 2). The significant SNPs detected in the three sets of data were all located on chromosome 2A. The two common SNPs (AX-111609105 and AX-111031595) on chromosome 7B were detected in E2 and B, while the two common SNPs (AX-109440948 and AX-94482751) on chromosome 7D were identified in E1 and B. The GMoC was significantly different among wheat accessions carrying different genotypes at the 52 SNPs (Figure 3; Supplementary Figure 2).

The genomic regions of 500 kb upstream and downstream of the 52 common significant SNPs were defined as candidate regions and used to explore the candidate genes for GMoC. Genes in three candidate regions on chromosomes 2A, 7B, and 7D were screened, which were located in the intervals of 726.26–728.63 Mb (2A), 611.44–612.46 Mb (7B), and 610.74–612.09 Mb (7D). The potential candidate genes were selected based on their annotation information and functions of the homologous genes in Arabidopsis and Oryza sativa. Finally, three potential candidate genes associated with GMoC were screened out, which were all distributed in the candidate region of chromosome 2A, including TraesCS2A02G496200, TraesCS2A02G496700, and TraesCS2A02G497200 at the position of 727,244,098–727,245,793 bp, 727,914,334–727,918,774 bp, and 728,067,465–728,072,317 bp, respectively (Table 3).

The gene TraesCS2A02G496200 encodes a molybdate transporter 1;2 protein associated with GMoC and was only 132 bp from the significant SNP marker AX-109360792 (Chr 2A: 727,245,925 bp). The gene TraesCS2A02G496700 is annotated as a molybdate transporter 1;1 associated with the GMoC and was close (37.1 kb) to the SNP marker AX-109290174 (Chr 2A: 727,955,888 bp). The SNP marker AX-111628949 (Chr 2A: 728,070,767 bp) on chromosome 2A was located in the intron of the TraesCS2A02G497200 gene encoding a molybdopterin biosynthesis protein CNX1, and the SNP marker AX-111044883 (Chr 2A: 728,072,564 bp) was close (247 bp) to the 5′-untranslated region.

We then analyzed the expression levels of these candidate genes in grains at 20 days after pollination in 207 wheat accessions. Surprisingly, TraesCS2A02G496700 showed extremely low expression in grains (Supplementary Figure 3) and exhibited no significant difference in expression among wheat accessions carrying different genotypes of the SNP AX-109290174 (Figure 4B; Supplementary Table 4). In contrast, both TraesCS2A02G496200 (AX-109360792) and TraesCS2A02G497200 (AX-111628949 and AX-111044883) showed significant differences in expression among different genotypes of the closely associated SNPs (Figure 4; Supplementary Table 4). GMoC tended to be consistent with TraesCS2A02G496200 expression at the AX-109360792 site and inverse to TraesCS2A02G497200 expression at the AX-111628949 and AX-111044883 site (Figure 4; Supplementary Figure 2).

When grown in acidic soils, wheat tends to have low grain yield, quality, and Mo content due to Mo efficiency (Yu et al., 1999; Chatterjee and Nautiyal, 2001; Modi, 2002; Kaiser et al., 2005). Compared with the application of Mo fertilizer, the application of wheat seeds with high Mo concentrations is more economical and environment-friendly to solve the problem of Mo deficiency in acidic soils (Brennan and Bolland, 2007). Therefore, breeding of wheat cultivars with high GMoC may be an effective approach to overcome the Mo deficiency for wheat in acidic soils. Our results showed that GMoC had great variations among the 207 wheat accessions, ranging from 293.36 to 1324.41 ng/g in 2019 and from 496.00 to 1941.58 ng/g in 2020 (Table 1). There was a moderate heritability (H² = 0.497) between the 2 years, suggesting that GMoC is affected by both genetic and environmental factors. This seems to be consistent with previous findings in Arabidopsis (Baxter et al., 2008). Norton et al. (2014) and Yang et al. (2018) reported an even higher heritability of Mo accumulation in rice grains. Our results revealed that the GMoC between 2019 and 2020 was significantly positively correlated with each other (r = 0.621; p < 0.001), indicating that GMoC is relatively stable across different years and genetic factors are important determinants. Therefore, the wheat varieties with stably high GMoC like Taishan 5, Xinmai 18, and Bainong 160 have the potential to be utilized in future wheat breeding programs.

The bi-parental QTL mapping in a previous study identified one QTL (Qgmo.tamu.3B.540) associated with GMoC on chromosome 3B (Yu et al., 2021). However, no GMoC-related genes have been identified by GWAS in common wheat so far. In the present study, we identified 52 common SNPs significantly associated with GMoC through GWAS, which are distributed on chromosome 2A, 7B, and 7D. Interestingly, 48 out of the 52 common SNPs were located on chromosome 2A and distributed in the region of 726,761,412–728,132,521 bp. Based on gene functional annotations, TraesCS2A02G496200, TraesCS2A02G496700, and TraesCS2A02G497200 on chromosome 2A were identified as potential candidate genes.

TraesCS2A02G496700 and TraesCS2A02G496200, which are annotated as molybdate transport 1;1 and molybdate transport 1;2, respectively, are two specific molybdate transporters and belong to the Molybdate Transporter 1 (MOT1) family. In eukaryotes, the MOT1 family mediates high-affinity and specific molybdate transport (Tejada-Jimenez et al., 2013). The orthologous genes of MOT1 have been identified in different species, such as Chlamydomonas reinhardtii (Tejada-Jimenez et al., 2007), A. thaliana (Tomatsu et al., 2007; Baxter et al., 2008; Gasber et al., 2011), rice (Norton et al., 2014; Yang et al., 2018; Huang et al., 2019; Wang et al., 2020; Ishikawa et al., 2021), maize (Asaro et al., 2016), strawberry (Liu et al., 2020b), Lotus japonicus (Gao et al., 2016; Duan et al., 2017), and Medicago truncatula (Tejada-Jimenez et al., 2017; Gil-Diez et al., 2019). The phylogenetic tree showed that TraesCS2A02G496700 is the homolog of LOC_Os08g01120 (OsMOT1;1) in rice, while TraesCS2A02G496200 is the homolog of LOC_Os01g45830 (OsMOT1;2) in rice and AT1G80310 (AtMOT1;2) in Arabidopsis (Supplementary Figure 4). MOT1 is responsible for molybdate uptake, translocation, and accumulation (Tomatsu et al., 2007; Baxter et al., 2008; Gasber et al., 2011; Tejada-Jimenez et al., 2013). Yang et al. (2018) attributed the variations of Mo accumulation in rice grains to changes in the expression level of OsMOT1;1. However, TraesCS2A02G496700 showed quite low expression in wheat grains (Supplementary Figure 3), which could be verified by the public transcript data from expVIP platform (Supplementary Figure 5). The public data demonstrate that TraesCS2A02G496700 is highly expressed in roots. Huang et al. (2019) revealed that the expression level of OsMOT1;1 in roots affects the Mo concentration in rice grains. Therefore, TraesCS2A02G496700 is considered as a candidate gene in this study.

In A. thaliana, AtMOT1;2 (formerly named AtMOT2) is a vacuolar molybdate transporter and involved in inter-organ Mo translocation as well as Mo accumulation in seeds (Gasber et al., 2011). Ishikawa et al. (2021) identified OsMOT1;2 as a vacuolar molybdate export protein that plays an important role in inter-organ Mo distribution in rice. The deletion of OsMOT1;2 decreased the grain Mo concentration in the rice mutant osmot1;2. In our study, TraesCS2A02G496200 is considered as the ortholog of MOT1;2 in wheat (Table 3; Supplementary Figure 4). At the AX-109360792 locus (Chr 2A: 727,245,925 bp; close to the upstream of TraesCS2A02G496200), the accessions with the AA genotype had both significantly higher GMoC and TraesCS2A02G496200 expression than those with the GG genotype (Figure 4A; Supplementary Figure 2), suggesting that it is a candidate gene of GMoC. However, the expression of TraesCS2A02G496200 was the highest in leaves and shoots, followed by roots, while the lowest in grains (Supplementary Figure 5). This is similar to the expression pattern of OsMOT1;2 in rice and AtMOT1;2 in A. thaliana, suggesting its role in inter-organ Mo translocation and distribution (Gasber et al., 2011; Ishikawa et al., 2021). Therefore, it is of great significance to explore the molecular mechanism for the regulatory effect of these two candidate genes (TraesCS2A02G496700 and TraesCS2A02G496200) on Mo accumulation in wheat grains in future studies.

The third potential candidate gene TraesCS2A02G497200 is annotated as molybdopterin biosynthesis protein CNX1 (cofactor for nitrate reductase and xanthine dehydrogenase 1), which is involved in Moco biosynthesis in plants by inserting Mo into molybdopterin (Schwarz et al., 2000; Kuper et al., 2003; Llamas et al., 2004, 2006; Tejada-Jimenez et al., 2018). Significant differences in phenotype and expression were observed for the genotypes at two SNPs (AX-111628949 and AX-111044883) closely associated with TraesCS2A02G497200. Interestingly, GMoC tended to have an inverse relationship with the gene expression at the two sites (Figure 4; Supplementary Figure 2). This phenomenon may be attributed to the conversion of more molybdate into Moco under the catalysis of CNX1, which leads to a decrease in GMoC. CNX1 has been found to be constitutively expressed in all organs of Arabidopsis plants (Schwarz et al., 2000). Similarly, the expression of TraesCS2A02G497200 was also found in almost all organs of wheat plants and was higher in roots than in grains (Supplementary Figure 5). This might be due to a large amount of Mo absorbed by roots from the soil, which is then used for MoCo synthesis. Thus, the above-mentioned three genes are considered as potential candidate genes for GMoC in wheat, which may be utilized in wheat breeding. However, further functional studies are required to verify their functions.