Genome-Wide Association Study Reveals Genomic Regions Associated With Molybdenum Accumulation in Wheat Grains

Molybdenum (Mo) is an essential micronutrient for almost all organisms. Wheat, a major staple crop worldwide, is one of the main dietary sources of Mo. However, the genetic basis for the variation of Mo content in wheat grains remains largely unknown. Here, a genome-wide association study (GWAS) was performed on the Mo concentration in the grains of 207 wheat accessions to dissect the genetic basis of Mo accumulation in wheat grains. As a result, 77 SNPs were found to be significantly associated with Mo concentration in wheat grains, among which 52 were detected in at least two sets of data and distributed on chromosome 2A, 7B, and 7D. Moreover, 48 out of the 52 common SNPs were distributed in the 726,761,412–728,132,521 bp genomic region of chromosome 2A. Three putative candidate genes, including molybdate transporter 1;2 (TraesCS2A02G496200), molybdate transporter 1;1 (TraesCS2A02G496700), and molybdopterin biosynthesis protein CNX1 (TraesCS2A02G497200), were identified in this region. These findings provide new insights into the genetic basis for Mo accumulation in wheat grains and important information for further functional characterization and breeding to improve wheat grain quality.


INTRODUCTION
As a critical component of enzymes that catalyze key reactions in nitrogen, carbon, and sulfur metabolism, molybdenum (Mo) is an essential micronutrient required for the growth and development of plants and animals (Arnon and Stout, 1939;Turnlund, 2002;Mendel and Bittner, 2006). Higher plants and animals absorb or take in Mo as oxyanion molybdate, which becomes biologically active by binding to pterin to form Mo cofactor (Moco; Mendel and Bittner, 2006;Mendel, 2013). Thereafter, Moco participates in the synthesis of molybdatedependent enzymes (Molybdoenzymes), including nitrate reductase, sulfite oxidase, xanthine oxidoreductase/dehydrogenase, aldehyde oxidase, and mitochondrial amidoxime reducing component (Hille et al., 2011;Bittner, 2014). To date, a number of molybdate transporters March 2022 | Volume 13 | Article 854966 Jin et al.
Wheat (Triticum aestivum L.) is a major staple food crop worldwide, providing ~25% of calories and nutrients in human diet (Zhou et al., 2021b). Food products of wheat grains are one of the primary sources of dietary Mo for adults (Novotny, 2011;Novotny and Peterson, 2018). However, Mo deficiency is a widespread problem of micronutrient deficiency in agriculture. In Australia, Mo deficiency generally occurs in large areas of cropland with acidic soils and has been identified as the second most common micronutrient deficiency after Zinc (Holloway et al., 2008). In China, more than 44.6 million ha of arable land is Mo-deficient (Zou et al., 2008;Nie et al., 2014). Wheat plants grown in these Mo-deficient soils exhibit severe Mo deficiency symptoms such as pale green leaves, chlorosis of seedlings, and lower yields (Yu et al., 2010). Therefore, dissecting the genetic basis and the molecular mechanism for Mo accumulation in wheat grains will greatly help to improve wheat yield and quality.
In previous studies, to understand the genetic mechanism for Mo accumulation in plants, a number of loci and genes related to Mo accumulation have been identified in several species using the genome-wide association study (GWAS) strategy. For example, molybdate transporter 1 (MOT1), which controls the natural variation of Mo concentration in A. thaliana leaves, was identified by GWA mapping (Shen et al., 2012;Forsberg et al., 2015). In rice, Norton et al. (2014) found a significant locus for Mo accumulation in grains on chromosome 8, which contains the MOT1 orthologue. Yang et al. (2018) analyzed the Mo concentration in rice grains by GWAS and demonstrated that variations of Mo concentration in rice grains can be attributed to variable expression of OsMOT1;1. Furthermore, Cobb et al. (2021) also found strong genetic signals for Mo concentration in the shoot in the genomic region of the OsMOT1;1 gene on rice chromosome 8 by GWAS analysis. In wheat, GWAS has been employed to analyze the concentrations of several minerals in the grains, such as zinc (Alomari et al., 2018;Velu et al., 2018;Cu et al., 2020;Zhou et al., 2020), iron (Alomari et al., 2019;Cu et al., 2020), and calcium (Alomari et al., 2017). However, no GWAS has been carried out on grain Mo concentration (GMoC) in wheat to our knowledge.
In this study, by using a panel of 207 wheat accessions, we aimed to: (i) explore the genetic variation of GMoC in wheat; (ii) identify the genomic regions associated with wheat GMoC using the GWAS approach; and (iii) identify the candidate genes for the variations of wheat GMoC. The results of the present study may facilitate the development of wheat varieties with improved nutritional quality.

Plant Materials and Growth Conditions
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.

Determination of Grain Mo Concentration
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 3 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).

Statistical Analysis
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 2 ) was calculated using the equation H g g e where σ 2 g is the variance of genotype, σ 2 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.

SNP Genotyping and GWAS
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).
Frontiers in Plant Science | www.frontiersin.org 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 1 was used to visualize the Manhattan and quantile-quantile (QQ) plots.

Candidate Gene Identification
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).

Gene Expression Level Analysis
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, 2 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.

Phenotype Variation of Grain Mo Concentration in Wheat
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 2 ) of GMoC across the 2 years was 0.497, and the Pearson's correlation coefficient was 0.62 (p < 0.001; Figure 1C).
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).

Prediction of Candidate Genes for Grain Mo Concentration in Wheat
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 ( 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).

DISCUSSION
As an essential element, Mo is indispensable for nearly all living organisms (Mendel and Bittner, 2006). Several genes have been identified to be responsible for natural variations of Mo level in plants (Huang and Salt, 2016;Yang et al., 2018;Huang et al., 2019;Whitt et al., 2020). However, the genetic basis for natural variation of Mo level in wheat is still poorly understood. Thus, we carried out a GWAS analysis on the Mo concentration in grains among 207 wheat accessions grown in two consecutive years to dissect the genetic basis for Mo accumulation in wheat grains.

Wide Variations of Mo Concentration in Wheat Grains
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 2 = 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.

Identification of Potential Candidate Genes for Grain Mo Concentration in Wheat
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., 2004Llamas et al., , 2006Tejada-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.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
ZfZ, ZL, and ZW designed the experiment. CL conducted the field experiment. ZjZ and YP investigated the phenotype and revised and edited the manuscript. XJ and SY conducted the GWAS analysis. XJ and YP drafted the manuscript. All authors have read and approved the final manuscript.