Genetic dissection of zinc, iron, copper, manganese and phosphorus in wheat (Triticum aestivum L.) grain and rachis at two developmental stages

The development of high-yielding wheat genotypes containing micronutrient-dense grains are the main priorities of biofortification programs. At the International Maize and Wheat Improvement Center, breeders have successfully crossed high zinc progenitors including synthetic hexaploid wheat, T. dicoccum, T. spelta and landraces to generate high-zinc varieties. In this study, we report a genome-wide association using a wheat diversity panel to dissect the genetics controlling zinc, iron, copper, manganese and phosphorus concentrations in the grain and rachis during grain development and at physiological maturity. Significant marker-trait associations (MTAs) were identified for each nutrient using multi-locus mixed model methodologies. For mature grain, markers that showed significant pleiotropic effects were found on chromosomes 1A, 3B and 5B, of which those on chromosome 5B at ∼95.5 cM were consistent over two growing seasons. Co-located MTAs were identified for the nutrient concentrations in developing grain, rachis and mature grain on multiple chromosomes. The identified genomic regions included putative candidate genes involved in metal uptake and transport and storage protein processing. These findings add to our understanding of the genetics of the five important nutrients in wheat grain and provide information on genetic markers for selecting high micronutrient genotypes.


Introduction
Cereal grains including wheat are inherently low in essential micronutrients, thus relying on cereals as a major staple can lead to micronutrient malnutrition; a problem affecting over three billion people worldwide [1,2]. The micronutrients most frequently lacking in the human diet of developing country populations are iron (Fe) and zinc (Zn) and these elements are critical for a wide range of metabolic processes [3]. Fe is important for oxygen transport and haemoglobin formation, whereas Zn plays a central role in growth, development and in the immune system. Deficiency in these micronutrients could result in poor health, increased morbidity and mortality rates, and lower worker productivity [4,5]. The impact is more serious in women of reproductive age (especially pregnant women) and children under the age of five, due to their increased micronutrient requirements [6,7]. Though copper (Cu) or manganese (Mn) deficiency in humans is not a widespread global issue, these trace elements also play pivotal roles in growth and development [5]. Both minerals function as cofactors or components of several key enzyme systems. Manganese deficiency in animals is associated with changes in carbohydrate and lipid metabolism and impaired growth and reproductive function [8]. Copper deficiency results in defects in connective tissue that lead to vascular and skeletal problems and anaemia related to defective iron metabolism [9]. The negative consequences of micronutrient deficiencies for health, productivity and mental development are devastating and often lifelong and these problems require unified global efforts to alleviate [4,10].
An intervention method to combat micronutrient deficiencies is to improve the nutritional values of crop plants and is known as biofortification [11]. Two main biofortification strategies include (i) genetic biofortification, which is based on traditional plant breeding and/ or genetic engineering for improved nutrient concentrations and (ii) agronomic biofortification, which is based on optimised fertilizer applications [12,13]. Genetic biofortification requires minimal ongoing investment once the crop lines with higher nutrient densities are developed, thus this approach is highly cost-effective and more useful in developing countries.
Breeding crops for improved micronutrient levels requires adequate genetic variation for exploitation. Modern wheat cultivars are generally not a good source of genes for high micronutrient levels, likely because modern cereal breeding has focused primarily on improving yield, disease resistance and starch quality rather than grain micronutrient contents [14]. There was a significant trend of decreasing micronutrient concentration with the date of variety release and the current cereal cultivars typically had low grain mineral concentrations [15,16]. In contrast, the wild and primitive wheats represent a promising genetic resource. The collections of Triticum boeoticum, Triticum monococcum, Triticum dicoccoides, Aegilops tauschii, Aegilops speltoides and synthetic hexaploid lines were found to have impressive genetic variation and superior grain Fe and Zn concentrations in comparison to cultivated wheat varieties [17][18][19]. This large genetic diversity has great potential for exploitation to breed for increasing Zn, Fe and possibly other mineral levels in biofortification programs.
However, little information is available about the genetic controls and molecular physiological mechanisms contributing to high accumulation of grain micronutrient levels of the exotic genetic materials to date [20,21]. Until recently, Velu et al. [22] published a report on genetic dissection of wheat grain Zn concentration in a diversity panel established at the International Maize and Wheat Improvement Center (CIMMYT). The wild relatives and landrace genomes known to have elevated grain micronutrient levels include T. dicoccon, T. spelta and Ae. Squarrosa. Segments of their genomes were introgressed into modern varieties through the development of synthetics and backcrossing to create a rich genetic diversity within this panel [22,23]. The authors reported 39 genetic markers (on nine genetic linkages) associated with the variation in grain Zn concentration and suggested seven potential candidate genes including zinc finger motif of transcription-factors and metal-ion binding. In continuation of the study by Velu et al. [22], we investigated the genetic basis of five essential nutrients including Zn, Fe, Cu, Mn and phosphorus (P) in wheat grain and rachis at an early stage of grain development and at physiological maturity.
Grain mineral supplies come from two sources, the direct uptake from the soil and/or the remobilization of stored minerals in vegetative tissues during grain filling [24]. In wheat, studies have shown that the majority of grain micronutrients came from the reserves in the vegetative tissues, with more than 70 % of the Zn reserves in the vegetative tissue of wheat plants remobilized during grain filling [25,26]. Up to 77 % and 62 % of shoot Fe and Cu was transported into grain [27]. In rice, the translocation of Zn from rachis to grain was most critical for grain Zn accumulation [28]. It is apparent that the available pool of minerals in vegetative tissue, especially the rachis during the reproductive growth stage is of critical relevance for their remobilization into grain. Our study aimed to elucidate the relationships between the key nutrients in the grain and between the grain and rachis at different development stages and to dissect the genetics underlying mechanisms that control mineral accumulation in both grain and vegetative tissues. The results would improve our understanding of the mechanisms that control micronutrient loading into the grain. The resulting marker-trait associations (MTAs) and potential candidate genes would assist breeding for increasing mineral concentrations in wheat varieties.

Plant materials and field trial
The HarvestPlus Association Mapping (AM) panel consisted of 330 wheat lines from CIMMYT's biofortification breeding program. This panel contained large genetic diversity and could be divided into five sub-populations: (1) progenies derived from landraces; (2) Triticum durum-based synthetic hexaploid derivatives; (3) Triticum dicocconbased synthetic hexaploid derivatives; (4) Triticum spelta derivatives; and (5) CIMMYT's pre-breeding derivatives of diverse progenitors including T. polonicum [22]. The panel was grown in a randomized complete block design in two replications with the plot size of 2 m 2 at CIMMYT's experimental station in Ciudad Obregon, Sonora, Mexico (27°24ˈN, 109°56ˈW)) during two successive crop seasons (2014-15 and 2015-16). During the growing seasons (April-November) in 2014-15 and 2015-16, the minimum temperature ranged from 1.1-7.6°C and 4.0-8.5°C and the maximum temperature ranged from 28.9-34.5°C and 30.6-38.5°C, respectively. The trials were irrigated five times throughout the crop cycle and fertilized at a rate of 200-50 (N-P), of which 50-50 was applied in pre-sowing and 150-00 at tillering stage. Diseases and pests were controlled chemically, whereas weeds were controlled manually and chemically according to CIMMYT's standard protocols. To minimize Zn heterogeneity or gradient, the soil of the field trials was previously enriched with 25 kg ha −1 of ZnSO 4 .7H 2 O over three crop cycles. Soil analysis of the experimental area showed an average Zn concentration of 1.2 ppm at soil depth of 0-30 cm, and 0.86 ppm at a soil depth of 30-60 cm. The average Fe concentration in the soil was 5.0 and 6.1 ppm, at 0-30 and 30-60 cm soil depth, respectively. Field observations were recorded for plant height (PH), days to heading (DH) and days to physiological maturity (DM), grain-filling period (GFP) and thousand kernel weight (TKW) according to CIMMYT standard protocols [29].

Sampling
Wheat head samples were collected at two developmental stages; namely immature and mature. For immature samples, five heads per genotype in each plot were collected at 10 days after anthesis in both growing seasons. For mature samples, five heads per genotype in each plot were collected at physiological maturity in 2014-15. In 2015-16, the plots were harvested at physiological maturity and 20 g of mature grain were sub-sampled from the bulk. All samples were oven dried at 40°C for 4 days and each head was analysed separately as one replicate.

Nutrient analysis
Dried immature and mature wheat heads were dissected to separate the grain, rachis, peduncle and floret. The immature grain and rachis and mature grain samples were subsequently subjected to nutrient analysis. Approximately 0.3 g of each sample (oven dried at 80°C for 4 h) was digested with a closed tube acid digestion as described in [30]. Mineral concentrations were determined using inductively coupled plasma mass spectrometry (ICP-MS 7500x; Agilent, Santa Clara, CA) following the method described in [31]. A blank and a certified reference material (CRM; NIST 1567a wheat flour) were included in each digestion batch for quality assurance. Mineral concentrations were expressed on a dry weight basis. Contamination with soil was monitored by analysis for aluminium and titanium.

Statistical analysis
Genstat 18 th software was employed to carry out statistical analyses [32]. The Pearson's correlation coefficient was used to evaluate the relationships among the measured parameters of data. The significance of the correlations was determined with two-sided test of the correlations against 0 at three levels of probability; 0.001, 0.01 and 0.05. The means of different seasons were compared using one-way ANOVA at a 0.05 level of probability. The frequency distributions of grain mineral concentrations and TKW were demonstrated using Histogram. Shapiro-Wilk test was used to assess whether the traits were normally distributed.

Genotypic data, population structure and linkage disequilibrium
DNA was extracted according to standard CTAB procedures and genotyped using the Illumina iSelect 90 K Infinitum SNP array developed in wheat [33]. The array revealed a total of 28,074 SNP markers, which we coded as 0 and 1 for homozygote and 0.5 for heterozygote allele scores for data analysis. Markers with less than 0.05 allele frequency, greater than 20 % missing data and greater than 10 % heterozygosity were removed giving a remaining total of 17, 900 SNP markers for subsequent analysis.
A kinship matrix was generated using Centered-Identity by State (IBS) and the basic population structure was analysed using Principal Component Analysis (PCA) in TASSEL 5.2. Linkage disequilibrium (LDs) between SNPs were estimated as the squared correlation coefficient (R 2 ) of alleles within a 10 cM window (Supp. Table S1, Figs. S1-S3).
The two single-locus models were performed using TASSEL 5.2. The GLM corrected only the population structure and the MLM corrected both population structure and kinship relationship among individuals. Significant levels of marker-trait association were set at an adjusted Pvalue of 0.05/n e where n e was the effective number of independent markers (modified Bonferroni correction for multiple tests). The n e number was determined using snpgdsLDpruning function in SNPRelate package of R software to prune SNPs which were in high LD with each other (LD threshold = 0.9) [38]. Manhattan plots and QQ plots were created using TASSEL5.2.
The multi-locus mixed models, which involved two steps were performed using the mrMLM.GUI software. In the first step, a singlelocus GWAS method was applied to scan the entire genome, and putative markers were selected with a low criterion of significance test. In the second step, a multiple-locus GWAS method was implemented for the markers that had passed the initial screening. Statistical tests and marker effect estimation were then based on the multiple-locus model. All parameters were set at default values. In the first step, the critical Pvalues were set to 0.01 for all models except for FASTmrEMMA, where it was set to 0.005. In the second step, the significant threshold was set to P value ≤ 2 × 10 −4 (-log(P value ) = 3.7) for all models [35]. Model comparison was performed using the Mean Square of the Difference (MSD), based on observed P values and expected P values, and Q-Q plot. (Supp. Figs. S3 and S4).
GWAS analysis was conducted for each growth season separately. Phenology traits (DH and DM), plant height and thousand kernel were used as covariates for MLMM where they had significant correlations with the traits of interests. The mean of the squared correlation coefficients (R 2 ) between observed phenotypic measurements and their predicted values was taken as an estimate of the proportion of phenotypic variance explained by the marker. The phenotypic-effect value of each allelic variation was calculated by the phenotypic data over the accessions with each type.

Data mining from the wheat genome sequencing
All markers associated with significant QTL were blasted on the Ensemble Plants website (http://plants. ensembl.org/ Triticum_aestivum/Tools/Blast?db = core) to find the Gene identifiers Table 1 Phenotypic variation in phenology traits, thousand kernel weight (TKW) and the concentration (mg kg −1 ) and content (μg seed −1 ) of micro-and macro-nutrients between two growing seasons. Mean values labelled with different letters (a, b) indicate the differences between two seasons are statistically significant (ANOVA, P < 0.05). r values are the Pearson's correlation coefficients between the two growth seasons (p ≤ 0.001, two-sided test).

Phenotypic variation between two growth seasons
Twenty-five traits were measured for the AM panel and the maximum, minimum, mean values and standard deviations (SD.) of the traits are presented in Table 1. All these traits showed continuous variation and followed approximate normal distributions. The frequency distribution of Zn, Fe, Cu, Mn, P concentrations in mature grain as well as thousand kernel weight (TKW) are shown in Fig. 1.
The differences between two seasons for DH, DM, GFP, PH and TKW were statistically significant (one-way ANOVA, P < 0.05). In 2015-16, the AM panel tended to have longer DH, DM and GFP, with shorter plants and smaller kernels than the previous season ( Table 1).
The mature grain results showed that the Zn concentration was significantly higher in 2015-16, while Fe, Cu and P concentrations were higher in 2014-15 whereas Mn concentration did not change between the two seasons. In immature grain and rachis, Zn, Cu and P concentrations were significantly higher in 2014-15 while Fe and Mn concentrations were higher in 2015-2016 (Table 1).

Correlations between nutrient concentrations
Correlation analysis was performed for phenotypic traits in each season. The resulting correlation coefficients (r) and their significance levels are shown in Table 2A (2014-15) and 2B (2015-16).
In mature grain, all five nutrient concentrations were positively correlated with each other (r = 0.28-0.78; P < 0.001). These correlations were relatively consistent over the two seasons for all nutrients.

Table 2
Pearson's correlation coefficients (r) between mineral concentrations in mature grain (MG), immature grain (IG) and immature rachis (IR) and days to heading (DH), days to maturity (DM), grain-filling period (GFP), plant height (PH) and thousand kernel weight (TKW). The correlations were calculated from 330 lines of the AM panel grown in Obregon, Mexico in two consecutive seasons (A) 2014-15 and (B) 2015-16. The significance thresholds for r values determined by two-sided test are: ***P < 0.001, **P < 0.01 and *P < 0.05.
Zn concentrations in mature grain, immature grain and immature rachis were positively correlated to each other in in both seasons (r = 0.18-0.35; P < 0.01). Similar correlations were observed for Mn and Cu concentrations. In the case of Fe and P, the correlations differed considerably between seasons. Correlation between Fe concentrations in mature grain and immature rachis was significantly only in 2014-15 while between mature grain and immature grain was significantly only in 2015-16. For P concentration, significant correlations between immature rachis and mature grain (r = 0.29; P < 0.001) and between immature rachis and immature grain (r = 0.60; P < 0.001) were found in 2015-16 but not in the previous year.

Effects of phenology traits
The relationships between nutrient concentrations in mature grain and days to heading (DH), days to maturity (DM), grain filling period (GFP), plant height (PH) and thousand kernel weight (TKW) were investigated. TKW was the only trait that had consistent positive correlations to all nutrient concentrations over two seasons (r = 0.12-0.25; P < 0.01). Days to heading (DH) and days to maturity (DM) were both negatively correlated to Zn, Fe, Cu and Mn concentrations in 2014-15 (r= -0.12 to -0.22; P < 0.05) while no significant relationship with any nutritional traits in 2015-16. Plant height was positively correlated to Zn, Cu and Mn concentrations only in 2015-16 (r = 0.15 to 0.18; P < 0.05). Grain-filling period (GFP) had no significant correlations to any nutrient in both seasons.

Marker-trait association
In total, 279, 379 and 481 significant MTAs were detected for five nutrient concentrations in mature grain, immature grain and immature rachis, respectively by the three MLMM models simultaneously (P < 2 × 10 −4 ). In contrast, the single-locus models detected only a couple of MTAs that qualified the critical threshold of P < 1.25 × 10 -5 (modified Bonferroni correction for multiple tests). The summaries of the genomic regions and SNPs associated with each trait detected by the MLMM models simultaneously are presented in Tables 3 and 4 and  Supp. Table S2.

Mature grain
In mature grain, 72 MTAs were significantly associated with Zn concentration, 65 with Fe, 30 with Cu, 47 with Mn and 46 with P (Table 3, Supp. Table S2). The MTAs were identified on 15 linkage groups and the highest numbers were on chr 5B (59), 1A (47) and 2A (32).
Zn concentration was most significantly associated with the markers on chr 1A and 2A (LOD scores > 5) in 2014-15 ( Table 3). Each of these MTAs could explain for approximately 3.7-4.2 % of the phenotypic variation and the estimated additive effects were around 4.1-4.8 mg/ kg. The highest additive effects estimated at ∼ 4.6-5.5 mg/kg came from the markers on chr 5B (95 cM). In 2015-16, the most significant MTAs were detected on 1D, 2B and 5B (LOD > 5). Each of these markers accounted for ∼4.2-5.2 % of the phenotypic variation and the higher alleles were estimated to increase Zn concentration by 4.6-5.9 mg/kg. There were five stable MTAs on chr 5B at ∼95 cM, detected in both seasons (Supp . Table S2).
For Fe concentration, the most significant MTAs were identified on chr 5A and 5B in 2014-15 and on chr 7B in 2015-16 (LOD > 4.9) ( Table 3) Table S2). The highest additive effects were ∼0.6 mg/kg in each season.
The most significant MTAs for Mn concentration were detected on chr 6D (87 cM) for both seasons (Table 3). Those markers were estimated to have additive effects of ∼ 4.9-5.4 mg/kg. Apart from those, five stable MTAs for Mn were also found on chr 5B (95 cM).
For P concentration, five stable MTAs were detected on chr 5B (95 cM). Those MTAs had the highest LOD scores in both seasons (5.

Immature grain and rachis
In immature grain, stable MTAs were identified for Zn concentration on chr 7B (75 cM); Cu on 2A (74 cM) and Mn on 5B (68 cM) and 7B (67−68 cM) ( Table 4). In immature rachis, stable MTAs were detected for Zn concentration on chr 5A; Fe on 2A, 2B and 2D; Cu on 5A and 6D and P on 2B. No stable MTAs were detected for Fe or P concentration in immature grain or Mn in immature rachis.
The genomic region on chr 7B at ∼67-77 cM had the highest numbers of MTAs for both immature rachis and immature grain; 56 and 66, respectively. In immature rachis, the genomic region was associated with four nutrient concentrations (Zn, Fe, Cu and Mn) in 2014-15. In immature grain, it was associated with Zn and Mn concentrations in both seasons and Cu concentration in 2014-15.

Co-located MTAs
In mature grain, multiple genomic regions harboured markers that were associated with more than one nutrient concentrations (Table 3, Supp. Table S2). The most significant one was on chr 5B at ∼95 cM, which was associated with all five nutrients in both seasons (except for Fe in 2015-16). Two genomic regions that were associated with three nutrients were on chr 1A at ∼113 cM (Zn, Fe and Mn) and 3B at ∼11-14 cM (Fe, Cu and Mn). Furthermore, Zn and Fe concentrations also had co-located MTAs on chr 2A (106-108 cM); Fe and Cu concentrations on chr 7B (102 cM). Except for those on chr 5B, the other co-located MTAs were identified in only one season.
Nutrient concentrations in mature grain shared multiple co-located MTAs with those in immature grain and rachis and (Tables 3 & 4). The region on chr 2A at ∼109 cM was associated with Zn concentration in all sample types. Furthermore, between mature and immature grain, Fe concentration shared co-located MTAs on chr 1B (99 cM), 3B (8-14 cM) and 5B (95 cM), Cu on chr 2B (69 cM), Mn on 1A (113 cM) and P on 2A (127 cM). Between mature grain and immature rachis, Fe concentrations had co-located MTAs on chr 1A (113 cM) and 2A (104-108 cM), P on 1D (67 cM) while none was found for Cu or Mn.

Days to heading, days to maturity, plant height and thousand kernel weight
Twelve significant MTAs were detected for DH; 17 for DM and seven for TKW (Supp . Table S3). There were no significant MTAs detected for plant height and grain-filling period. The MTAs were compared to those found for nutrient concentrations and genomic regions harbouring colocated MTAs were identified. The region at 68-70 cM on chr 7B associated with DM was in close proximity to the MTAs detected for Mn concentration (2014-15). This genomic region also had the largest number of MTAs found in developing grain (Zn and Mn concentrations) and immature rachis (Zn, Fe, Cu and Mn concentrations). The genomic regions on chr 1B (158 cM) and 5B (144 cM) were associated with DH and DM and also with Fe and P concentrations in immature grain.

Genes encoding observed MTAs
The identified markers their physical positions were used to locate genes that were either linked to or in close proximity (within 1 cM) by using the annotated wheat reference sequence (RefSeq V1.0). Table 5 lists 116 potential candidate genes linked to key MTAs on chr 1A, 2A, 3B and 5B. The genes were grouped based on their putative functionality. The largest groups were: (i) transporter families including Zinc transporter, Yellow stripe-like (YSL) 1, 9 and 12 transporter, Vacuolar iron transporter (VIP), copper, potassium, sulphate, phosphate and sugar transporters, (ii) protein processing enzymes (proteinases, proteases and peptidases), (iii) storage proteins (late embryogenesis abundant protein, dehydrin, lipid transfer protein), (iv) metalloprotein and (v) antioxidant enzymes (Zn-Cu dismutases, Glutathione-S-transferase, Glutaredoxin). There were also genes involved in grain development (sucrose synthases, starch synthases), plant development and stress responses.

Multi-locus mixed models
The implementation of GWAS is useful to identify genomic regions and/or genes associated with traits of interest and to utilize the information of associated markers in breeding programs to efficiently incorporate particular loci in elite germplasm [40]. To date, there have been many studies on genetic bases of complicated traits in wheat and Table 4 Summary of the genomic regions significantly associated with Zn, Fe, Cu, Mn and P concentrations (mg kg −1 ) in immature grain (A) and immature rachis (B) harvested at 10 days after anthesis from the AM panel in two seasons using MLMM models. Chr, chromosome; Pos, centimorgan position. Genomic regions detected in both seasons were shaded. Serine/threonine-protein phosphatase TraesCS1A01G374200 Lipid transfer protein TraesCS1A01G374800 Zinc finger protein-like TraesCS1A01G375300 Ubiquitin carboxyl-terminal hydrolase 12 TraesCS1A01G377400 Sucrose synthase 6 TraesCS1A01G377500 Sucrose synthase 6 TraesCS1A01G379900 Thioredoxin superfamily protein TraesCS1A01G380100 sucrose-proton symporter 9 TraesCS1A01G380300 E3 ubiquitin-protein ligase SINA-like 10 TraesCS1A01G382600 Protease inhibitor/seed storage/lipid transfer protein TraesCS1A01G382800 Asparagine synthetase TraesCS1A01G383200 2Fe-2S ferredoxin-type iron-sulfur binding TraesCS1A01G394100 Calcium-binding family protein TraesCS1A01G394400 Structural constituent of cell wall TraesCS1A01G394500 3-deoxy-manno-octulosonate cytidylyltransferase TraesCS1A01G394600 Purple acid phosphatase TraesCS1A01G394800 Calcium-dependent lipid-binding family protein TraesCS1A01G394900 Cytochrome P450, putative TraesCS1A01G395200 Protein DEHYDRATION-INDUCED 19  TraesCS1A01G396200 Xyloglucan endotransglucosylase/hydrolase TraesCS1A01G396300 Cyst nematode resistance protein-like TraesCS1A01G396400 Auxin response factor TraesCS1A01G396200 Xyloglucan endotransglucosylase/hydrolase TraesCS1A01G396900 Glutaredoxin, putative TraesCS1A01G397200 Glutaredoxin, putative TraesCS1A01G398100 Beta purothionin TraesCS1A01G398200 Beta purothionin 2A 103-108 Zn TraesCS2A01G143400
Pos (cM) Trait Gene ID Putative functionality 2A 103-108 Zn TraesCS2A01G287100 Carboxypeptidase Fe TraesCS2A01G287200 Heavy-metal-associated family TraesCS2A01G287300 Chloroplast ATP-dependent Clp protease chaperone TraesCS2A01G293400 1,4-alpha-glucan-branching enzyme TraesCS2A01G293800 Late embryogenesis abundant protein TraesCS2A01G294200 Purple acid phosphatase TraesCS2A01G294700 Glucan endo-1,3-beta-glucosidase TraesCS2A01G295000 Glutaredoxin family protein TraesCS2A01G295700 Senescence regulator TraesCS2A01G297000 Copper transporter family protein TraesCS2A01G297100 Protease inhibitor/seed storage/lipid transfer protein TraesCS2A01G373200 Sugar transporter, putative TraesCS2A01G373300 Alpha 1,4-glycosyltransferase family protein TraesCS2A01G373600 Starch synthase, chloroplastic/amyloplastic TraesCS2A01G376800 Serpin (serine protease inhibitor) TraesCS2A01G377200 Yellow stripe-like transporter 12 TraesCS2A01G377400 Yellow stripe-like transporter 12 TraesCS2A01G377600 Yellow stripe-like transporter 12 TraesCS2A01G386300 NAC domain protein, TraesCS2A01G386600 Vacuolar-processing enzyme TraesCS2A01G386700 Vacuolar-processing enzyme TraesCS2A01G387100 Vacuolar iron transporter TraesCS2A01G390800 Yellow stripe-like transporter 1 TraesCS2A01G391000 Transporter family protein TraesCS2A01G391100 Yellow stripe-like 9 transporter TraesCS2A01G394500 Heavy metal transport/detoxification superfamily TraesCS2A01G395500 RING-H2 zinc finger protein TraesCS2A01G552000 Vacuolar protein sorting-associated protein TraesCS2A01G552200 Pyruvate decarboxylase TraesCS2A01G553300 Yellow most of them used single-locus GWAS models (GLM and/or MLM) [22,[41][42][43][44][45][46][47][48]. These models require multiple testing corrections to control for false positives. The corrections, however, are often too stringent and leading to many important MTAs not being detected, especially in the cases of small-effect MTAs or large errors in phenotypic measurement from field experiments [35,39]. To overcome this issue, MLMM methodologies have been developed and shown to be much more powerful and accurate [35][36][37]39,49,50]. In this study, two single-locus and three multi-locus GWAS models were tested with Zn and Fe concentrations. The GLM had a high false positive rate, and the MLM detected only two MTAs that qualified for the stringent P threshold (P = 0.05/n e where n e was the effective number of independent markers) (Supp. Fig. S3). The multi-locus models demonstrated a clear advantage of higher detection power and accuracy compared to the single-locus models (Supp. Fig. S4). This is most likely because the control of Zn and other metal transport into the grain is a complex process involved a large number of genes and/or gene families with each having a small effect that could not be detected by the single-locus models. The MLMM models, therefore were chosen to analyse all traits in this study. Theoretically, correction for multiple tests is unnecessary in multi-locus GWAS because all the potential genes or loci for complex traits are fitted to a single linear model and their effects are estimated and tested simultaneously. The critical threshold of P value ≤ 2 × 10 −4 was set for declaring significant MTAs. This threshold has been shown to sufficiently control false positives while maintaining high power of MTA detection in multi-locus GWAS [35,39]. To maximise the credibility of the MTAs, only those detected by all three multi-GWAS models were reported.

Variation in nutrient concentrations
Breeding wheat varieties with high micronutrient levels requires source materials that feature adequate genetic variation in micronutrient concentrations. In our study, wide ranges of grain Zn, Fe, Cu and Mn concentrations were observed in both seasons (Table 1, Fig. 1). Grain Zn concentrations ranged from 38 to 67 and 39−76 mg/kg in 2014-15 and 2015-16, respectively. Similarly, Fe, Cu and Mn concentrations also varied by approximately two folds between the maximum and minimum values. The variation in the grain micronutrient concentrations were likely due to the diverse genetic resources preserved in CIMMYT's wheat germplasm with wild relatives introgressed into elite varieties [23,51].

Correlations between nutrients in mature grain and their co-located MTAs
Correlation analysis showed significant positive correlations between mature grain Zn, Fe, Cu, Mn and P concentrations in both growth seasons (r = 0.48-0.78, p < 0.001). Our results are in agreement with the strong correlations between those elements reported in previous studies using various germplasm containing wild emmer wheats, durum wheats or cultivated hexaploid wheats [17,[52][53][54]. These strong correlations likely resulted from common genetic factors controlling the accumulation of the minerals in mature grain. The underlying genetics was further elucidated with our GWAS results, which identified multiple genomic regions harbouring multi-trait MTAs (Table 3). There were three genomic regions that associated with three or more nutrients in mature grain; on chr 1A, 3B and 5B. These multi-trait MTA regions have been reported in previous QTL studies [43,[55][56][57]. However, due to differences in the genetic markers used in each study (e.g DArT, SSR, GBS or Illumina 90k SNP markers), the comparisons between genomic regions were only approximate estimations. The region around 95 cM on chr 5B was found controlling the grain Zn and Fe concentrations in mapping populations developed from diverse genetic backgrounds such as synthetic hexaploid wheat and a T. spelta L. derived line, durum wheat cultivar and a tetraploid wild emmer wheat (T. dicoccon) [55,58]. Water soluble carbohydrate concentration, kernel length and width and TKW were also significantly associated with this genomic  TraesCS3B01G017100  ABC transporter G family member  TraesCS3B01G018000  E3 ubiquitin-protein ligase  TraesCS3B01G018700 Glutathione S-transferase T3 TraesCS3B01G023800 Cysteine/Histidine-rich C1 domain family TraesCS3B01G024400 Cysteine-rich receptor kinase TraesCS3B01G025000 Glutathione S-transferase TraesCS3B01G025100 Xyloglucan galactosyltransferase TraesCS3B01G026900 E3 ubiquitin-protein ligase TraesCS3B01G027600 SKP1-like protein TraesCS3B01G028000 E3 ubiquitin-protein ligase TraesCS3B01G028100 Cell wall invertase TraesCS3B01G028200 Beta-fructofuranosidase, insoluble protein TraesCS3B01G029000 Ubiquitin-conjugating enzyme E2 Ulp1 protease family TraesCS5B01G374300 Glutaredoxin family protein, putative TraesCS5B01G376000 Ubiquitin carboxyl-terminal hydrolase family TraesCS5B01G376800 Glutathione S-transferase TraesCS5B01G377700 Cysteine proteinase TraesCS5B01G377800 Cysteine protease TraesCS5B01G378800 ATP-dependent zinc metalloprotease TraesCS5B01G379000 Copper ion-binding protein TraesCS5B01G379900 Eukaryotic aspartyl protease family protein region [59,60]. The region ∼113-115 cM on chr 1A was associated with the grain phosphorus, nitrogen and copper concentrations as well as grain weight per ear [56,61]. The region on top of chr 3B was associated with grain-filling duration, grain yield, grain number per ear and plant height [43,56,57]. To date there has not been any QTL for grain Zn or Fe concentration reported for the 1A and 3B regions. The co-located MTAs reported in our study demonstrate that these loci controlling Zn, Fe and other micronutrients were either pleiotropic or closely linked suggesting that the identified markers could be used in MAS to potentially improve the concentrations of more than one important micronutrient at the same time. To date there has not been any yield QTL associated with the 1A or 5B thus marker-assisted selection in this area to improve multiple mineral concentrations can potentially be achieved without yield penalty.
At the 5B loci (95 cM), selection for increasing concentrations of Zn and other micronutrients would also lead to higher P concentrations, which may affect the bioavailability of those micronutrients. In mature grain, P is mainly stored as phytate (myo-inositol-1,2,3,4,5,6-hexakisphosphate, InsP6), which has the ability to complex essential micronutrients forming insoluble complexes that cannot be digested or absorbed by humans [62]. Thus, further studies on bioavailability would be required before those markers can be used for biofortification programs.

Prediction of candidate genes
The most recently published wheat genome sequence and gene annotation enabled further investigation into candidate genes located within those genomic regions, potentially responsible for the variation of the mature grain nutrient concentrations. The list of candidate genes located within the genomic regions controlling multiple elements is presented in Table 5. The markers on chr 5B (95.5 cM) that were associated with all five nutrients were all linked to a gene coding for a cysteine proteinase belonging to the papain-like cysteine protease family C1A (Gene ID: TraesCS5B01G377700). Cysteine proteinase genes were found being expressed during seed development and responsible for controlling the correct processing of storage proteins deposited in the endosperm of cereals [63,64]. The processing of storage proteins such as would affect its packing and density, which in turns would affect the amount of zinc and possibly other protein-bound metals. Several genes coding for other proteases located in this genomic region included ATP-dependent zinc metalloprotease, aspartyl protease, Ulp1 protease family. Other potential candidate genes located closely to the markers included Potassium transporter, E3 ubiquitin protein ligase, Copper ion-binding protein, Glutathione S-transferase, Glutaredoxin. The biological functions of the potential candidate genes in relation to improving grain micronutrient concentrations, however, requires further experimental validations to confirm.
The identified markers on chr 1A (∼113 cM) (associated with grain Zn, Fe and Mn concentration in 2014-15) were linked to the genes coding for Lipid phosphate phosphatase, Phospholipase D, Zinc finger protein-like, Ubiquitin carboxyl-terminal hydrolase 12 and Purple acid phosphatase. The genes coding for Sucrose synthase 6 (homologs of the Arabidopsis AT1G73370 gene) were also located closely to the markers.
The region on top of chr 3B (11-14 cM) harboured genes controlling flowering time (Flowering locus T), which explains why it was associated with grain filling duration in previous study [43]. Since this region may control some plant development traits and grain yield [57], selection for higher micronutrient densities at those loci may suffer yield penalty and should be avoided.

Relationship between grain and rachis nutrient concentrations
In addition to mature grain, we were set out to investigate the genetics controlling mineral accumulation in developing grain and rachis at an early stage of grain development (10d after anthesis). It is important to understand the mechanism of Zn and other micronutrient loading into the grain to increase its efficiency. Grain mineral supplies come from the direct uptake from the soil and/or the remobilization of stored minerals in vegetative tissues during grain filling [65]. The relative importance of these processes differs depending on the availability of the micronutrients in the soil, fertilization schemes and plant varieties [24,25,66]. In our study, positive correlations were found between mature grain, immature grain and immature rachis for three micronutrient concentrations; Zn, Cu and Mn in both seasons (P < 0.05). For Fe and P, the correlations between sample types differed between growth seasons. The majority of the correlation coefficients between sample types were rather moderate (r ≤ 0.35) with the only exception of Cu concentrations between mature grain and immature gain ( Table 2). These correlation results indicate that the rachis was likely an important source for mineral translocation into the grain during the reproductive stage, particularly for Zn, Cu and Mn. However, there may be additional important factors controlling the final amount of these minerals in the grain. A possible controlling factor could be the transport barriers for minerals between the rachis and the grain during grain filling. Previous studies showed that there were two possible barriers two possible barriers exist to control Zn transport into wheat grain: (i) between rachis and grain and (ii) between endosperm and crease vascular tissue [67]. The control of Zn and other metal transport is a complex process involved a large number of genes and gene families. The gene families involved in the uptake, transport, and accumulation of Zn, Fe, Cu and Mn in plants included (but not limited to) ZIP (Zinc-regulated transporter (ZRT), Iron-regulated transporter (IRT)like protein) family, YSL (yellow stripe-like) proteins, HMA (heavy metal transporting ATPase) family, MTPs (metal tolerance proteins), COPT/Ctr (Copper transporter) family, metallotheioneins [68][69][70]. Another possible contributing factor could be the capacity of the seed acting as a sink for mineral storage. In cereal grain, trace elements (Zn, Fe, Mn and Cu) are bound to various molecules/peptides/proteins with high metal binding affinity such as phytate (or phytic acid), ferritin, nicotianamine, phytosiderophores, metal-containing enzymes (e.g. Cu/ Zn superoxide dismutase), metalloproteins and storage proteins [71][72][73][74][75][76]. The de novo synthesis of these sink molecules in grain tissue during grain filling would directly affect the amount of minerals that can be stored in the grain.
This study identified multiple genomic regions harbouring co-located MTAs between mature grain and immature grain and rachis nutrient concentrations. Interestingly, the three genomic regions harbouring multi-trait MTAs in mature grain on chr 1A (113 cM), 3B (11-14 cM) and 5B (95 cM) were all detected in immature grain. Those on 3B and 5B were also found in immature rachis. Another important region was on chr 2A (104−109 cM), which was associated with Zn concentration in all sample types and with Fe in mature grain and immature rachis. The markers on chr 2A were linked to a gene coding for a Zn transporter belonging to the Cation efflux family (homologous to the Arabidopsis At2g04620 gene). In Arabidopsis, it was found responsible for sequestration of excess Zn in the cytoplasm into vacuoles to maintain Zn homeostasis [77]. Several genes coding for Yellow stripe-like transporter 1, 9 and 12 families, Vacuolar iron transporter (VIP), Copper transporter were also located in this region. It's likely that those genes could play important roles in transporting metals such as Zn, Fe, Cu and Mn up to the rachis and into grain. Further validation studies on those genes would be required to confirm the significance of the candidate genes in improving grain nutrient concentrations.

Stable MTAs over seasons
The MTAs on 5B at 95 cM were found associated with four nutrient concentrations (Zn, Mn, Cu, P) in the grain in both growing seasons. At this locus, the allele associated with increased mineral concentrations was present in 35 accessions (∼11 % of the panel) indicating this was a rare allele, probably originating from an uncommon genetic pool. The results demonstrate the success of using exotic genetic materials to improve mineral concentrations in elite wheat varieties. In addition, the markers on chr 6D at ∼87 cM was consistently associated with grain Mn over seasons. The consistency of those marker effects indicates that they would be useful for marker-assisted breeding for the improvement of micronutrients in wheat grain.
The number of MTAs that were consistently detected over two seasons accounted for a small percentage of the total MTAs. It is likely attributed to environmental effects. The 2014-15 growing season tended to have lower average temperatures around anthesis and higher rainfall through seed maturation phase when compared to the 2015-16 season (Fig. S5). These factors could have considerable effects on plant growth, development and grain nutrition. On average the panel flowered 10 d earlier and had significantly higher grain Fe, Cu and P concentrations in 2014-15 compared to 2015-16 (Table 1). Having higher rainfall in 2014-15 could increase water availability between irrigations, which in turn influencing on nutrient availability from the soil for plant uptake, longdistance transport and remobilization within the plant. Temperature and daylength are two important factors influencing flowering time in wheat [78], which in turn could affect plant nutrient-use efficiency, grain yield and grain nutrient levels [79,80]. Earlier senescence was found to be associated with enhanced grain protein and micronutrient contents including Fe, Zn and Mn [81,82]. Our results show significant environmental effects on the genetics controlling grain nutrient levels, which have also been reported in previous studies in wheat and rice [42,83,84]. Although consensus MTAs can generally be considered as more favourable for marker-assisted selection, some MTAs detected in one environment may lead to important discoveries [39,42]. These MTAs may be used to mine candidate genes through network analysis using bioinformatics analysis and/or experimental validation.

Conclusion
Results presented here are the first to report on the genetic basis controlling multiple micronutrient levels in the grain and vegetative tissues at two different stages of grain development using multi-locus GWAS models. The positive outcomes were attained through the exploitation of an association mapping panel containing rich genetic diversity. Our results signified the importance of the genomic region on chr 5B in controlling multiple grain micronutrients. The markers consistently associated with Zn, Fe, Mn and Cu were identified together with the beneficial alleles and their estimated effects. These results could potentially benefit biofortification programs to breed for high micronutrient concentrations in wheat varieties without affecting yield. Multiple candidate genes were identified, which play various roles in controlling mineral accumulations in wheat grain, especially Zn transporters, cysteine proteases and pectin lyases. This shows that increasing both metal transport efficiency and the strength of the protein sink are essential to improve micronutrition in wheat grain. Further gene functionality studies would be helpful to validate the significance of the candidate genes in breeding for higher micronutrients.

Author contribution statement
Suong T. Cu: carried out data collection and analyses, performed genetic mapping, wrote and edited the manuscript.
Georgia Guild: collected samples from the field, prepared samples for nutrient analysis and reviewed the manuscript.
Alison Nicolson: dissected and prepared samples for nutrient analysis samples.
Velu Govindan: oversaw field trial design and management, carried out the genotyping and genetic map construction and reviewed the manuscript.
Ravi Singh: oversaw field trial design and management. James Stangoulis: supervised sample collection, preparation and analyses, reviewed the manuscript.

Funding
This work was supported by the Australian Research Council and HarvestPlus through an ARC-Linkage Project (grant number #LP150101242).

Declaration of Competing Interest
The authors declare that there are no conflicts of interest.