Genome Wide Association Studies for Grain Micronutrient Concentration in Spring Wheat (Triticum aestivum L.)

1. Introduction

Wheat is (Triticum aestivum L.) one of the world's leading food crops and the cradle for the human civilization, marks an annual production of 797.86 million tons in about 223.01 million hectares [1]. It plays a fundamental role in providing global food security with an average of 19% calories and 21% of the protein required in daily human diet [2]. However, reports have shown that worldwide one out of 3 people suffer from some form of malnutrition. Almost 50% of deaths below the age of 5 years results from malnutrition, which majorly occurs in poor or middle-income nations [3]. Combating malnutrition is a global health challenge, for which United Nations declared 2016 to 2025 as United Nations Decade of Action on Nutrition. Since it is a decade to address malnutrition in every form and achieve the Sustainable Development Goals, studies and sustainable approach are required to remove this micronutrient malnutrition.

Calcium (Ca) in human body forms the structure of bones and is responsible for the normal body movements including tissue rigidity, strength, and flexibility, with undergoing continuous remodeling through constant deposition and resorption of Ca into the bone [4]. However, this starts declining in women due to hormone fluctuations and menopause. Moreover, children suffering with rickets does not have normal mineralization of growth cartilage, leading to irreversible deformities in the skeleton. Ca deficiency results in reduction of bone strength and osteoporosis. Another problem arising due to chronic Ca deficiency is osteomalacia and bone softening, which occurs in children and adults. Since Ca binds to fatty acids in the body, it reduces lipid absorption and therefore involved in lowering the risk of cardiovascular disease [5].

Another important micronutrient in human body, Magnesium (Mg), is a cofactor for more than 300 enzyme systems regulating various biochemical reactions, which includes blood glucose regulation, nerve and muscle functions, blood pressure regulation and protein synthesis [6]. It is involved in active transportation of calcium and potassium ions across cell membrane, which is required for nerve impulse conduction, muscle contractions, and normal heart rhythm. Essential for energy generation, glycolysis, and oxidative phosphorylation, it also plays an important role in DNA and RNA synthesis and development of bones by influencing the activities of osteoblasts and osteoclasts [7]. Initial symptoms of Mg deficiency include nausea, weakness, and loss of appetite. With deficiency being grim resulting into muscle contractions and cramps, abnormal heart rhythms, numbness seizures [8]. Since Mg influences the concentrations of parathyroid hormone and the active form of vitamin D, it is highly impactful in bone homeostasis. According to WHO, recommended Ca allowances range from 200-400 mg/day for infants to 1300 mg/day in children of 9-18 years and 1000-1200 mg/day in older adults. Similarly, for Mg it ranges from 30-75 mg in infants and around 400 mg in children and adults [3].

Molecular breeding is an efficient method for facilitating crop improvement targeted for micronutrient-related traits. Marker-trait associations (MTAs) associated with micronutrients and various other traits in wheat have been reported using Genome Wide Association Studies (GWAS) [9,10,11,12,13,14]. In GWAS, differentiation between true associations and false-positive MTAs resulting from population structure and kinship has been a big hurdle [15]. To overcome such associations, various statistical models have been developed, namely, single-locus mixed linear model and multi-locus models. In single-locus model incorporation of two confounding factors are considered for analysis as covariates [15]. However, this may lead to false-negative MTAs due to overfits, which results in missing of chances to uncover the loci associated with targeted traits [16]. For surpassing such incidences of false-negative MTAs, variety of multi-locus models have been developed such as, fixed, and random model circulating probability unification (FarmCPU) and Bayesian-information and linkage disequilibrium iteratively nested keyway (BLINK). FarmCPU model includes multiple markers simultaneously as the covariates, which eliminates the confounding effect of markers and kinship and therefore removing the false-positive MTAs without any compromise for true associations. It utilizes fixed-effect and random effect model for complete removal of cofounding. However, FarmCPU is time consuming where number of markers and individuals is high. In comparison, BLINK is statistically better and time-efficient using Bayesian information criterion. Unlike FarmCPU, BLINK utilizes linkage disequilibrium (LD) to eliminate the assumption that causal genes are evenly distributed across the genome [17,18]. Proper execution in GWAS reveal novel genomic regions associated with targeted traits and hence facilitate the identification of genes that control such traits. In the present investigation, 105 diverse wheat genotypes from the Wheat Association Mapping Initiative (WAMI) panel has been analyzed for grain Ca and Mg content using FarmCPU and BLINK models. Identification of genomic regions associated with controlling the two micronutrients was done, from which markers could be designed for their introgression in future breeding program.

2. Results

2.1. Variability for Grain Ca and Mg Concentration

Estimation of micronutrient concentrations of Ca and Mg in all 105 grain samples from the cropping season 2021-22 and 2022-23 was performed. Variations were observed in the concentration of the estimated micronutrients. Ca and Mg in two consecutive cropping season is presented in Figure 1.

Table 1. Combined ANOVA table for grain micronutrient in WAMI panel constituting of 105 wheat accessions.

Source	Df	Ca_SS	Mg_SS
Year	1	42.2*	476*
Entry	104	8.79*	9.68 *
Replication	2	10.3	60.5
Year× Entry	104	5.35*	7.61*
SS: Sum of square; Df: degrees of freedom; significant at ** p <0.001.

Table 1. Combined ANOVA table for grain micronutrient in WAMI panel constituting of 105 wheat accessions.

Source	Df	Ca_SS	Mg_SS
Year	1	42.2*	476*
Entry	104	8.79*	9.68 *
Replication	2	10.3	60.5
Year× Entry	104	5.35*	7.61*
SS: Sum of square; Df: degrees of freedom; significant at ** p <0.001.

Based on the samples obtained during 2021-2022, the mean GCaC for all three replicates was 398.62 ppm and, for 2022-2023, it was 380.42 ppm, with the final mean of 389.52 ppm (Table S1). The mean GMgC for the year 2021-2022 for all three replicates was 1211.87 ppm and, for 2022-2023, it was 1118.01ppm, with the final mean of 1164.94 ppm. Frequency distribution of variation for

Table 2. The statistical analysis for the GCaC and GMgC in the WAMI panel constituting of two consecutive years 2021-22 and 2022-23.

Statistics	2021-22		2022-23
	BLUP Ca	BLUP Mg	BLUP Ca	BLUP Mg
Heritability	0.66	0.73	0.68	0.78
Genotypic variance	2312.00	10528.4	2786.74	4323.17
Grand Mean	398.62	1211.88	380.42	1118.02
LSD	49.83	73.2	49.36	85.21
CV	8.36	3.88	8.55	5.34
LSD = Least significant difference; CV = Coefficient of variation

Table 2. The statistical analysis for the GCaC and GMgC in the WAMI panel constituting of two consecutive years 2021-22 and 2022-23.

Statistics	2021-22		2022-23
	BLUP Ca	BLUP Mg	BLUP Ca	BLUP Mg
Heritability	0.66	0.73	0.68	0.78
Genotypic variance	2312.00	10528.4	2786.74	4323.17
Grand Mean	398.62	1211.88	380.42	1118.02
LSD	49.83	73.2	49.36	85.21
CV	8.36	3.88	8.55	5.34
LSD = Least significant difference; CV = Coefficient of variation

2.3. Population Structure Analysis

The allele frequencies in the WAMI panel were found to be distributed evenly by principal components analysis (Figure 2A). PC1 explained 17.9% of the variance, whereas PC2 explained 5.5% and PC3 explained 3.9%. These findings imply that PC1, PC2 and PC3 collectively account for a significant amount of the panel's underlying genetic variation, making it possible to recognise and distinguish the subgroups that make up the WAMI population.

2.4. MTAs for the Target Trait Utilizing GWAS

To find out the significant associations of Ca and Mg content and SNP markers in the lines, GWAS analyses were carried out. A scatter plot of r² values of paired markers was developed, demonstrating genome-wide linkage disequilibrium (LD) degradation for 105 WAMI genotypes assessed using the Hill and Weir formula [20]. LD analysis was performed using Trait Analysis by aSSociation, Evolution and Linkage (TASSEL) data with a sliding window of 100 SNPs. The average genome-wide r² was 0.11, with LD decline starting at r² = 0.46 and ending at r² = 0.23 (Figure 2B). The LD decay curve intersected the half-decay and standard critical (r² = 0.3) lines at 23754184 and 11877092 basepair (bp), respectively. This establishes 11877092 bp as the genome-wide threshold distance for detecting linkage, as done in previous works [9,19].

The MTAs were investigated utilizing two multi-locus models, FarmCPU and BLINK, which identified 86 SNP loci with prominent associations with the trait studied and utilization of Manhattan plots and Quantile-Quantile plot (Q-Q) plots (Figure 3). Genotypic data was available for 26,814 SNPs, but only 21,132 SNPs were utilized for GWAS. We have found 49 SNPs significantly associated with Ca and located on the 13 chromosomes (1A, 1B, 1D, 2A, 2B, 2D, 3A, 3B, 3D, 5A, 6B and 7A), which is confirmed in the data provided for the consecutive years of cropping season using both the models (Table S2, S3, S4 and S5). About 37 SNPs were significantly found for GMgC located significantly on chromosomes 1B, 2A, 2B, 3A, 4A, 4B, 4D, 5B, 6A, 6B, 7A, 7B and 7D (Table S6, S7, S8 and S9).

2.4.1. MTAs for GCaC

For the initial crop year 2021–22, the MTAs were found to be distributed on multiple chromosomes, which were found notably on the chromosome 2A (4 MTAs), 3D (13 MTAs), 5A (2 MTAs), 6A (3 MTAs), 7B (6 MTAs), and 7A (5 MTAs). Interestingly, a dense concentration of 13 MTAs were observed between the physical position of 604.6-746.6 Mb on chromosome 3D, exhibiting a significant and rich association. This striking presence underscores potentially a genomic hotspot associated with Ca content. Presence of 4 MTAs on 7B was observed within the range of 35.4-37 Mb. Another profound presence of MTAs was on 7A, with three MTAs co-located on 52.9 Mb and one on 49.8 Mb. For the year 2022-23, these 7 MTAs were identified on 2B and 5A. Notably, these four MTAs were consistent in both the cropping season of 2021-22 and 2022-23, with each located on chromosome 2A, 6B and 2 MTAs on chromosome 5A (Table 3). There were four loci for GCaC, which showed presence in both the cropping years, Kukri_c11327_977 on chromosome 2A, wsnp_Ex_c2718_5038582 and RAC875_c9984_1003 both on 5A, and wsnp_Ex_c34597_42879718 on 6B.

2.4.2. MTAs for GMgC

In study for the year 2021-22, MTAs for GMgC extended on various chromosomes with notable presence on 4A (6 MTAs), 6A (3 MTAs) and 7A (4 MTAs). At a specific position of 99.2 Mb on chromosome 4A, two MTAs were found to be co-located. Moreover, two pair of co-located MTAs on 7A with positions at 159.5 Mb and 696.9 Mb were identified. In the cropping year 2022-23, prominent number of MTAs were observed on 2B (6MTAs), 4B (2 MTAs), 5B (3 MTAs) and 6A (4MTAs). The most prominent observation for MTAs associated with Mg, which were consistent in both the cropping season were on 6A and 6B. There were four loci for GMgC which showed presence in both the cropping years. Interestingly out of these, three of them were located on 6A, at 597.8 Mb i.e. Tdurum_contig28802_213, wsnp_Ex_c34597_42879693 and RFL_Contig6053_3082, while the fourth one, wsnp_Ex_c34597_42879718 was found on 6B at 597.8 Mb.

2.4.3. Multi-Effect MTA Locus for Ca and Mg

A multi-effect MTA locus, represented by wsnp_Ex_c34597_42879718 on chromosome 6B at 597.8 Mb and was found to be responsible for controlling both the traits simultaneously, and was consistent for both the consecutive crop seasons.

3. Discussion

One of the biggest global health concerns which needs immediate attention is micronutrient malnutrition. High throughput genotyping technologies coupled with statistics provide us vast information and better understanding of genomes. Genome-wide association mapping is an effective strategy to facilitate the identification of genes that regulate traits of interest, whose efficiency is dependent on the genetic diversity within the germplasm that are utilized as association panels. Since genome sequencing and genotyping technologies have advanced so quickly, GWAS has been widely used in wheat and many other crops [23,24]. Gene discovery and cloning have been made easier by the wheat reference genome and a wealth of transcriptomic resources [25,26]. Additionally, this has made it much easier to investigate the QTLs for traits related to wheat yield and quality. Since wheat is one of the most important crop in the world, breeding it to accumulate more nutrients will help reduce nutrient deficiencies. For fulfilling such objectives GWAS is used to identify genomic regions associated with micronutrient trait. Pyramiding of different micronutrients related markers into the high yield genotypes can remove micronutrient deficiency in population. The availability of multiple micronutrient related QTLs/MTAs in wheat germplasm utilized from different countries has also been studied previously [28,29]. Preliminary work on WAMI population was undertaken for investigating these advanced wheat lines for SNPs associated with complex traits without the confounding effects of phenology and plant height [27]. Various panels for investigating QTLs associated with grain Ca have been undertaken earlier, one of which utilized a European wheat diversity panel of 353 varieties with 90k and 35k SNP markers [28]. Moreover, for the mitigation of false positive association, it is important to consider population structure in GWAS [30].

In present investigation, for deciphering the novel GCaC and GMgC associated MTAs, the crop grown in the years 2021-22 and 2022-23 were undertaken, utilizing diverse wheat panel of WAMI. The study decoded involvement of multiple QTLs that contributed for the desired traits and highlighted a total of 86 MTAs (49 for GCaC and 37 for GMgC). Many MTAs identified in our studies have not been reported earlier and therefore, they are potentially novel MTAs controlling GCaC and GMgC. It may be due to the different origin of the panels/populations and different methods used for their detection. Here, chromosome 1A, 1B, 1D, 2A,2B, 2D, 3A, 3B, 3D, 5A, 6B, 7A and 7B have been found to carry QTLs for GCaC in wheat as reported earlier [28,31,32,33,34,35,36,37,38,39], indicating potential role of these chromosomes in different populations for GCaC. The QTLs wsnp_Ex_c2718_5038582, Kukri_c11327_977, and RAC875_c9984_1003 have been found associated with GCaC earlier [24], are consistent with the QTLs detected in the present investigation for both the cropping seasons, making them the most significant loci for GCaC and their accumulation. QTLs Excalibur_c23906_303 and wsnp_Ra_c193_406396 were found to be in the proximity of another QTL reported [40] on 1D for GCaC under hydroponics condition indicating their utility as a potential marker under different conditions.

For GMgC, chromosome 1B, 2A, 2B, 3A, 4A, 4B, 4D, 5B, 6A, 6B, 7A, 7B and 7D were found to carry QTLs for the said trait as found involved in previous studies [32,38,41,42,43,44,45]. The QTLs Tdurum_contig28802_213, wsnp_Ex_c34597_42879693 and RFL_Contig6053_3082 associated with GMgC and reported in present investigation for both the cropping seasons have been identified [22] in a study but under different environmental conditions and thus promises to be genuine for GMgC accumulation. QTL wsnp_Ex_c34597_42879718 located on 6B at 597.8 Mb had been associated with GCaC [22], was also found to be present in both the cropping season and thus promises be a genuine locus for GCaC.

4. Materials and Methods

4.1. Genetic Material and Experimental Conditions

The grain micronutrient concentrations (Ca and Mg) were investigated using spring wheat lines from the Wheat Association Mapping Initiative (WAMI) population from the International Maize and Wheat Improvement Center (CIMMYT) [27]. It consisted of 105 elite and genetically diverse wheat lines, distribution of which was done via International Wheat Improvement Network (IWIN) by CIMMYT. Due to presence of narrow range in variation of days to heading and plant height, this is suitable for gene discovery with no confounding effects of plant height and phenology. The crop was grown for two consecutive cropping seasons in the year 2021-22 and 2022-2023 at the research farm of Eternal University, Baru Sahib, Himachal Pradesh, India using three replications. Each block of a 2 m row with a 0.10 m plant spacing was designated to represent a genotype. Fertilizers used was 120 kg N, 40 kg K₂O and 60 kg P₂O₅, with thorough use of all fertilizers at the time of sowing, except nitrogen. Dose of N was used in three parts: half at sowing, one-fourth on initial irrigation i.e. 21 days after sowing, and rest one-fourth during second irrigation i.e. 45 days after sowing. Distance in plantation was 20 cm and 5 cm between the plants.

4.2. DNA Extraction and Genotyping

Extraction of DNA, genotyping the samples and data processing were followed as previously [46]. Further, genotyping was carried out by USDA-ARS Small Grain Genotyping Center, Fargo, using Illumina 90K Infinium iSelect assay (Illumina Inc., San Diego, CA, USA) [47]. Process of SNP calling used the default clustering algorithm integrated in Genome Studio v2011.1 (Illumina Inc., San Diego, CA, USA), which resulted in identifying total 26,814 bi-allelic SNPs [46,48]. For upholding standards in data quality, SNPs characterized with minor allele frequency (MAF) lower than0.05 were eliminated from the analysis, including monomorphic and low-quality SNPs. Such thorough filtration resulted in retention of around 21,132 polymorphic SNPs, which were utilized for GWAS in our study [46,49].

4.3. Elemental Analysis for GCaC and GMgC

Approximately 0.5 g of homogenised seed material was measured into a digestion tube and subjected to digestion with suprapure nitric acid with a microwave digester. The digestate was filtered using Whatman® filter paper no. 42 and diluted to a volume of 10 ml using ultrapure water in a calibrated volumetric flask.

Concentrations of analytes were quantified utilising inductively coupled plasma–optical emission spectroscopy (ICP-OES; Manufacturer: Perkin Elmer, Model: 7300DV). Calibration was conducted using blanks and five matrix-matched standards. The calibration curve was constructed by linear regression, requiring a minimum correlation coefficient of 0.995. GCaC and GMgC were examined in radial mode, using the standard approach (approach 984.27).

4.3. Statistical Analyses

Analysis of Pearson’s correlation coefficients was conducted utilizing Agricolae (version 1.2–4) package of the R (version 4.0.3) software [50]. Components of variance were studied by restricted maximum likelihood (REML) method implemented in META-R software v6.0.4 [51]. Parameters of covariance were used as component unbiased genetic variance for broad-sense heritability (H2) as Alemu et al. [9]. Data of BLUP was extracted utilizing ‘lme40 package in R [52]. Statistical analysis like, mean, coefficient of variation (CV), and standard deviation were studied utilizing SPSS v. 17.0 (SPSS Inc., Chicago, IL, USA, 2008).

4.4. Population Structure, Kinship Matrix, and Principal Components Analyses (PCA)

Population structure matrix or Q matrix were modelled utilizing PCA for the genotypic data of total 21,132 high-quality SNPs, as earlier [15]. The Kinship matrix (relatedness or K matrix) was evaluated using R (version 4.0.3) software, utilizing the parameters given by VanRaden and Yin et al. [53,54]. Utilizing Bayesian information criterion (BIC), the optimum numbers of PCAs were identified [55]. Scatter plot was formed using first two principal components, which demonstrated the genotype distribution.

4.5. Genome-Wide Association Analyses

GWAS was performed utilizing 21,132 high quality SNPs available from the CIMMYT, Mexico website (https://data.cimmyt.org/dataset.xhtml?persistentId=hdl:11529/10714; accessed on 22 December, 2023). Among the SNP markers, pairwise squared allele-frequency correlations (r2) were determined through TASSEL software (Trait Analysis by Association, Evolution, and Linkage) ver. 5.2.9, with sliding window size of 100. Plot between r2 values and genetic distance (cM) was made for assessing LD between the loci. Forming a smoothing spline regression line at the genome level, LD decay curve utilizing Hill and Weir method [19] in the R environment was achieved with script utilized by Marroni et al. [56]. Following the earlier studies for identification of significant MTAs, significant threshold of p < 0.001 (-log10(P) > 3.0) was considered [9,57,58]. In present study, this QTL has also been identified for GMgC for both the cropping seasons proving its non-redundancy and its ability for detecting multiple micronutrient content in wheat. Breeding for higher GCaC and GMgC content, attention may be paid to these MTAs. These identified MTAs provide valuable insights into the molecular mechanisms corresponding to grain micronutrients in wheat. However, it is prerequisite to consider the environmental factors while making breeding decisions, with proper validation using independent populations which plays crucial role in the performance of the traits.

5. Conclusions

The MTAs so reported in this study were found significant in both the cropping seasons and the non-redundancy was further confirmed by the two GWAS models. This makes them robust enough to be used for breeding for higher GCaC and GMgC and therefore attention may be paid to these MTAs. These identified MTAs provide valuable insights into the molecular mechanisms corresponding to grain micronutrients in wheat. However, it is prerequisite to consider the environmental factors while making breeding decisions, with proper validation using independent populations which plays crucial role in the performance of the traits.