A Major Facilitator Superfamily Transporter is a putative candidate gene for nutrient mineral (Ca, K, Mg, Mn, P and S) accumulation in bread wheat grains

Highlight: A multi-locus GWAS for a diverse wheat panel, revealed several putative candidate genes for mineral accumulation in wheat grains with a transmembrane transporter as the candidate gene. Abstract Here we report a multi-locus genome-wide association study for a set of 369 diverse wheat ( Triticum aestivum L.) genotypes, which were genotyped by 90k iSELECT Infinium and 35k Affymetrix arrays and yielded 15,523 SNPs. The panel was grown under field conditions for three consecutive years: 2015, 2016 and 2017. ICP-OES (Inductively coupled plasma atomic emission spectroscopy) was used to measure the concentration of six nutrient minerals in wheat grains including: Ca, K, Mg, Mn, P and S. Wide ranges of natural variation among the genotypes in nutrient minerals concentrations were detected. The phenotypic correlation showed a strong positive correlation among the nutrient minerals except for K that showed opposite correlation trends with other nutrient minerals. The genetic association analysis detected eighty-six significant marker-trait associations (MTAs) underlying the natural variation in nutrient minerals concentration in grains. The major MTA was detected on the long arm of chromosome 5A at 698,510,027 bp that showed a pleiotropic effect on Ca, K, Mg, Mn, and S. Further significant MTAs were distributed among the whole genome except chromosomes 3D and 6D. We identified putative candidate genes which are potentially involved in metal uptake, transport and assimilation. TraesCS5A02G542600 gene at chromosome 5A (698,507,247- 698,511,217 bp) annotated as transmembrane transporter activity and belonging to major facilitator superfamily transporter is a putative candidate gene for Ca, K, Mg, Mn, and S grain concentrations. The allelic variation at this gene showed that the T allele increased the concentration of nutrient minerals in grain except for K . The expression analysis revealed that the TraesCS5A02G542600 gene is highly expressed in seed coat followed by peduncle, awns, and lemma . Furthermore, the genomic prediction findings indicate that genomic selection may be useful for the genetic improvement of nutrient minerals accumulation in wheat. Our study provides crucial insights into the genetic basis of nutrient minerals variation in wheat and serves as an important foundation for boosting nutritional value for further genetic and molecular mechanisms studies controlling nutrient minerals accumulation in wheat grains.


Introduction
Cereal-based foods represent the largest proportion of the daily diet consumed worldwide, while wheat (Triticum aestivum L.) is the primary protein source in the developing countries with 2.5 billion consumers in 89 countries (https://wheat.org/). However, wheat is inherently containing a low proportion of nutrients. The term "nutrients" refers to a broad list of minerals and vitamins that play an important role in the human body including the biological functions and cell physiology where the body needs these nutrients in the appropriate amount to function properly (Pujar et al., 2020;Tapiero et al., 2003).
Most of the wheat research programs focus on increasing the yield production while ignoring grain quality that leads to depletion in the nutrient composition and nutritional value. Thus, a high incidence of nutrient deficiencies and malnutrition is occurring particularly in the countries relying on cereal diet where bread wheat forms 60% of the daily calorie's intake. For example, the World Health Organization (WHO) records that in 2012, 162 and 99 million children were stunted and underweight, mainly because of insufficient intake of essential nutrients (WHO, 2013). Therefore, to develop high concentration of nutrients in wheat grains also is an imperative need to overcome malnutrition in the next generations.
Wheat is the second most important staple crop worldwide; it covers nearly two hundred million hectares, providing one-fifth of the total calories and feeding more than one-third of the world's population (Godfray et al., 2010). It is expected that the food demand particularly for wheat will be 35%-40% higher by 2050 (Ray et al., 2013). Therefore, to meet the increasing food demand for the growing population and achieving food security, recent research starts to focus on identifying the varieties having high quality including the nutritional aspects using the natural variations, (Alomari et al., 2017;Yu and Tian, 2018) that offer promising strategies to develop varieties meeting the nutritionally limited populations and to sustain nutrition security.
Recent developments in molecular markers and genome sequencing technologies, as well as the recent release of the high-quality wheat reference genome sequence, allow plant researchers to characterize the genetic basis of complex phenotypic traits by using hundreds of thousands of genetic markers in association mapping and quantitative trait loci detection (Appels et al., 2018).
Applying such high-density single nucleotide polymorphic (SNP) marker arrays with a suitable approach like genome-wide association study (GWAS) can allocate the robust QTLs in order to detect genes underlying complex phenotypic traits (Alqudah et al., 2020;Rasheed and Xia, 2019). Nutrients accumulation in wheat grains are complex inherited traits that are controlled by several factors such as nutrients absorption from the soil, uptake by the roots, translocation, assimilation, and remobilization to the grains (Sperotto et al., 2014). Thus, the identification of the loci underlying high calcium (Ca), magnesium (Mg), manganese (Mn), phosphorus (P), potassium (K) and Sulphur (S) in wheat grains will assist in genetic biofortification through marker-assisted selection. Recently, another approach based on using molecular markers covering the whole genome, the so-called genomic prediction (GP) approach has become popular in plant breeding. This approach is based on using the genome-wide marker information to predict the breeding value of complex traits in order to accelerate the breeding programs (Meuwissen et al., 2001). GP ability and accuracy may vary among genomic prediction methods that use different assumptions and treatments of marker effects and models (Desta and Ortiz, 2014).
In the current study, we aimed to investigate the natural variation in nutrient mineral accumulation in grains of elite European hexaploid winter wheat. Using high-density SNP arrays, we also targeted to detect stable genomic regions associated with the natural variation of studied minerals based on multi-locus GWAS and to predicate the genomic prediction ability and breeding value. The ultimate goal is to uncover the most relevant candidate genes potentially involved in controlling mineral accumulation in wheat grains. The information which comes from this scientific work will help wheat breeders for genetic biofortification and the development of new wheat varieties with improved nutritional values and assist in ensuring nutritional security.

Plant material
Wheat germplasm comprised 369 elite European wheat registered varieties including 355 genotypes of winter wheat and 14 genotypes of spring wheat (Supplementary Table S1). The genotypes were mostly from Germany and France in addition to other European countries. Field experiments were conducted at IPK, Gatersleben, Germany during the 2014-2017 year (2014/2015, 2015/2016, and 2016/2017) where the whole set of the germplasm was sawn at each year. A plot with a size of 2 × 2 m was used for each genotype with six rows spaced 0.20 m apart and more details were described in a previous study by (Alomari et al., 2017). Standard agronomic wheat management practices were subjected to the soil.

Grain digestion and nutrients measurement
Mineral concentration measurements were done for the whole set of wheat genotypes for each year. For each genotype sample, 50 kernels were randomly selected and counted using a digital seed analyzer/counter Marvin (GTA Sensorik GmbH, Neubrandenburg, Germany) and a thousand-grain weight (TGW) was estimated as well. The samples were grinded using a Retsch mill (MM300, Germany) then the ground wheat samples were dried at 40°C in the incubator and left overnight. Minerals such as Ca, K, Mg, Mn, and S were measured by inductively coupled plasma optical emission spectrometry (ICP-OES, iCAP 6000, Thermo Fisher Scientific, Germany) combined with a CETAC ASXPRESSTM PLUS rapid sample introduction system and a CETAC autosampler (CETAC Technologies, Omaha, United States). Each wheat genotype ground sample was prepared for wet digestion in 2 ml nitric acid (HNO3, 69%, Bernd Kraft GmbH, Germany) using a high-performance microwave reactor (UltraClave IV, MLS, Germany). Then the samples were filled up to 15 ml final volume with de-ionized distilled water (Milli-Q® Reference System, Merck, Germany). Element standards were prepared from Bernd Kraft multi-element standard solution A c c e p t e d M a n u s c r i p t 5 (Germany). An external standard for the measured minerals and Y (ICP Standard Certipur®, Merck, Germany) were used as internal standards for matrix correction.

Statistical analysis
Descriptive statistical analysis for each mineral among the years was calculated and the significant differences among genotypes and years were determined by calculating the analysis of variance (ANOVA) at a probability level of P ≤ 0.05. Pearson correlation coefficient was used to evaluate the relationships among the measured parameters of data at p-value 0.05. ANOVA and Pearson's correlation coefficient were calculated using GenStat v19 software (VSN International, Hemel Hempstead, Hertfordshire, United Kingdom). Through using the software GenStat v19 software (VSN International, Hemel Hempstead, Hertfordshire, United Kingdom), best linear unbiased estimates (BLUEs) over three years were calculated by restricted maximum likelihood (REML) analysis with a mixed linear model (MLM) where the genotype considered as a fixed effect and the environment as a random effect.. BLUEs was calculated for each genotype of each trait across the years (2015, 2016, and 2017).
Broad-sense Heritability (H 2 ) was calculated for each trait using the formula: Where VG is the variance of the genotype, Ge represents the variance of the residual and nE is the number of years.

Genotyping and markers quality control
Our wheat population was genotyped using two marker arrays: a 90K iSELECT Infinium array (Wang et al., 2014) including 7761 markers and a 35K Affymetrix SNP array including 7762 markers (Wang et al., 2014) (Axiom ® Wheat Breeder's Genotyping Array, https://www.cerealsdb.uk.net/cerealgenomics/) and these two arrays were genotyped by SGS-TraitGenetics GmbH company, Gatersleben, Germany (www.traitgenetics.com) and described for this germplasm by Alomari et al. (2019). ITMI-DH population was used as a reference map (Poland et al., 2012;Sorrells et al., 2011) to anchor the SNP markers of the 90K and 35K arrays. To obtain high-quality makers, the SNPs in chips underwent a quality control and filtration process, by removing the markers with ≤3% missing values, a minor allele frequency (MAF) of ≤3%, and markers with unknown chromosomal positions. Then, we used the physical position of wheat genome sequence RefSeq v1.0 for the SNPs.

GWAS and GP
GWAS analysis was calculated by using Genomic Association and Prediction Integrated Tool (GAPIT) in R software (Lipka et al., 2012). First, GWAS analysis was computed by the mixed linear model (MLM) which took into account variance-covariance kinship matrix and PCA and accomplished by incorporating the phenotypic and genotypic dataset. Moreover, the kinship matrix was calculated using the VanRaden method (VanRaden, 2008) to determine A c c e p t e d M a n u s c r i p t 6 relative kinship among the sampled individuals. Both PCA and kinship matrix were used for population correction and stratification.
Another recent powerful GAPIT model known as Fixed and random model Circulating Probability Unification (FarmCPU) was applied to our data analysis. FarmCPU model was applied by considering the incorporates of multiple markers simultaneously as a covariate in a fixed-effectect model and optimization on the associated covariate markers in a random effect model separately that empowered us to avoid any false-negative and control the falsepositive associations by preventing model overfitting (Liu et al., 2016). Thus, FarmCPU is a powerful tool with less false-positives than MLM. The selection of an appropriate model and threshold are important steps in identifying markers that are truly associated with specific traits and which could be located within or very close to genes that control the trait variation while controlling both false-positive and false-negative associations. To determine which of the tested models best fit the data, we plotted the quantile-quantile (Q-Q) plot which was drawn based on the observed and expected -log 10 (P) values. Then, based on the GWAS output of three models (GLM, MLM, and FarmCPU), the number of significant associations, and the resulted QQ-plot, we selected the FarmCPU. A threshold P-value 0.001 equal tolog 10 (P) ≥ 3 was used to determine the significance of marker-trait associations (MTAs), and then Bonferroni correction at P < 0.05 was used to adjust the -log 10 (P) threshold equal to 5.49 value for the studied traits.
Marker effects (positive/negative) and phenotypic variance explained by the associated markers (R 2 ) were taken out from GWAS results.
Ridge regression best linear unbiased prediction (RR-BLUP) and Bayes-Cπ methods were used to evaluate the genomic prediction (Habier et al., 2011;Meuwissen et al., 2001). Both methods are implemented in R software. Fivefold cross-validation was applied to the complete set where the set was randomly divided into five subsets, in which four of the five were used as estimation sets, and the remaining ones were used as the test set. The whole process was repeated 100 times and then the mean value was used as the final prediction ability. The prediction ability was calculated from the correlation between observed and predicted values.

Genes' identification, annotation, and expression analysis
Significant markers and the markers located within the LD region (r 2 ≥0.2) were considered for BLAST. The sequences of the identified makers were obtained from the wheat 90k (Wang et al., 2014) and 35k database (Allen et al., 2017). Markers sequences were BLASTed against the recently released IWGSC RefSeq v1.1 genome by Ensemble Plants (http://plants.ensembl.org/Triticum_aestivum) to get their gene annotation. The expression profile of all the putative candidate genes associated with the identified SNPs was checked using the published RNA-seq expression database of wheat in WheatGmap web tool (https://www.wheatgmap.org) .
A c c e p t e d M a n u s c r i p t 7

Genetic variation and correlation
Wheat genotypes showed a wide variation through the three years and BLUEs for the studied traits which followed approximate normal distributions (Figure 1a and Figure S1). Statistical analysis for the studied minerals including the maximum, minimum, and mean values were calculated and presented in Supplementary Table (S2). The variation ranged from 208.5 to 797.2 μg.g -1 for grain Ca, 3495 to 6727 μg.g -1 for grain K, 963.9 to 1988 μg.g -1 for grain Mg, 23.37 to 62.2 μg.g -1 for grain Mn, 2943 to 5807 μg.g -1 for grain P, and 974.8 to 2368 μg.g -1 for grain S concentration (Supplementary Table S2).The nutrients grain results showed that Ca and S concentrations were significantly higher in 2015, while K, P and Mn concentrations were higher in 2016 whereas, Mn concentration did not change among the three years and this variation among years could be explained by environmental influences (Supplementary  Table S2). Correlation analysis was performed for nutrient traits in addition to thousand kernel weight (TKW) based on BLUE values ( Figure 1B). The correlation coefficients (r) and their significance levels showed a strong positive correlation between Mg with P and between Mg with Mn (r= 0.69 and 0.63 respectively, P < 0.05) while, a moderate positive correlation was detected between Ca with Mg, Mn, P and S (0.25, 0.43, 0.25 and 0.27 respectively, P < 0.01) as shown in Figure 1b. A negative correlation was detected between K with Ca, Mg, Mn, P and S while this correlation was not significant for both Ca and P ( Figure 1b). All the nutrient minerals showed a neutral correlation with TKW. Our results showed that the five genotypes with highest mineral concentrations are: Isengrain, Inoui, Nirvana, Exotic and Lona, indicating that these genotypes are produced high nutritive grains over years.
The ANOVA revealed significant effects among genotypes (P < 0.001) and the years (P < 0.001) for all traits except Mn (Supplementary Table S3). The broad-sense heritability (H 2 ) values were high and ranged from 0.72 for Mn to 0.87 for Ca (Supplementary Table S3). The phenotypic results suggesting that there is ample natural variation in mineral concentration which are predominantly genetically controlled with low influence by environments that might lead to detect shared genomic regions among the studied traits.

Association analysis and Genomic Prediction
GWAS analysis was performed by FarmCPU using 15, 523 SNP markers to reveal the genetic basis of the accumulation of Ca, Mg, Mn, P, S, and K in wheat grains. GWAS output revealed a quantitative genetic nature of our traits of interest. Most of the markers were mapped on the B-genome followed by A-and D-genomes. Briefly, the MTAs were distributed among the whole genome except 3D and 6D where 86 significant [-log 10 (Pvalue) ≥3] associations were detected, where the highest number of significant markers were found on chromosome 2B (12 markers) followed by chromosome 5B (9 markers) and 3B (8 markers). The MTAs were identified using the estimated BLUE values across the three years (Supplementary Table S4A). Seventeen MTAs were above the Bonferroni threshold which is equal -log 10 (P-value) 5.42 (Figure 2A and Supplementary Table S4B Table S4A). The phenotypic variance explained by each SNP ranged from 0.04% to 10.54% (Supplementary Table S3A). Markedly, RAC875_c8642_231 that showed the highest -log 10 (P-value) on the studied traits ( Figure 2A and Figure 3A) and had a positive effect on Ca, Mn, Mg, and S (~ +26, 1.8, 49, 28 μg.g -1 for Ca) and negative effect on (~ -117 μg.g -1 for grain K). Meanwhile, this marker could explain ~10% of the variation in Ca and ~ 3.4% in K. This indicated that there is allelic variation at this marker among the genotypes that lead to variation in mineral concentration (Supplementary Table S5).
The Q-Q plots illustrating observed associations between SNPs and grain nutrient concentrations compared to expected associations are presented in Figure 2B.
Statistical methods for predicting the breeding value of the studied traits such as ridge regression best linear unbiased prediction (RR-BLUP) and various Bayesian models are used for developing prediction models. In our study, the mean of prediction ability values that resulted from using two different statistical methods (RR-BLUP and Bayes-Cπ) were almost in the same range. The values using Bayes-Cπ were ranged between 0.31 -0.59 and for RR-BLUP were ranged between 0.30 -0.57. The highest value was for Mg followed by Mn ( Figure 4).

Candidate gene detection and gene expression
Significant SNPs identified in the current study were used to predict gene models located on the wheat genome using Chinese Spring RefSeqv1.1. Candidate gene transcripts and their corresponding annotation information were obtained from the website of Wheat Ensemble (http://plants.ensembl.org/Triticum_aestivum). The candidate genes were annotated as transmembrane transporter activity, protein kinase, ATPase-coupled cation transmembrane transporter, metal ion binding and magnesium ion binding.
All of the potential candidate genes and their corresponding annotations within all detected loci were listed in Table 1. Most notably from the candidate gene list, TraesCS5A02G542600 with annotation as transmembrane transporter activity and underlying the natural variation of Ca, Mn, Mg, S, K (Table 1) is located on chromosome 5A (698,507,247-698,511,217 bp, Figure 3A). Interestingly, the SNP-marker RAC875_c8642_231 that had the highest association value among all SNPs is located on chromosome 5A (698,510,016 bp) inside the gene TraesCS5A02G542600 ( Figure 3A). Based on the gene structure, the RAC875_c8642_231 marker is physically located inside the exon 3 of TraesCS5A02G542600 ( Figure 3C). The polymorphism analysis of RAC875_c8642_231 alleles showed that the genotypes of the population are either carrying the C allele (337 genotypes) or carrying the T allele (32 genotypes) as shown in Figure 3C and Table S1. The genotypes with the T allele have significantly (P < 0.05) higher nutrient concentration values in wheat grain except for K ( Figure 3D). This finding showed a positive effect of the T allele on TKW as well, confirming the effect of the gene-based alleles on the nutrient concentration and TKW in wheat grains ( Figure 3D). The genotypes which are carrying the T allele are originally coming from France ( Figure 3E).
A c c e p t e d M a n u s c r i p t 9 The expression analysis of candidate genes in different grain tissues and developmental stages showed a wide range of expression for the genes ( Figure 5). Generally, gene TraesCS5B02G012300 at 5B, and gene TraesCS3B02G006700 at 3B, showed the highest expression in most of the organs and tissues, while the expression of gene TraesCS5A02G542600 was very high in seed coat followed by peduncle, awns, lemma, ovary, stem, spikelets, grain, leaf ligule, shoots, rachis, endosperm, glumes, flag leaf blade, internode, spike, aleurone layer, and roots ( Figure 5). Their expression showed that they play a vital role during growth, filling and development. TraesCS3B02G006700 gene is linked with Ca and it is highly expressed in lemma and seed coat in addition to flag leaf blade, leaf ligule and grain. One of the significantly associated genes was very low expressed (TraesCS6B02G002500) in the organs compared with the aforementioned genes and removed from the expression analysis output. RNA-seq are presented as log 2 TPM where TPM is transcripts per million.

Phenotypic
Breeding wheat varieties with essential nutrients through genomics breeding tools is a promising approach for the production of high nutritional quality wheat grains. Therefore, detection of natural wide phenotypic variation in the source material is the primary step toward genetic biofortification (Yu and Tian, 2018). Overall, a wide natural variation in wheat grain mineral concentration was observed in the current study which also is in line with other studies on wheat grain Ca, Mg, Mn, P, K, and S (Bhatta et al., 2018). High H 2 values were found in the studied traits with H 2 more than 70% for Ca, Mg, Mn, Mn, P and K demonstrated that the major part of the variation in the mineral is genetically controlled. The same trend of findings was found in recombinant inbred lines (RILs) of tetraploid wheat by Peleg et al. (2009) in Ca, Mg and S, while a moderate H 2 was found for Mn, P, K in a wild barley NAM population (Herzig et al., 2019). The positive correlation among minerals in the current study suggested the presence of common genetic factors affecting the accumulation of these micronutrients in grains. Similar trends of correlations were observed by Shi (Shi et al., 2013) in hexaploid wheat. The lack of correlation between TKW and other traits implied that there is no clear effect of TKW on the mineral concentrations in grains.

GWAS and GP
Most of the genetic studies in wheat are using QTL mapping to investigate the genetic basis of nutrients accumulation (Pu et al., 2014;Shi et al., 2013). In these previous studies, biparental crosses were used to identify QTLs and genes associated with the mineral concentration. However, the resulted QTL mapping is not in high resolution and therefore, there is a tremendous interest in using GWAS as an alternative to QTL mapping that can identify alleles represented in a broader set of germplasm with high resolution (Alomari et al., 2018;Alomari et al., 2017;Alomari et al., 2019;Alqudah et al., 2020). Furthermore, the implementation of GWAS analysis is useful to identify molecular markers that are tightly linked to genomic regions underlying natural variation of nutrients such as Ca, Mg, Mn, P, K A c c e p t e d M a n u s c r i p t 10 and S that can be used for enhancing the efficiency of biofortification using genomics assisted breeding (Collard et al., 2005).
Moreover, most of the association mapping studies undertaken for different complex traits are using single-locus GWAS models (GLM/ MLM) and these models require further correction e.g. multiple testing corrections to control false positives (Price et al., 2010;Pujar et al., 2020;Yu et al., 2006). However, the corrections are often not stringent enough or too stringent and overcompensate for population structure and kinship, which could lead to overcorrection where some potentially important MTAs can be missed (Alqudah et al., 2020). To overcome this issue, FarmCPU has been developed and shown to be much more powerful, efficient and superior in reducing false positive/negative associations, particularly in complex traits by incorporating multiple markers simultaneously as covariates in a stepwise MLM to partially remove the confounding between testing markers and kinship (Liu et al., 2016). Therefore, in our study, we applied the FarmCPU model to use its advantages in GWAS analysis.
In the current study, many genomic regions harbored markers that were associated with nutrient concentration were observed. In total, 19 MTAs were detected for S, 15 MTAs for both of K and P, 14, 13 MTAs for Mg and Ca respectively while 10 MTAs were found for Mn. The most significant association was detected at the end of the long arm of chromosome 5A and linked with RAC875_c8642_231 marker at the position of ∼114.5 cM (698,5100,16 bp), which was linked with all five nutrients except P. P is positively correlated with phytate content in plant and it is known as an anti-nutrient compound that negatively influenced other minerals absorption in the human body (Stangoulis et al., 2007). Therefore, this MTA is a good candidate in breeding wheat to improve grain many mineral concentrations simultaneously. The present study also observed that grain K concentration is negatively correlated with grain Ca, Mn, Mg and S ( Figure 3D, E). MTA for grain K concentration also negatively linked with the grain Ca, Mn, Mg and S MTA at the 5A locus (Table S3).
The QTL co-localization on chromosome 5A for some of the minerals such as nitrogen, Fe, Cu, Mg, and K (Shi et al., 2013) and Fe, Zn, Cu and Mg but not P have previously been reported in wheat, as well Cu et al., (2020) detected another co-located locus for Zn, Fe, Cu, P and Mn which is on chromosome 5B ~95 cM (Cu et al., 2020). Therefore, it seems that the QTL on chromosome 5A plays a vital role in mineral accumulation in wheat which needs further functional validation. The results showed that most of the highly nutritive genotypes over the years are originated from France and characterized by having awns of which the highest five genotypes are: Isengrain, Inoui, Nirvana, Exotic and Lona. Therefore, including highly nutritive genotypes in breeding programs may help in enhancing the minerals accumulation in wheat grain. We extended our analysis to include the GP method which has a practical role in improving the breeding efficiency of quantitative and complex traits. Our wheat panel predictability results showed low to moderate values (Figure 4), which agrees with another report that obtained low to moderate predictability values for the macro and micro-nutrients in wheat (Manickavelu et al., 2017). Relying on these findings, GP may be considered as a promising approach for enhancing nutrient minerals in wheat especially when A c c e p t e d M a n u s c r i p t 11 larger size germplasm panels with increasing the numbers of markers are used to have more accurate estimates of breeding values.

Candidate genes and expression
The International Wheat Genome Sequence Consortium (IWGSC) and gene annotation V1.1 enabled further investigation into candidate genes potentially responsible for the variation of the grain nutrient concentrations. The list of candidate genes located within the genomic regions controlling Ca, K, Mg, Mn, P and S are presented in Table 1. As mentioned above, we focused on transport-related genes which are potentially involved in nutrient mineral uptake, translocation, and/or homeostasis in grains. Only the high-confidence (HC) candidate genes that are linked to the transcription factor regulators, transporters, and grain development will be demonstrated.

Ca
The most significant MTA which was detected in the GWAS output found to be associated with five nutrients (Ca, K, Mg, Mn, and S) and located on chromosome 5A (114.5 cM) and linked to a gene coding for a transmembrane transporter belongs to the major facilitator superfamily (MSF) (Gene ID: TraesCS5A02G542600). This gene family is known as one of the two largest superfamilies of membrane transporters and they are acting as uniporters or symporter for different substances (Niño-González et al., 2019). Interestingly, three recent studies in wheat explored the role of TraesCS5A02G542600 in the inhibition of awn formation; while TraesCS5A02G542600 is located at the genetic locus of awn suppression, a closely linked gene was considered as the most likely candidate for awn formation (DeWitt et al., 2020;Huang et al., 2020;Würschum et al., 2020). Further studies are needed to shed light on its role underlying nutrient mineral uptake and accumulation in wheat grains. The expression of gene TraesCS5A02G542600 was very high in the grain layers such as seed coat, lemma, aleurone and within the grain in addition to other grain parts like endosperm, peduncle, spikelets and spikes (Table S5) which demonstrated our hypothesis of the involvement of this gene in mineral accumulation in wheat grains.
Another candidate gene for Ca accumulation (TraesCS7A02G169100, (125,712,948-125,715,936bp)) has an annotation as transmembrane transporter and was encoded as WALLS ARE THIN1 (WAT1-related protein). This gene was found to be involved in secondary cell wall formation (Kaur et al., 2017) and it has been demonstrated that Ca 2+ plays a major role in determining the structural rigidity of the cell wall where high Ca 2+ concentrations rigidify the wall and make it less plastic and where low Ca 2+ levels make the cell wall more pliable and easily ruptured (Hepler, 2005). The TraesCS7A02G169100 gene is highly expressed in the aleurone layer, lemma, and within the grain (Table S5). A Diacylglycerol O-acyltransferase gene (TraesCS3B02G006700, (360,145,0 -360,720,8bp)) which annotated as transferase activity, catalyses the final step of the Triacylglycerol (TAG) synthesis. It showed significantly higher accumulation in Arabidopsis seedlings during nitrogen deprivation (Yang et al., 2011). Therefore, it may have a role in other nutrient A c c e p t e d M a n u s c r i p t 12 minerals accumulation such as Ca in wheat grains. The expression analysis for the TraesCS3B02G006700 gene showed high values for lemma seed coat, flag leaf blade, leaf ligule and grain that potentially validate our findings.
TraesCS5B02G403400 (580,100,440 -580,104,246bp) encoded the gene of semialdehyde dehydrogenase which is one of three enzymes constituting the gamma-aminobutyric acid shunt, a metabolic pathway that was associated with abiotic stress responses in durum wheat (Carillo, 2018). The pathway is controlled by three enzymes, the first of which is Glu decarboxylase (GAD) which was found to be a calcium/calmodulin-binding protein (Baum et al., 1996;Busch and Fromm, 1999).

K
Two Homeobox superfamily genes (TraesCS7B02G478200 (733,527,123-733,530,917bp)) and TraesCS7D02G540700 (630,527,242-630,530,005bp)) were linked with K and this superfamily is one of the transcription factor families that are involved in plant development, growth, and in the response to diverse stresses (Wei et al., 2019). Both of these two genes showed almost the same expression trend where the highest value is for the peduncle followed by stem, rachis, aleurone layer, seed coat, leaf ligule and spikelets, which support our data.
TraesCS3B02G590500 (815,476,561 -815,480,015bp) and TraesCS2D02G190600 (134,638,904 -134,645,251bp) were annotated as protein kinase which is a large superfamily that plays vital roles in plant development and stress tolerance, whereas only a limited number of protein kinases have been functionally studied in wheat (Wei and Li, 2019). Therefore, they may have a significant role in K transportation within the wheat plant. Several studies reported about the E3 ubiquitin-protein ligase gene (TraesCS4A02G352200, ( 627,813,928-627,816,836bp)) and its association with various minerals uptake and accumulation such as Mn, Zn and P in wheat and pearl millet (Cu et al., 2020;Pujar et al., 2020).

Mg
For Mg, we detected a P-type ATPase with gene TraesCS4B02G293600 (579,391,172 -579,398,644bp) on chromosome 4B and interestingly this gene was considered as transmembrane protein which play a crucial role in the transport of a wide variety of cations across membranes and which are vital for ion homeostasis and heavy metal detoxification (Axelsen and Palmgren, 2001). Very high expression values were detected for the TraesCS4B02G293600 gene in these parts: Peduncle, seed coat, spikelets, leaf ligule, and stem. On chromosome 7A, Ubiquitin specific protease (TraesCS7A02G498400, (688,960,269 -688,968,295bp)) was detected and annotated as metal ion binding.  A  gene  encoded  as  Tetratricopeptide-like  helical  domain  superfamily  (TraesCS4B02G024300 (173,621,92 -173,688,28bp), TraesCS5B02G042900 (475,842,53 -475,880,87bp)) was observed on both chromosomes, 4B and 5B; this gene is involved in root development and auxin polar transport as well as gibberellin signal transduction in Arabidopsis (Jacobsen et al., 1996;Zhang et al., 2015) which may indicate its involvement in Mn accumulation in wheat plants. These two genes were highly expressed in endosperm, flag leaf blade, lemma, grain, and glumes.
Another two genes were identified on chromosome 5B encoded as Phosphopyruvate hydratase (TraesCS5B02G012300, (123,242,87 -123,288,41bp)) and annotated as magnesium ion binding. (Silva et al., 2018) reported about the association between the Phosphopyruvate hydratase gene and enhancing the nitrogen metabolism in maize seedlings.

P
On chromosome 2A, a gene encoded as a Haloacid Dehalogenase (HAD) (TraesCS2A02G130200, (779,384,66 -779,418,46bp)) was observed in our panel. HAD gene enhances Phosphate accumulation where LePS2 was the first low Pi inducible HAD superfamily gene characterized in tomato (Baldwin et al., 2008;Baldwin et al., 2001). Another HAD gene, PvHAD1 showed low P specific induction and encodes a functional serine/threonine phosphatase (Liu et al., 2012). Two responsive HAD genes, AtPPsPase1 and AtPECP1 were reported to encode functional pyrophosphatase and phosphoethanolamine phosphatase in Arabidopsis. In rice, only two HAD genes, OsACP1 (Hur et al., 2007) and OsHAD1 (Pandey et al., 2017) were shown to be upregulated under P deficiency. These studies suggest important functions of HAD superfamily members in P accumulation in different plants, while there are no available studies on wheat.
A c c e p t e d M a n u s c r i p t 14 As well, another gene has a significant role in protein phosphorylation/dephosphorylation was observed and called as 191,195,889bp)) (Bergey et al., 2014).

Conclusions
There have been very few GWAS studies identifying the significant loci associated with nutrient mineral accumulation in wheat grains. Therefore, the MTAs identified in this study and candidate genes will be useful for the genetic improvement of wheat nutritional quality through marker-assisted selection. Moreover, this study also provides details on the range of phenotypic variation encountered within the European wheat germplasm. The results of GWAS show promising prospects in applied breeding to improve genetic gains for breeding of high nutrients. These results implied the possibility of simultaneous improvement of these traits in wheat grain by marker-assisted selection. However, further studies are needed for a better understanding of the relationship between each mineral nutrient in wheat and to supply more information for an efficient breeding.
M a n u s c r i p t 21 Tables   Table 1: Identification of candidate genes for the significant MTAs associated with nutrient mineral concentration in a germplasm set of 369 European wheat genotypes. Figure 1A: Grain nutrient minerals concentration (µg·g −1 ) in a germplasm set of 369 wheat genotypes.