The association panel consists of 166 representative wheat cultivars, such as 144 accessions from Yellow and Huai Valley, the major wheat zone in China, and 22 from five other countries (Supplementary Table 1). The panel was used in our previous studies on black point (Liu et al., 2017), flour color and polyphenol oxidase (PPO) activity (Zhai et al., 2018, 2020), grain yield related traits (Li et al., 2019), and water-soluble carbohydrate contents (Fu et al., 2020). In this study, the panel was grown in randomized complete blocks with two replications in four environments comprising Beijing (39°56′N, 116°20′E), Gaoyi (37°33′N, 114°26′E), and Shijiazhuang (37°27′N, 113°30′E) in Hebei province, and Urumqi (42°45′N, 86°37′E) in Xinjiang province during the 2019–2020 cropping season. These environments were designated as 20BJ, 20GY, 20SJZ, and 20XJ hereafter, respectively. Each entry with approximately 40–50 grains was grown in 1.0 m long row with an inter-row spacing of 20.0 cm. Standard agronomic practices were performed at each location, along with soil application of 25 kg/ha ZnSO₄.7H₂O granular fertilizer (Sinochem Group Co., Ltd., Beijing, China) in 20GY and 20SJZ over three crop cycles prior to this experiment to enrich the available soil Zn and minimize heterogeneity. The green manure returning was routinely practiced in Beijing and field soil held sufficient Zn intrinsically, so that additional Zn fertilizer was not applied in 20BJ. To investigate GZnC and GFeC in a low-Zn environment, 20XJ was set without soil Zn fertilizer application.

Plant materials were hand-harvested in the field and completely dried grains were hand-threshed and cleaned carefully to avoid the potential contamination of mineral element. Approximately 15 g of each sample was subjected to GZnC and GFeC analyses with an X-ray fluorescence spectrometry (EDXRF) instrument (model X-Supreme 8000, Oxford Instruments plc.), following the protocol of high-throughput screening of micronutrients in whole-wheat grain (Paltridge et al., 2012). PH, TKW, and GA in all environments were investigated following Xu et al. (2019).

ANOVA, Pearson’s correlation analysis, and Student’s t-test were conducted using the SAS 9.4 software (SAS institute, Cary, NC, United States). Broad-sense heritabilities (H²) were calculated using the following equation (Holland et al., 2003):

where σ²_g indicates the variance of genotypes, σ²_ge and σ²_ε represent the variances of genotype × environment interaction and errors, and e and r are the numbers of environments and replicates in each environment, respectively. For each trait, a best linear unbiased estimation (BLUE) using a linear model for each accession was calculated across environments by the QTL IciMapping v4.1 software (Yin et al., 2015).

Genomic DNA was extracted following the CTAB method (Murray and Thompson, 1980) and the population was genotyped by the Illumina wheat 90K (such as, 81,587 SNPs) and Affymetrix wheat 660K (such as, 630,517 SNPs) SNP arrays by CapitalBio Technology Co., Ltd, Beijing, China.¹ Markers with minor allele frequency (MAF) < 5.0% and missing data containing heterozygous genotypes > 20.0% were filtered out and the same final 373,106 high-quality polymorphic SNPs were obtained (Fu et al., 2020). These markers included 359,760 (96.4%) from the 660K and 13,346 (3.6%) from 90K SNP array, indicating SNPs in wheat 90K array have relatively low polymorphisms in Chinese cultivars. Flanking sequences of SNP markers were used for BLAST analysis against IWGSC RefSeq v.1.0,² to obtain physical positions in accordance with the best BLAST hit results. Eventually these markers were integrated into one consensus map, covering a total physical distance of 13.7 Gb, accounting for about 95% of wheat whole genome sequences, and were further utilized for GWAS. Relevant information is available in our previous studies (Liu et al., 2017; Fu et al., 2020).

The 166 accessions were classified into three subpopulations by the population structure using Structure v2.3.4 (Liu et al., 2017). The average linkage disequilibrium (LD) decay distances for A, B, D, and whole genomes were approximately 6, 4, 11, and 8 Mb, respectively (Liu et al., 2017). In this study, GWAS was conducted for each environment and BLUE values of GZnC across 20GY, 20SJZ, and 20BJ. The site 20XJ was excluded from the analysis owning to relatively low Pearson’s correlation coefficients with other environments (data not shown). The low Pearson’s correlation was possibly attributed to soil Zn content heterogeneity in 20XJ with null Zn fertilizer application (Velu et al., 2018). As for GFeC, GWAS was carried out across all four environments and BLUE values. Three models, MLM (Q + K), FarmCPU, and MLMM, were employed for the GWAS in GAPIT a software package operating in R v3.5.1³ (Tang et al., 2016). To maximize chances of identifying all possible QTL, a threshold of p = 1.0 × 10^–3 (−log₁₀ (p) = 3.0) was set for the significance of marker-trait associations (MTAs). To control the false discovery rate (FDR) at an appropriate level, those detected in three or more environments (BLUE value was regarded as one environment hereafter) by any of the three models were considered as reliable MTAs (Fu et al., 2020). MTAs detected in the same LD block or physically closely linked were grouped into a single QTL, and the distance between the two very far flanking markers was seen as the QTL interval. The most significant MTA with the lowest value of p across environments in this interval was selected as the representative marker, and its R² output by MLM was used to reflect the proportion of phenotypic variance explained (Fu et al., 2020). Manhattan and quantile–quantile (Q–Q) plots were generated using the CMplot package in R v3.5.1 software.⁴

According to the gene annotations from IWGSC RefSeq v.1.0 and putative homologs in the UniProt database,⁵ the genes located in or adjacent to the physical intervals of QTL identified in this study were subjected to be screened as follows. (1) Those related to the molecular processes in Zn/Fe homeostasis, such as metal uptake in the root, transport in the xylem and phloem, storage in seeds, and regulatory factors, were predicted to be potential candidates; (2) analysis of orthologs between wheat and model plants was carried out to obtain candidate genes in Triticeae-Gene Tribe browser⁶ (Chen et al., 2020) with most of the genes collected in Tong et al. (2020); (3) haplotype analysis was performed within Wheat SnpHub Portal⁷ for 641 accessions with known genome sequences (Hao et al., 2020; Wang et al., 2020), and the genes without sequence variation were removed from candidacy. The most promising candidate genes were eventually selected. The database expVIP⁸ was used to investigate spatio-temporal transcriptional dynamics of candidate genes, providing expression profiles of these genes across tissues and developmental stages (Ramírez-González et al., 2018).

The representative markers for 3AL and 7AL QTL were chosen and converted to KASP markers. Chromosome-specific primers were designed using Polymarker⁹ (Ramirez-Gonzalez et al., 2015) and KASP assays were performed following Xu et al. (2020). The KASP markers were successfully converted when their genotypic results were the same as the original SNP arrays. The 146 F₆ recombinant inbred lines (RILs) developed from a cross between two modern wheat cultivars Zhongmai 175 and Lunxuan 987 (shorten as ZM175/LX987) with significantly different GZnC were used as validation population for the target QTL. GZnC was determined in the bi-parental population in the same way as mentioned earlier in the natural population. Student’s t-tests were conducted to verify allelic effects based on the phenotypic data from the bi-parental population. For the candidate gene of QTL, its variations between ZM175 and LX987 were extracted from the genome resequencing database at Wheat SnpHub Portal (see text footnote 7; Hao et al., 2020; Wang et al., 2020).

A wide range of continuous variations of GZnC and GFeC were observed for the 166 accessions across all environments with near-normal distributions (Supplementary Figures 1, 2). The BLUE values for GZnC and GFeC were 29.25–50.98 mg/kg and 39.86–54.66 mg/kg with mean values of 38.03 and 45.97 mg/kg, respectively (Table 1). In all the accessions, six cultivars with the highest values for GZnC and GFeC were selected with the BLUE values greater than 50 mg/kg and stable performance across the environments (Supplementary Figure 3). Yumai2 and Xinong1376 were the cultivars with both high GZnC and GFeC. Significant correlations (r) between environments were observed for both GZnC and GFeC (Supplementary Figure 4). The results of ANOVA revealed that the factors of genotypes, environments, and G × E interactions markedly affect the GZnC and GFeC (Supplementary Table 2). Broad sense heritabilities (H²) of GZnC and GFeC were 0.71 and 0.72, respectively, indicating a determinant role of genetic factors for phenotypic variations (Table 1). A significant and positive correlation was found between GZnC and GFeC based on BLUE values (r = 0.45, p < 0.0001). GFeC was positively correlated with PH (r = 0.35, p < 0.0001), whereas GZnC was not (Supplementary Figure 4). TKW was positively correlated with GA and negatively correlated with PH, and appeared not to be correlated with either GZnC or GFeC.

The significant associations were identified between SNPs and GZnC and GFeC using MLM, FarmCPU, and MLMM based on the BLUE values (Figure 1) and the values in each environment (Supplementary Figures 5A, 6A). The Q–Q plots showed that the observed values of p were close to the expected distributions indicating the proper control of false positive in GWAS (Supplementary Figures 5B, 6B). Among the 2,214 and 1,340 significant MTAs for GZnC and GFeC, respectively, 154 and 72 corresponding to 25 and 16 different loci were identified in at least three environments. These 41 loci were considered as stable QTL and are summarized in Table 2. In total, 9, 11, and 8 QTL were only detected by MLM, MLMM, and FarmCPU, respectively; five QTL were identified by two models; eight were simultaneously detected through all the three models (Supplementary Figure 7), indicating the complementarity and reliability of the three models we employed.

For GZnC, 25 stable loci located on chromosomes 1A (2), 2A, 3A (4), 3B (2), 5A, 5D, 6A (2), 6B, 6D (3), 7A, 7B (4), and 7D (3) explained the phenotypic variation (R²-value) ranging from 7.73 to 13.57% (Table 2). Among these, six loci on chromosomal arms 3AS (AX_111528452, 19.89 Mb), 6AS (AX_111556928, 30.88 Mb), 7BL (Tdurum_contig65979_289, 539.22 Mb; BS00022045_51, 626.07 Mb; AX_11046452, 687.84 Mb), and 7DL (AX_108866365, 605.17 Mb) were the most stable and identified in all environments. The most significant marker was AX_111528452 on 3AS with a p at 1.5 × 10^–5 and R²-value at 13.57%. For GFeC, we detected 16 stable QTL on chromosomes 1A, 1B (5), 5A, 5B (4), 7A, 7B, and 7D (3) with R²-values ranging from 7.56 to 14.49% (Table 2). One stable locus on the long arm of chromosome 7D (AX_95151824, 614.54 Mb) was detected in all five environments comprising BLUE values. The most significant QTL was identified on 1BS (AX_111633663) with a p at 1.1 × 10^–5 in the genomic region of 26.17–26.37 Mb.

Four pleiotropic loci were identified by comparing their physical positions of stable MTAs for GZnC, GFeC, and PH (Table 2 and Supplementary Table 3). One locus on 7AL (AX_94741862, 706.91 Mb) was simultaneously associated with GZnC and GFeC while independent of PH and TKW, indicating its potential value in breeding for high grain zinc and iron. Two loci on 5AL (AX_109311262, 495.92 Mb) and 7BL (Tdurum_contig61856_900, 706.45 Mb) significantly increased GFeC and PH, suggesting that they might be useful when plant height is not an issue. The fourth pleiotropic locus located on 3AS (AX_111528452, 20.33 Mb) increased the GZnC but decreased the TKW. All the other 37 MTAs identified have no obvious pleiotropic effects and can be easily used in breeding for genetic improvement of either GZnC or GFeC (Table 2).

Twenty-eight promising candidate genes, located in or adjacent to the physical intervals of the QTL, were identified to be potentially involved in Zn/Fe uptake, transport, storage, and regulations (Table 3). They showed polymorphisms in 641 cultivars which were re-sequenced (data now shown). Among them, 17 genes were located within 2 Mb proximity of the representative markers, and nine genes were located less than 1 Mb from the marker (Table 3). The remaining 11 genes, mainly cloned wheat genes or wheat orthologs of known genes in model plants (7 out of 11), were likely to be causal genes for Zn/Fe homeostasis.

For Zn/Fe uptake, generally proton ATPases, nicotianamine synthases, and phenolic compound transporter played key roles. At least eight candidate genes were identified, and two of them, TraesCS3A01G499300 and TraesCS7D01G098100, were wheat orthologs of rice genes OsPEZ2 (Bashir et al., 2011; Ishimaru et al., 2011) and OsIRT2 (Ishimaru et al., 2006) with known functions for Zn/Fe homeostasis. The expression profiles showed that all the above genes expressed at high levels in roots, further supporting their potential roles in Zn/Fe acquisition (Supplementary Figure 8). For Zn/Fe translocation, key transporters, such as citrate efflux transporter, yellow stripe-like (YSL) transporter, zinc-regulated transporter, and iron-regulated transporter-like protein (ZIP) transporter, were predominantly responsible for this process, and 11 putative genes were identified accordingly. Four genes were wheat orthologs of YSL genes (OsYSL1, OsYSL3, OsYSL7, and OsYSL8) in rice (Waters et al., 2006; Chu et al., 2010; Grillet and Schmidt, 2019), one was cloned gene (TaCNR5) in wheat (Qiao et al., 2019) related to zinc homeostasis, and the remaining six were possible wheat specific or yet to be discovered in model plants. These genes showed diverse expression patterns in different tissues and ages (Supplementary Figure 8). None or very little expression was detected in wheat grain, which partially supported their Zn/Fe transport roles in vegetative tissues (Kobayashi and Nishizawa, 2012).

The storage of metal into vacuoles in seeds was mainly mediated by vacuolar iron transporter (VIT) family members (Zhang et al., 2012), and the efflux of metal from the vacuolar to cytosol was mainly controlled by NRAMP family (Segond et al., 2009), which made TraesCS5D01G209900 (TaVIT2-5D) and TraesCS7B01G364800 (ortholog of AtNRAMP1) promising candidate genes responsible for the Zn/Fe storage in seeds (Segond et al., 2009; Connorton et al., 2017). TraesCS7B01G429600 encoding a zinc binding protein showed high expression in wheat grain and was also considered as a promising candidate gene involved in the Zn/Fe storage (Supplementary Figure 8). Additionally, six regulator genes, such as Vrn-A1, mainly known as transcriptional factors (TFs) regulating the plant development, were identified. TraesCS7B01G299200 and TraesCS7D01G392800, the bZIP TFs, TraesCS7D01G113100, the bHLH TF, TraesCS6A01G0517000, the NAC domain protein, and TraesCS1A01G425100, the wheat ortholog of jasmonic acid-amido synthetase (JAR1), were all putatively related to Zn/Fe homeostasis in wheat. Their spatio-temporal expression profiles are provided in Supplementary Figure 8.

Among all the MTAs, we chose the major QTL on 3AL for GZnC (AX_109875082, 724.58 Mb) with high R²-value (13.53%) and the pleiotropic locus on 7AL for GZnC and GFeC (AX_94741862, 706.91 Mb) to develop high-throughput KASP markers (Table 2 and Supplementary Table 4). Both QTL were independent of TKW and PH (Table 2). Their representative SNPs were successfully converted to KASP assays after genotyping the association panel and the same genotypic data were obtained from the wheat 660K SNP array (data not shown). The marker information was presented in Supplementary Table 4. The two KASP markers were further used to test the ZM175/LX987 bi-parental population. For K_AX_109875082 on 3AL, a t-test showed a significant difference of GZnC between lines carrying favorable allele and lines with opposite allele across environments (Figure 2B). No polymorphism was detected for K_AX_94741862 on 7AL in this population. Those markers closely linked to these QTL could be further detected in the ZM175/LX987 bi-parental population.

For the 3AL QTL, the most significant markers were located in the physical interval of 724,305,362–724,582,163 bp, where nine genes were annotated (Figures 2A,C). Sequence similarity analysis showed that TraesCS3A01G499300, annotated as a multidrug and toxic compound extrusion (MATE) family gene responsible for protein detoxification, had high amino acid sequence identity (87%) with OsPEZ2 (Tong et al., 2020). The gene OsPEZ2 played an important role in metal uptake and translocation in rice and pez2 mutant showed significantly reduced Fe concentration in root tips and xylem sap (Bashir et al., 2011; Ishimaru et al., 2011). TraesCS3A01G499300 is 5,562 bp long and consists of seven exons in coding region (Figure 2C). Further sequence analysis revealed that TraesCS3A01G499300 had abundant SNPs in its flanking and coding regions between ZM175 and LX987 (Supplementary Table 5), especially four missense SNPs within exons (Figure 2C), indicating TraesCS3A01G499300 is likely to be the causal gene underlying 3AL QTL. Gene expression data indicate that the gene is constitutively expressed and is uniformly highly expressed in wheat root throughout developmental stages, which further supports its candidacy (Figure 2C).

In previous GWAS studies on wheat Zn/Fe, large gaps in the genetic map were generally found due to the low marker number and density as well as relatively long LD decays. In this study, GWAS was conducted using high-density markers, such as 373,106 SNPs, making it possible to change the LD pattern, thereby resulting in efficient identification of MTAs in many low-recombination regions on wheat chromosomes (Kim and Yoo, 2016). Forty-one loci for GZnC/GFeC were detected, and the number is much higher than the QTL detected in previous literature (Gupta et al., 2021). Out of the 41 representative significant MTAs, only five markers came from the 90K SNP array, and 36 were from the 660K array, indicating that high-resolution GWAS was achieved by high density markers (Fu et al., 2020; Pang et al., 2020).

In terms of population size, a relatively small collection was used in this study. Nevertheless, the 166 core accessions were selected from over 400 elite cultivars and represented much diverse wheat accessions in China’s major wheat growing regions. The same population has been used to conduct GWAS on many yield or quality-related traits and performed highly informative in our previous studies (Liu et al., 2017; Zhai et al., 2018, 2020; Li et al., 2019; Fu et al., 2020). MLM using Q + K as covariant might be over-fitted possibly due to the strict control of population structure and kinship, leading to false negatives (Korte and Farlow, 2013). The multiple testing corrections tended to be too stringent to detect sufficient both major and minor MTAs. Whereas, the model FarmCPU as a complement, could avoid over-corrected population structure to some extent by utilizing fixed and random effect models iteratively (Kusmec and Schnable, 2018). Using three models together, reliable loci for GZnC and GFeC detected were almost doubled than using an MLM alone in the present research.

To date, hundreds of Fe/Zn-related QTL have been mapped on all wheat chromosomes except 6D by bi-parental linkage mapping and GWAS (Tong et al., 2020). For the 25 loci associated with GZnC in this study, eight have similar physical positions with known Zn-related QTL on 1AL (2), 3AL (2), 6BS, 7AL, and 7BL (2), respectively (Table 2 and Figure 3). As for GFeC, half of the 16 loci were identified to coincide with the documented QTL on 1BS, 5AL, 5BL (3), 7AL, 7BL, and 7DL, respectively. Major overlapping regions were observed on 5BL, where previous literatures reported a large number of Fe-related QTL by bi-parental mapping, confirming its importance for the gene discovery (Tong et al., 2020).

For the other 25 loci, comprising 17 loci for GZnC and 8 for GFeC, they were located in the chromosomal regions that were different from the QTL previously reported and are probably new loci (Table 2 and Figure 3). For example, three QTL on 6D (AX_108846745, 16.81 Mb; AX_95220141, 27.52 Mb; AX_110431664, 357.01 Mb) are highly likely to be novel since no Zn/Fe-related QTL has been previously mapped on this chromosome (Tong et al., 2020). Similarly, on chromosome 7D, there were only two Fe related QTL reported at 13.9 Mb and 616.6 Mb, so the three loci for GZnC (AX_110717434, 203.16 Mb; GENE_3452_1116, 506.11 Mb; and AX_108866365, 605.17 Mb) and two loci for GFeC (AX_108920250, 54.99 and AX_111359934, 69.31 Mb) were likely to be new. The new loci identified in this research have significantly enriched our understanding of the genetic basis of the complicated Zn and Fe traits in wheat.

To date, numerous genes involved in Zn and Fe homeostasis in model plants have been extensively studied (Tong et al., 2020; Gupta et al., 2021). Wheat gene annotations from the IWGSC RefSeq v.1.0 provided a useful tool to further investigate the candidate genes located in or adjacent to the QTL identified. In the current study, 28 candidate genes were identified for Zn/Fe homeostasis, such as uptake, transport, storage, and regulations (Figure 4). The spatio-temporal expression profiles were highly consistent with their putative functions, indicating their possible roles in wheat Zn/Fe homeostasis. Of them, 11 genes were putative wheat orthologs of known Arabidopsis and rice genes related to Zn/Fe homeostasis, partially supporting the authenticity of the detected loci.

TraesCS3A01G499300, the candidate gene for 3AL QTL, was an ortholog of OsPEZ2, an efflux transporter for phenolics that served as chelators to facilitate metal uptake by root in rice (Bashir et al., 2011; Ishimaru et al., 2011). A very recent study demonstrated the pivotal roles of MATE family genes of protein detoxification for metal uptake and transport in rice roots and stems, further supporting the potential role of TraesCS3A01G499300 on Zn/Fe homeostasis (Ren et al., 2021). Abundant sequence variations of the gene between ZM175 and LX987 were observed and provided strong evidence that it is possibly the causal gene and deserves further studies. In addition, we successfully developed KASP markers for the 7AL QTL but unfortunately it returned monomorphic in ZM175/LX987 RILs. However, it was interesting to find that the candidate gene TraesCS7A01G527900 (706.91 Mb) underlying 7AL QTL was homologous with TraesCS7B01G444600 (708.11 Mb) on 7BL and TraesCS7D01G515800 (615.91 Mb) on 7DL for the three QTL identified in this study. The three candidate genes annotated as magnesium transporters were exactly the three homeologs for the same gene in each sub-genome. These loci are potentially important and we are now developing KASP markers for 7BL and 7DL QTL for validation. If any loci are validated in ZM175/LX987 or other mapping populations, the roles of all three loci in metal homeostasis can be presumed. The candidate genes identified in this study provided genetic bases for further elucidating the mechanisms of Zn and Fe homeostasis in wheat.

Since most common wheat had suboptimal GZnC and GFeC, it is imperative to identify Zn/Fe-related QTL or genes as many as possible and to re-introduce or to pyramid them into current wheat gene pools. In this study, 12 accessions with stable high GZnC or GFeC were identified with BLUE values over 50 mg/kg (Supplementary Figure 3). The elite germplasm, such as Xiaoyan54 and Xinong1376 are currently cultivars in wheat production and can be used immediately as donor parents for wheat biofortification.

Although no significant correlation was observed between TKW and GZnC/GFeC in the present study, one pleiotropic locus on 3AS (AX_111528452, 20.33 Mb) was identified to increase GZnC but decrease TKW in all four environments. This QTL may be useful in certain crosses for biofortification when TKW is high enough. In addition, when plant height does not significantly affect the overall performance, pleiotropic loci on 5AL (AX_109311262, 495.92 Mb) and 7BL (Tdurum_contig61856_900, 706.45 Mb) could exert effects on Zn/Fe enrichment. It should be noted that favorable allele frequencies for the above three representative markers were 0.31, 0.09, and 0.36, respectively, indicating quite low frequency in our association panel (Supplementary Table 3). Considering only very few GZnC/GFeC loci had negative pleiotropic effects on agronomic traits such as TKW and PH, we speculated that there were still high chances for the improvement of GZnC/GFeC without yield penalty (Velu et al., 2018).

The fast release of wheat reference genome sequences and pan-genomes will undoubtedly speed up the process of marker development and gene discovery for wheat biofortification. Cloning important genes involved in Zn/Fe uptake, transport, storage, and regulations and pyramiding favorable alleles are promising avenues to increase the GZnC and GFeC in wheat cultivars. The 3AL locus was validated using a bi-parental population and its candidate gene was proposed as an uptake related gene, highlighting its potential use in wheat breeding. High-throughput and breeder-friendly KASP markers will pave the way for MAS in breeding and accelerate the release of biofortified wheat.