Multi-Environment Quantitative Trait Loci Mapping for Grain Iron and Zinc Content Using Bi-parental Recombinant Inbred Line Mapping Population in Pearl Millet

Pearl millet is a climate-resilient, nutritious crop with low input requirements that could provide economic returns in marginal agro-ecologies. In this study, we report quantitative trait loci (QTLs) for iron (Fe) and zinc (Zn) content from three distinct production environments. We generated a genetic linkage map using 210 F6 recombinant inbred line (RIL) population derived from the (PPMI 683 × PPMI 627) cross using genome-wide simple sequence repeats (SSRs). The molecular linkage map (seven linkage groups) of 151 loci was 3,273.1 cM length (Kosambi). The content of grain Fe in the RIL population ranged between 36 and 114 mg/Kg, and that of Zn from 20 to 106 mg/Kg across the 3 years (2014–2016) at over the three locations (Delhi, Dharwad, and Jodhpur). QTL analysis revealed a total of 22 QTLs for grain Fe and Zn, of which 14 were for Fe and eight were for Zn on three consecutive years at all locations. The observed phenotypic variance (R2) explained by different QTLs for grain Fe and Zn content ranged from 2.85 (QGFe.E3.2014–2016_Q3) to 19.66% (QGFe.E1.2014–2016_Q3) and from 2.93 (QGZn.E3.2014–2016_Q3) to 25. 95% (QGZn.E1.2014–2016_Q1), respectively. Two constitutive expressing QTLs for both Fe and Zn co-mapped in this population, one on LG 2 and second one on LG 3. Inside the QTLs candidate genes such as Ferritin gene, Al3+ Transporter, K+ Transporters, Zn2+ transporters and Mg2+ transporters were identified using bioinformatics approaches. The identified QTLs and candidate genes could be useful in pearl millet population improvement programs, seed, restorer parents, and marker-assisted selection programs.


INTRODUCTION
Pearl millet [Pennisetum glaucum (L). R. Br.] is a climate-resilient and health-promoting nutritious crop of the semi-arid tropics of Africa and Asia (Satyavathi et al., 2015). It is a staple food crop for approximately 90 million people that live in Asia and Africa and practice low-input subsistence farming and livestock production systems (Gangashetty et al., 2016). Pearl millet is an ideal crop in the context of global climate change due to its inherent tolerance to heat stress, salinity, and drought (Shivhare and Lata, 2016). Deficiencies or imbalanced intakes of energy or nutrients, particularly vitamins and minerals, cause a number of dysfunctions and diseases in humans, that are collectively referred to as "micronutrient malnutrition or hidden hunger." Currently more than 2 billion people in the world are in risk of micronutrient malnutrition (Athar et al., 2020). Iron (Fe) and Zinc (Zn) deficiency in humans is caused by lack of purchasing power for highly nutritious food, reduced dietary intake and less bioavailability and bioutilization, especially among resource poor women, infants, and children in the developing countries (Shariatipour and Heidari, 2020). Nutrition supplementation, dietary diversification, and food fortification are among the targeted strategies available to address micronutrient deficiencies. Among, the above strategies, improving the nutrient profile of agricultural crops through genetic means is more pleasing as it is a one-time investment and becomes self-sustaining over a long period. Breeding of pearl millet varieties/hybrids with improved nutritional quality is one of the priority areas for providing nutritional security in the developing countries. Many breeders (Govindaraj et al., 2011(Govindaraj et al., , 2016Rai et al., 2015;Satyavathi et al., 2015;Kumar et al., 2016) have identified huge variability for the trait and utilizing them in development of a handful of varieties and hybrids rich in micronutrients, which now playing a major role in nutritional security in dry lands. However, hybrid released so far in India is having as mean Fe content at par or less than HarvestPlus target level of 77 mg kg −1 (HarvestPlus, 2014). Underutilization of germplasm material with higher nutritional value may attribute due to poor agronomic qualities. Hence there is a great need to make constant effort to breed more and more parental lines which can combine both for yield and micronutrient content along with resistance/tolerance to drought, downy mildew and blast diseases, good keeping quality, etc., and showing a stable performance under different agrogeographical regions.
Genetic maps provide important information for detailed genetic analysis of quantitative traits and have been proven to be significant method for crop improvement (Mohan et al., 1997;Doerge, 2002). Identifying molecular markers which are closely associated with the nutritional trait such as Fe and Zn enable rapid introgression of such traits into elite backgrounds. Among the different marker systems, simple sequence repeat (SSRs) is considered as a powerful and practically useful marker system for marker-assisted breeding (MAB) because of its cheap and more user-friendly along with co-dominant nature, abundance, reproducibility, and variability (Wang et al., 2018). In addition, the size and type of the mapping population determine the gene effect for an economically important trait (Kumar et al., 2018). Recombinant Inbred Lines (RILs) are the best choice among different types of mapping population as they are homozygous as well as have undergone multiple cycles of recombination which helps in mapping of tightly linked markers. Further use of codominant markers along with large RIL population helps in the construction accurate high-resolution map.
There are so many examples that SSR markers being employed for marker assisted biofortification of crop varieties such as QPM and pro-vitamin rich maize hybrids (Gupta et al., 2009;Muthusamy et al., 2014) and development of high yielding basmati varieties with disease resistance (Singh et al., 2012).
There are two different strategies that the plant can able to uptake iron from the rhizosphere: Strategy I and Strategy II (Connorton et al., 2017). Most of the cereals like rice, maize, and wheat follow Strategy II through secretion of phytosiderophores (PSs) (Roberts et al., 2004). TOM1 (Transporter of Mugineic acid family phytosiderophores)/ZIFL4 is a member of major facilitator superfamily (MFS) which exports Fe 3+ -PS chelates in rice and barley (Nozoye et al., 2011). In addition to Strategy II, rice and barley both have a functional Iron-Regulated Transporter 1 (IRT1) homolog that allows direct Fe 2+ uptake from the rhizosphere. Metal transporters, such as the P1B-ATPase family, zinc-regulated transporter (ZRT), iron-regulated transporter (IRT)-like protein (ZIP), natural resistance-associated macrophage protein (NRAMP) family, and cation diffusion facilitator (CDF) family, have been reported in plants, where these two metals play critical roles (Colangelo and Guerinot, 2006). Several genes/ gene families involved in Fe and Zn homeostasis such as yellow stripe-like (YSL), ZRT, ZIP, NRAMP, nicotianamine synthase (NAS), nicotianamine aminotransferase (NAAT), heavy metal ATPases (HMA), zinc-induced facilitator-like (ZIFL), metal tolerant protein (MTP), and others have been characterized in cereal crops such as rice (Anuradha et al., 2012), sorghum (Kotla et al., 2019). In pearl millet, PglYSL, PglZIP, and PglNRAMP and a single gene for PglFER and PglNAS were identified (Mahendrakar et al., 2020). Hence, the aim of this study was to identify the QTL associated with grain iron and zinc content through construct of linkage map using SSRs markers in a large-sized RIL mapping population.

Plant Materials
The experiment materials comprised of 210 F 6:7 RILs mapping population derived from an intra specific cross between PPMI 683 (P 1 , Female) and PPMI 627 (P 2 , Male) made during rainy season, 2009. The RIL population advanced by single seed descent method from F 2 to F 6 from single plants selected during F 2 , and was segregating for grain Fe and Zn contents, respectively. Parent (P 1 ), i.e., PPMI 683 having high grain iron and zinc content and parent (P 2 ), i.e., PPMI 627 having low grain iron and zinc content with checks (ICMB 98222 for high grain Fe and Zn content, ICMB 92111 for low grain Fe and Zn content) ( Table 1). Both checks were planted after every 20 th row of the population.
The average grain yield of five randomly selected plants was taken as grain yield per plant and weight of randomly counted 1,000 seeds in each genotype was taken as 1,000 Seed Weight. This study is mainly focused on grain Fe and Zn content. Hence, only the two most important traits viz., seed yield per plant and 1,000 seed weight were considered for the study.

Environments
A total of 210 RILs along with both parents and checks were evaluated in an alpha lattice design with two replicates (Patterson and Williams, 1976), under field conditions at three different geographical areas, showing all three pearl millet growing agro-climatic zones of India, (i) E1: New Delhi (28 • 382 N,77 • 802 E) representing Zone A with annual rainfall of more than 400 mm, (ii) E2: ICAR-IARI Regional Centre farm, Dharwad (15 • 212 N, 75 • 052 E) from zone B (covering the southern peninsular India), and (iii) E3: Agricultural Research Station, All India coordinated pearl millet project, Mandor,Jodhpur (26 • 252 N,72 • 992 E) falling in zone A 1 with annual rainfall <400 mm. The trials were conducted at all locations during three consecutive (south-west monsoon season, 2014, 2015, and 2016) growing seasons. Each RIL was sown on a four-rows plot of 4 m long with an inter-row spacing of 65 cm and plant-to-plant spacing of 10 cm. For a normal, healthy crop, standard pearl millet cultivation practices were strictly followed. Table  2 shows the weather and soil parameters at each location throughout the year. All climatic parameters except rainfall are presented as means over the crop-growing season (June-October). Rainfall is cumulative rainfall received during the crop-growing season. The relative humidity was calculated as the average of measurements taken in every morning and afternoon on each day.

Sample Preparation and Extraction of Mineral Elements
At physiological maturity, panicles from five plants per plot were randomly selected from the central two rows and harvested separately and threshed using wooden mallet. Seeds from individual plants within the plot were mixed and washed with deionized double-distilled water to remove debris and other possible contaminations in order to obtain a representative sample for trait measurements. The grain samples were sealed in individual paper bags and placed in an oven at 75 • C overnight. About 15 g of air-dried sample from each genotype was further used. Mineral (Fe and Zn) content in the grains during initial year (2014) was measured with EDXRF instrument (a non-destructive, bench-top, energy-dispersive X-ray fluorescence spectrometry, Oxford Instruments X-Supreme 8000) at all three locations (Anuradha et al., 2017). Later, in 2015-2016, a more precise method, di-acid digestion of ground floor samples was used (Dhyan et al., 2005), followed by readings taken on an atomic absorption spectrometer (AAS, ZEE nit 700 tech Analytikjena). To determine whether higher grain mineral content is due to contamination with dust, soil, or is inherent,

Statistical Analysis
Descriptive statistics and Analysis of Variance (ANOVA) was carried out using PBTools v1.4. Histograms and correlations between pairs of traits were estimated through Pearson correlation co-efficient using R software. The model used for ANOVA was: Where µ is the grand mean; yi is the fixed effect of year i; ej is the fixed effect of location j; ye ij is the fixed effect of interaction between year i and location j; gm is the random effect of genotype m and is ∼NID (0, σ 2 g); r ijk is the random effect of replication in location j and year i and is ∼NID (0, σ 2 r); b ijkl is the random effect of block l nested with replication k in location j and year i and is ∼NID (0, σ 2 b ); (yg) im is the random effect of the interaction between genotype m and year i and is ∼NID (0, σ 2 yg ); (eg) jm is the random effect of the interaction between accession m in location j and ∼NID (0, 0, σ 2 eg ); (yeg) ijm is the random effect of the interaction effect of the genotype m in year i and location j and ∼NID (0, 0, σ 2 yg ); and ε ijklm is the random residual effect and ∼NID (0, 0, σ 2 ε ).
Analysis of variance was also conducted using data from each environment for grain Fe and Zn content.
Heritability (H 2 , broad sense) at individual environment was estimated from analysis of variance. The formula used was- Whereas Heritability (H 2 ) estimates across environments were estimated by the formula- Where r, y, l denotes the number of replicates, years and environments, respectively. The REML model also produced best linear unbiased predictors (BLUPs) of each genotype, thereby adjusting the influence of the neighboring rows. These BLUPs were used for downstream analysis.

Genomic DNA Extraction and Genotyping
The genomic DNA was extracted from fresh and young leaf tissue of each RIL and the parents by using modified cetyltrimethylammonium bromide (CTAB) method (Aboul-Maaty and Oraby, 2019). Polymorphism was assessed using 372 SSR markers, EST-SSRs (ICMP and IPES series), gSSRs (PSMP and CTM series), and gene-based markers. SSRs were amplified in 10 µl PCR reaction mixture containing 25-30 ng genomic DNA, 100 µM of each of the four dNTPs, 1.5 mM MgCl 2 , 5 pmol/µl of forward and reverse primer, 1 U of Taq DNA polymerase (Bangloregenei, Merck) and 1× PCR buffer in Biorad My Cycler Thermal Cycler, United States. PCR conditions were used as described by Anuradha et al. (2017). The amplified product was checked using 1.2% agarose gel and PCR products were separated by 10% polyacrylamide gel electrophoresis at 120 V (C.B.S. Scientific vertical electrophoresis unit). The SSR marker alleles were scored manually by comparing the position of bands with that of 100 bp ladder from bottom to top (smallest to largest). Based on the position, the allele size in base pairs (bp) was recorded for all polymorphic SSRs in all 210 RILs studied. The parent P 1 (PPMI 683) was scored as "A" and parent P 2 (PPMI 627) was scored as "B" and both type band present on single locus were designed as "H". Diffused or ambiguous bands were recorded as missing bands with "0" notation. The chi-square method was conducted to test whether the inherited alleles of the mapping population were consistent with Mendelian segregation ratios. The segregation ratio across the mapping population was tested against a 1:1 ratio. The segregation pattern that did not fit either ratio (p < 0.05) were treated as distorted.

Linkage Analysis and QTL Mapping
The linkage map was constructed using 151 polymorphic SSR markers. The markers that were not taken for the linkage analysis were those with: a minor allele frequency of less than 5%, more than 20% missing data, and those that were monomorphic between the parents. The MAPMAKER/EXP3.0 software was used to develop a genetic linkage map (Lander et al., 1987). The maximum recombination fraction (θ) was set at 0.49 and the critical logarithm of odds (LOD) score for the test of marker pair independence was set at 3. Markers with high segregation distortion (p ≤ 0.01) were removed in a preliminary analysis of the data. Subsequently, the linkage analysis of genotyping data for all the markers including those with distorted segregation was used. Using MapChart software, a graphical map was generated (Voorrips, 2002). The consensus map developed by Supriya et al. (2011) and Rajaram et al. (2013) was used as a reference map in the current study for grouping markers and constructing the map. QTL analysis was carried out using the Inclusive Composite Interval Mapping method (ICIM). The ICIM approach uses a strategy in which a stepwise regression is firstly performed, so markers with significant effect on QTL are selected. ICIM-ADD method was used with the multi-environmental model builtin QTL IciMapping (Li et al., 2015). Significant LOD thresholds were set at 5% tail of the null distribution in a 1,000 permutations test (Luciano Da Costa et al., 2012).

Quantitative Trait Loci Validation
All QTL intervals identified for grain Fe and Zn content were analyzed in-silico with reference genome, Tift23D 2 B 1 -P 1 -P 5 (PROJECT ID: PRJNA294988.) in order to navigate the presence of candidate gene associated with Fe and Zn content in grain. The physical position of flanked markers of identified QTLs Frontiers in Plant Science | www.frontiersin.org were obtained from NCBI-BLASTn 1 . The sequence similarities searched against related cereals (Setaria italica) by Ab initio gene prediction method by using Fgenesh software for the identification of potential protein coding regions (Salamov and Solovyev, 2000) 2 . For each predicted exon similarity searches were run using the BLAST program on protein and EST databases (Altschul et al., 1997).

Phenotypic Evaluations
The mean performance of parents and range of trait values of RILs during three consecutive years (2014, 2015, and 2016) at E1, E2, and E3 are summarized in Table 3. Both traits (Fe and Zn) were approximately normally distributed when analyzed across the years at their respective locations (Figure 1). The levels of grain Fe in RILs ranged from 43.27 to 108.33 mg kg −1 , 33.34-118.52 and 32.88-115.47 mg kg −1 , while the range of grain Zn content was 22.4-94, 20.97-120.87, and 17.33-105.17 mg kg −1 across the years at E1, E2, and E3, respectively. The mean BLUP value of grain Fe and Zn content in both parents and RILs were higher across the described years at E1 when compared to E2 and E3, whereas mean BLUPs of grain Fe and Zn for both parents was lowest at E3 across the years. For grain Fe and Zn studied in different environments, both parents and RILs showed a wider range of variation.
The seed yield per plant and 1,000 seed weight of the 210 RILs along with both parents varied greatly. At E1 during three consecutive years (2014-2016), 1,000 seed weight ranged from 6.75 to 12.17 g with a mean of 9.05 g, and seed yield per plant ranged from 14.09 to 36.21 g with the mean of 24.05 g. And at E2 from 2014 to 2016, 1,000 seed weight ranged from 6.77 to 12.33 g with a mean of 9.09 g, and seed yield per plant ranged from 19.05 to 35.37 g with a mean of 26.81 g. From 2014 to 2016, 1,000 seed weights in E3 ranged from 6.48 to 11.41 g, with an average of 8.38 g, and seed yield per plant was 18.20-39.8 g, with an average of 26.75 g.
The ANOVA indicated a highly significant G × Y interaction [p < 0.0001, F (66.45) = 1.112 for Fe] and [p < 0.0001, F (62.43) = 1.112 for Zn]. The pooled mean, range, coefficient of variation (CV), genotypic variance (σ 2 g ), residual variance, broad sense heritability (H 2 ), least square difference (LSD) and G × E variance in RILs evaluated during 2014-2016 cropping season at E1, E2, and E3 has shown in Table 4. In the RIL population, the genotypic components of variance (σ 2 g ) for both traits were significant, and the operational heritability (H 2 ) estimates ranged from 0.83 (Zn) to 0.88 (Fe) ( Table 4). The high heritability shows that much of the phenotypic variance in the population is genetically controlled and QTL can be mapped with a high degree of reliability. Broad-sense heritability was usually high enough to allow effective QTL mapping, indicating that the traits studied have moderately high to high proportions of genetic 1 https://blast.ncbi.nlm.nih.gov/ 2 http://www.softberry.com/berry.phtml?topic=fgenesh&group=programs& subgroup=gfind  Significance at **p < 0.01.
variance. Correlation coefficients (r) for grain Fe and Zn content were found statistically significant (p < 0.001) for each year of evaluations for each location. Homogeneity of variance tests based on Bartlett's test indicated homogeneous error variance for both the traits which allowed for a combined analysis across the years for all locations. When performing a combined ANOVA, the year was considered as a random and genotypes were considered as fixed effects. Combined ANOVA showed significant Genotype × Environment Interaction (GEI) (p < 0.001), exhibiting the influence of changes in environment on the grain micronutrient contents of genotypes. It was observed that Genotypes (G) and GEI effects were highly significant (p < 0.01) and the contribution of year is very less for both the traits ( Table 5).
Correlation of Grain Fe, Zn, 1,000 Seed Weight and Yield Grain Fe and Zn content were highly (0.711) and significantly (p < 0.01) correlated ( Table 6). Thousand seed weight was also observed to be positively correlated with seed yield per plant (0.579, p < 0.01). For successful incorporation of high levels of both grain Fe and Zn content, association between these micronutrients and their association with thousand seed weight and yield is important to plan a successful breeding approach so as to develop new improved lines or varieties having higher grain yield along with enhanced levels of grain micronutrients.

Linkage Map Construction
For the polymorphism survey of both parental lines, 366 SSR primer pairs (IPES, PSMP, and ICMP series) and six gene-based primers were used. Of them, 151 (40.59%) SSRs detected polymorphism between the two parents. Polymorphic SSR markers included 94 IPES, 42 PSMP, and 11 ICMP series primer pairs and four newly synthesized primer pair named as Ppmsb (Pusa pearl millet SSR biofortification) (Figure 2). These polymorphic primers were further used for linkage map construction with LOD score 3.0.
The total length of the linkage map was 3273.1cM (Kosambi). The individual linkage group ranged from 910 cM for LG 2 with the highest number of markers (31) to 135.3 cM for LG 3 with 14 markers, while the least number of markers were assigned on LG4 (10) followed by LG 3 (14) and LG 5 (14). The average interval size was 8.79 cM. The average LG length was 467.5 cM with an average of 21.57 loci. The average adjacentmarker interval lengths ranged from 3.07 cM (LG 3) to 14.91 cM (LG 2) followed by 11.34 cM (LG 1). Map distance between adjacent markers varied from 0.9 to 49.2 cM (LG 1) and 8.33% of the intervals (31 out of 372 intervals) were smaller than 10 cM. Details of all linkage group, total markers, skewed marker and their percentage, total map length (cM) and average distance (cM) were given in Table 7.

Quantitative Trait Loci Analysis
The multi-environmental QTL analysis exhibited the association of various genomic regions with grain Fe and Zn content at all the three locations during three consecutive years, 2014, 2015, and 2016 (  Fe, SY, TSW, and Zn are Grain iron content in mg/kg, Seed yield per plant in g, thousand seed weight in g, and Grain zinc content in mg/kg. **Significance at p < 0.01. QGFe.E3.2014-2016_Q4 for Jodhpur location (E3) (Q = QTL; G = grain; Fe = iron). In addition, two more QTLs for GFe were also found at E3, flanked by IPES0101-IPES0139 and IPES0027-IPES0236 and having very less PVE of 2.92 and 2.85% on LG 1 and LG 2, respectively. The genome wide significant threshold was LOD = 3.5 for E1, E2, and E3. For GZn, one common QTL was found in E1, E2, and E3 on LG 3 at significant threshold LOD value of 3.5 which was flanked by IPES0142-IPES0180 with total phenotypic variance of 14.39, 16.07, and 22.26% at E1, E2, and E3, respectively. And one more common QTL was also found at E1 and E2 on LG 2, flanked by Ppmsb001-Ppmsb002 having total PVE of 25.95 and 21.88%, respectively. Whereas, at E3, one QTL with flanked markers Ppmsb002-IPES0181was mapped on LG 2 with comparatively low total PVE of 11.18%. Two more un-repetitive QTLs between IPES0166-PSMP2249 (PVE 6.2%) on LG 3 and IPES0141-IPES0151 (PVE 2.93%) on LG 6 were detected in E1 and E3, respectively.
The proportion of the total phenotypic variance (PVE) explained by the QTLs for GFe ranged from 14.99 to 19.66%, 16.26-18.83, and 2.85-16.85% in E1, E2, and E3, respectively (   LG 5 over the locations, which was also displayed 3.6 additive effect of QTL (Table 8). Whereas, the QTL for GZn having largest PVE (25.95%) was QGZn.E1.2014-2016_Q1 on LG 2 over the locations (Table 8). Similarly, one more QTL for grain Zn (QGZn.E2.2014-2016_Q1) was also found with high PVE 21.88% on LG 2 at E2. The additive effect was detected for all the QTLs on which the alleles from PPMI 683 (donor parent) had positive additive effect.  (Table 9) were selected and their DNA was used for amplification with ten consistently associated primers, IPES0142, IPES0180, IPES 0160, IPES0093, IPES0206, IPES0015, IPES0166, PSMP2249, Ppmsb001, and Ppmsb002. The amplicons of these genotypes were obtained by PCR amplification of these primers. The similar regions of reference sequence obtained were thus compared with the amplified sequences of high and low grain iron of expressing pearl millet genotypes. These markers were compared with reference genome Tift 23 D 2 B 1.

In-silico Identification of Candidate Gene
Candidate genes identified for each QTL was shown in Table 10.
In the study, we identified ten metal transporter genes. Out of the ten genes seven were identified in the genome of Seteria italica whereas remaining three genes were identified in the genome of Zea mays. Among the genes from the Seteria italica,

Variance Components and Heritability
Pearl millet has a great variation of mineral and agronomic traits. For both traits, the parents of this RIL population exhibited  statistically significant divergent phenotypes. The BLUP means of RILs were significant for Fe and Zn content. A very wide range of variation was detected among the population for grain Fe and Zn content ( Table 3).
Wider range of variation existed in the present mapping population for grain yield per plant (18.89-36.70 g) and thousand seed weight (6.84-11.79 g). Wide range of variation was earlier reported by many workers in pearl millet for yield and micronutrient traits (Govindaraj et al., 2011;Sankar et al., 2013;Vinodhana et al., 2013;Bashir et al., 2014;Sabiel et al., 2014;Anuradha et al., 2017). Greater amount of variation in the present population show that there is a lot of scope for selection of better lines to increase the grain yield in pearl millet.
The operational heritability (H 2 ) estimates were very high and the genotypic components of variance (σ 2 g) for both traits were significant. Similar study was also done by Kumar et al. (2016Kumar et al. ( , 2018. The trait heritability information is helpful to choose the selection procedure to be followed to improve the trait. Broad-sense heritability was high enough to allow effective QTL mapping, with moderately high to high proportions of genetic variance for the traits investigated. Higher estimates of heritability with genetic advance as a percent of mean were observed, suggesting the presence of additive gene action efficacy of selection for both traits. The traits which expressed high heritability and low genetic advance showed non additive gene action, hence heterosis breeding would be recommended for these traits. Sumathi et al. (2010) and Govindaraj et al. (2011) were also found similar result.

Frequency Distributions
The frequency distributions of the overall means for grain iron and zinc content in the mapping population are represented by histogram (Figure 1). Likewise, previously Kumar et al. (2016Kumar et al. ( , 2018 were also concluded the normal frequency distribution for both traits in pearl millet. The normal Frequency distribution for both traits was also observed in other crop like wheat (Crespo-Herrera et al., 2017). Both traits in the RIL population have phenotypic normal distributions and transgressive segregations, indicating polygenic inheritance. The presence of transgressive segregation also suggests genetic recombination (Falconer, 1996), which means that both favorable and unfavorable alleles for the traits are dispersed between the parents.

Correlation Analysis
The understanding of the association between nutritional and agronomic traits will support the breeders to select suitable selection/breeding program for genetic improvement of associated traits such as mineral content. Several crops have been studied to determine the correlation between grain Fe and Zn contents. Iron and zinc didn't record any significant association with grain yield and thousand seed weight albeit they in turn are correlated. In many previous studies there was negative association of Fe and Zn with yield (Kanatti et al., 2014), which hampers the effective selection of all these traits. This study is in accordance with Anuradha et al. (2017), where they didn't notice any associations. Hence, there is a possibility for selection of high Fe and Zn content along with high yield. In current study, positive correlation between iron and zinc was observed as also reported in many previous studies (Rai et al., 2012;Govindaraj et al., 2013;Kanatti et al., 2014;Kumar et al., 2016Kumar et al., , 2018Anuradha et al., 2017Anuradha et al., , 2018Velu et al., 2017). This could be due to common molecular mechanisms that regulate mineral uptake and metabolism in grains, or to common transporters controlling mineral movement within plants (Vreugdenhil et al., 2004;Ghandilyan et al., 2006;Magallanes-López et al., 2017). The strong association between these minerals in populations may be due to the co-segregation of genes for these traits. The direction and intensity of the association suggest that cotransferring superior alleles regulating these traits into the genetic backgrounds of elite lines could provide good opportunities for simultaneous genetic improvement of both micronutrients (Allouis et al., 2001). Genotypic correlations were found to be higher for all traits, implying that there is less interaction between the genetic makeup of the traits and environment.

Linkage Map
Linkage mapping is essential for QTL mapping and markerassisted breeding programs. The most preferred mapping population is RIL, but the accuracy of the mapping resolution is dependent on population size. The linkage maps have been predominantly constructed using F 2 populations and up to certain extent RILs in pearl millet (Ferreira et al., 2006). Only few previous studies were done on linkage map and QTLs in pearl millet like Supriya et al. (2011) had constructed the linkage map based on 321 loci (258 DArTs and 63 SSRs) for 140 RILs. Similarly, Kumar et al. (2016Kumar et al. ( , 2018 had also used DArTs and SSRs for construction of linkage map and used a smaller number of SSRs for genotyping of RILs and map development. In our study, we found 151 polymorphic SSRs out of 366 and used for linkage map development for 210 RILs. The mapping population size and number of SSRs was high compared to the Supriya et al. (2011) and Kumar et al. (2016Kumar et al. ( , 2018. The present study recorded lesser polymorphism for SSRs compared to that of Senthilvel et al. (2008) who reported 74%. It may be because the parents used here for the development of RILs are closely related than that of Senthilvel et al. (2008). The distribution of SSRs across the linkage group was uniform in the current study. In terms of SSR marker positioning on the genetic map, all SSRs mapped on the similar linkage group (LG) as previously reported pearl millet maps by various researchers (Yadav et al., 2004;Senthilvel et al., 2008;Supriya et al., 2011;Kumar et al., 2016Kumar et al., , 2018. One hundred fiftyone polymorphic SSR markers were assigned to seven previously established pearl millet linkage groups (Figure 2). The present map had a larger length than the one reported by Yadav et al. (2004); Supriya et al. (2011) and Kumar et al. (2016Kumar et al. ( , 2018. Recombinant inbred lines usually exhibit segregation distortion, because, many recessive lethal genes become homozygous during the RIL development process (Kumar et al., 2016). Segregation distortion was found in 31.57% of the loci on the current linkage map. Whereas Kumar et al. (2018) found 60% segregation distortion of DArts and SSRs. Pollen abortion is most common in pearl millet than abnormalities in female gametes, leading to a comparatively greater loss of male alleles and the subsequent skewness toward the female parent (Kumar et al., 2018). Distortion from expected Mendelian segregation has been observed previously in maize (Wendel et al., 1987;Lu et al., 2002), barley (Graner et al., 1991;Devaux et al., 1995), rice (Causse et al., 1994;Xu et al., 1997), wheat (Blanco et al., 2004;Quarrie et al., 2005), and pearl millet (Supriya et al., 2011). The protogynous nature of pearl millet also contributes to segregation distortion (Liu et al., 1994). According to Cloutier et al. (1997) and Livingstone et al. (1999), residual heterozygosity and inbreeding depression during inbred development may also contribute to segregation distortion. The residual heterozygosity in some RIL may be advantageous at the same time, because deleterious genetic combinations in the form of reduced fitness or lethality can be avoided. The segregation of nearly all loci of LG 3 was distorted in comparison to earlier reports, which may explain part of the increase in LG 3 map length.
The presence of large gaps between the centromere and telomere is a characteristic feature of pearl millet linkage maps. Senthilvel et al. (2008) found a large gap in LG 4. The previously built framework map for the cross ICMB 841-P3 × 863B-P2 had big gaps in LG 2 and LG 7 (Senthilvel et al., 2008;Yadav et al., 2011). This new population also revealed a large gap (>25 cM) in LG 1, LG 2, LG 4, and LG 6, possibly due to extreme recombination localization at the ends of LG. The interval distance between the markers in LG 4 was high suggesting additional polymorphic markers are needed to fill-up the gap. According to many researchers like Devos et al. (2000); Varshney et al. (2006), andSenthilvel et al. (2008), large gaps in the distal regions reflect regions of high recombination, rather than a lack of markers in these regions. However, these linkage groups are still incomplete, and need the addition of new markers that are located on the distal regions of the linkage groups. The number of large gaps has tried to decrease in the present study, although there is still the possibility to map more markers to fill these gaps.

Quantitative Trait Loci for Grain Iron (Fe) and Zinc (Zn) Content and Environment Interaction
The use of QTL analysis to map populations is useful not only for identifying genomic regions associated with traits of interest, but also for using the associated marker information in breeding programs to integrate particular loci in elite germplasm. The application of QTL by environment interaction with composite interval mapping was used in the present study (Li et al., 2015), it was possible to locate several genomic regions associated with GFe and GZn during the three consecutive years over three diverse locations. Our findings support other authors' findings that GFe and GZn are quantitative traits (Kumar et al., 2016;Anuradha et al., 2017). These previous QTL studies have also mapped QTL or genomic region for GFe and GZn in different linkage group of pearl millet, including LG 3 and LG 5, with PVE 13.6 and 10.0%, respectively, for iron content (Kumar et al., 2016). on the other hand, Anuradha et al. (2017), identified different markers on different linkage groups for both high grain iron and zinc content, namely PSMP 2261 (LG 5), IPES0096 (LG 7), IPES 0180 (LG3), and Xsinramp 6 (unmapped). Kumar et al. (2016) reported the co-localization of high grain iron and zinc content alleles/ QTLs in pearl millet. It may happen because starting from uptake to final deposition into grain of iron and zinc may share some common pathways (Grotz and Guerinot, 2006). In present study we found four repetitive QTL on LG 2, LG 3, LG 5, and LG 7 with PVE ranging 16. E2,and E3,respectively,and the largest PVE (19.63%) was displayed by QGFe.E1.2014-2016_Q3 in LG 5. Kumar et al. (2016) has mapped one SSR marker (IPES142) on LG 3 and Anuradha et al. (2017) has mapped IPES180 on LG3. In the present study one QTL was mapped between IPES142-IPES180 on LG 3 which is in accordance with that of Kumar et al. (2016) and Anuradha et al. (2017). Further, Anuradha et al. (2017) mapped a marker ICMP 3092 on LG 7 for Fe. Similarly, we also identified one repetitive QTL between IPES0206-IPES0015 on LG 7 and interestingly the nearest marker to this QTL is ICMP 3092 as per Rajaram et al. (2013). However, the genomic regions of this study and other studies are not comparable since Kumar et al. (2016) used DArT markers while this study used SSR markers. Similarly, the mapping method used by Anuradha et al. (2017) was association mapping while it is QTL in the present study.
The genotype by environment interaction, particularly the cross-over type interactions, has significant implications for crop performance and breeding. The ideal case for pearl millet biofortification is to obtain stable pearl millet genotypes that perform well without cross-over interaction when tested in different environments or years in a particular geographical area. The analysis of the phenotypic data in our study revealed the presence of a significant G × Y interaction. QTL analysis revealed that most of the LOD scores for the additive average effect were higher than the LOD score for the interaction ( Table 7), showing that QTL with higher LOD (Add) are more stable than those with higher LOD (G Y) (Li et al., 2015). The additive effect was detected for all the QTLs for grain Fe and Zn, on which the alleles from donor parent had positive additive effect. Gaoh et al. (2020) found the additive × additive gene effects had the most important effects for grain Fe content, while additive × dominance gene effects were significant for grain Zn content in pearl millet. Similarly, Garcia-Oliveira et al. (2018) and Kumar et al. (2018) reported that epistasis plays a role in regulating grain iron and zinc content in pearl millet and cereals, respectively. Lu et al. (2008) reported 28 genome-wide additive × additive interactions for mineral elements in rice. Anuradha et al. (2012) also reported epistatic interactions between loci on chromosomes 1 and 5 for grain Zn concentration. Jin et al. (2013) reported the additive effect for Zn and partial dominant for Fe in maize.

VALIDATION OF IDENTIFIED QTLS IN A DIVERSE SET OF GERMPLASM
Validation of marker trait associations can reveal their worthiness for further utilization in breeding. Hence, DNA from eleven high, nine low grain iron and four check (two for high and two for low) expressing genotypes were amplified with ten SSR markers which were consistently associated with both the traits. Flanked SSR markers of identified QTLs were compared with the reference Tift 23 D 2 B 1 genome.

In silico Identification of Candidate Gene
The availability of the pearl millet genome sequence has allowed researchers to search for candidate genes involved in Fe and Zn accumulation in grains. The presence of candidate genes was investigated in genomic regions encompassing all QTLs contributing to Fe and Zn in pearl millet. QTLs identified for Fe and Zn content using composite interval mapping were analyzed In-silico for the presence of putative candidate genes. In the study, we identified ten metal related genes within the obtained QTLs.
Foods with higher iron content can be produced using modern genetic and molecular technologies. To increase the amount of iron and zinc in wheat, ferritin is over-expressed . Other mineral concentrations were also affected by ferritin expression, as reported by several studies (Drakakaki et al., 2005;Qu et al., 2005). The promoter selected to control ferritin expression is critical to iron accumulation in specific tissues. The Fe and Zn concentrations in transgenic indica rice grain improved by endosperm-specific glutelin gene promoter of soybean ferritin (Vasconcelos et al., 2003). According to Liu et al. (2016), over-expression of sickle alfalfa ferritin, which is controlled by the seed-storage protein glutelin GluB-1 gene promoter, increases grain Fe and Zn concentrations but also affects mineral homeostasis in transgenic wheat grain. Wheat TaALMT1 (ALMT, for Al-activated Malate Transporter) encoding a malate transporter involved in Al tolerance (Raman et al., 2005). Over expression of TaALMT1 in wheat, barley, and tobacco-cell suspension increases the efflux of Al-activated malate and increases the tolerance to Al stress (Pereira et al., 2010). Root exudation of organic acid such as malate/citrate may play an important role in providing iron to plants (Brown, 1983). Citrate efflux is crucial for iron translocation, and this is mediated by the efflux transporter FRD3 in Arabidopsis (Green and Rogers, 2004) and its orthologue FRDL1 in rice (Yokosho et al., 2016). Barak and Chen (1984) reported a positive effect of K on iron absorption which was associated with acidification of the rhizosphere. Sakaguchi et al. (1999), observed secretion of mugineic acid family phytosiderophores (MAs) from the use of potassium gradient, and content of potassium in barley roots increased with iron-deficiency. Grotz et al. (1998) identified the zinc transporter genes in Arabidopsis thaliana. Yang et al. (2020) studied on activity of zinc transporter (OsZIP9) mediates zinc uptake in rice. Plants acquire Mg from the environment and distribute in the ionic form via Mg 2+ -permeable transporters. In Arabidopsis, MGT6, (plasma membrane-localized magnesium transporter) mediates Mg 2+ uptake under Mg-limited conditions (Yan et al., 2018). The Magnesium/proton exchanger protein is a vacuolar exchanger of protons with cytosolic Mg 2+ and Zn 2+ in Arabidopsis. Hence, it may be possible that Mg and Zn share a common pathway of transport in plants.
Ferritin Gene (Ppmsb001-Ppmsb002) Plant ferritins are important iron storage proteins that serve as both an acceptor and a donor of iron in metabolic processes (Briat et al., 2010). It shares structurally and functionally similarities with animal ferritins. Plant ferritin is generally observed in cells which are inactive in photosynthesis such as seed, root, root nodules, and others. Plant ferritins are nuclear encoded and located primarily in plastids (Zielińska-Dawidziak, 2015). This protein also controls the cellular concentration of transition metals, metal ions apart from iron in the mineral core (including cadmium, beryllium, zinc, aluminum, and lead) (Kumar and Prasad, 1999). Many gene families like FER, ZIP, NRAMP, YSL, and NAS, were involved for mineral translocation from soil to the grain filling site. These gene families have been identified in plants such as Arabidopsis thaliana (Weber et al., 2004), Setaria italica (Alagarasan et al., 2017), and Oryza sativa (Neeraja et al., 2018). In the current report, the ferritin gene was reported by in silico method, which is similar to the Ferritin-like (PglFER1) gene family previously reported by Mahendrakar et al. (2020). With reference to previous studies it is clear that ferritin gene regulates the iron metabolism in plants.

Al 3+ Transporter in Plants (IPES0142-IPES0180)
Aluminum is a potentially phytotoxic metal and Al tolerant plants excrete organic acids such as malate, citrate, and oxalate from root tips, depending on the plant (Ryan et al., 2011). The organic acid exudation is vital for plants tolerating metal and nutritional stress at the root-soil interface. Aluminum interferes the uptake or transport of nutrients such as Ca, B, Fe, Mn, P, Mg, and Cu or K (Keltjens and Tan, 1993;Keltjens, 1995;Lukaszewski and Blevins, 1996;Ślaski et al., 1996). Root elongation and nutrient (B, Fe, Mg, Ca, and P) absorption may also be inhibited by Al or low pH soils. Root exudation of organic acid such as malate/citrate may play a key role in supplying Fe to plants (Brown, 1983). Hence, with the consideration of previous studies, we assumed that the aluminum activated malate transporter gene is responsible to cease the aluminum absorption by plant root and also induce the iron absorption.
K + Transporters in Plants (IPES0142-IPES0180, IPES0157-IPES0093, and IPES0206-IPES0015) Oertli and Opoku (1974), as well as Barak and Chen (1984), identified a positive effect of K on Fe nutrition, which was linked to rhizosphere acidification (excessive K + uptake and consequent release of H + ions by roots to maintain a cation/anion balance). The activity of H+ ATPase pumps located at the plasma membrane is responsible for rhizosphere acidification, which is essential for nutrient acquisition. The liberation of proton is in the favor of uptake of Fe and other macro and micronutrients, especially under the deficiency conditions (Houmani et al., 2015). Hughes et al. (1992) found a significant effect of K + on proton efflux and ferric reduction mechanisms (Strategy I of iron uptaking by roots). They also identified that K + play specific physiological functions in the biosynthetic pathway of mugineic acid production and in the transport of Fe 3+ -phytosiderophore complex (Strategy II). So, we can assume that K + transporter also regulate the iron absorption in plants under iron stress condition.
Zn 2+ Transporters in Plants (XIPES0157-XIPES0093) In present study, we have notified one zinc transporter gene of ZIP superfamily that was earlier reported in pearl millet by Mahendrakar et al. (2020). In higher plants, members of the ZIP family are involved in metal uptake, transport, and accumulation of iron and zinc in plant cells. The first identified zinc transporter genes were ZIP1, ZIP2, and ZIP3 in Arabidopsis thaliana (Grotz et al., 1998). Several ZIPs, including OsIRT1, OsIRT2, OsZIP1, OsZIP3, OsZIP4, OsZIP5, OsZIP7, and OsZIP8, have been identified as being responsible for Zn uptake from soil, translocation within the root and from the root to the shoot, and storage in rice seeds (Yang et al., 2009;Lee et al., 2010).

Mg 2+ Transporters in Plants (IPES0157-IPES0093)
Mg 2+ transporters are used to maintain Mg in plant organs. Among the proteins potentially involved in Mg 2+ transport, Magnesium transporter (MGT) family Mg 2+ transporters plant was the most well-investigated protein (Kobayashi and Tanoi, 2015). Meanwhile, the other transporters of Mg 2+ are also possible such as Non-selective cation channels (NSCCs) (Guo et al., 2003). According to Demidchik and Maathuis (2007), Voltage-independent NSCC (VI-NSCC) reported the uptake of several cations including Mg 2+ , Ca 2+ , and Zn 2+ at the resting membrane potentials. The Arabidopsis Magnesium/proton exchanger (MHX) protein is a vacuolar exchanger of protons with cytosolic Mg 2+ and Zn 2+ . Hence, it may be possible that Mg and Zn share a common pathway of transport in plants.
We have identified one of the gene in ferritin gene family that was already reported in pearl millet via insilico studies and validated by qRT-PCR gene expression (Mahendrakar et al., 2020). The transporter genes like Al, Mg, and K transporter are reporting for the first time in pearl millet. Whereas these transporter gene family were found previously in many different crops such as Al transporters gene in wheat (Li T. et al., 2017) and Rice (Sonoda et al., 2003); Mg transporter gene in Maize (Li et al., 2016 and Arabidopsis (Li H. et al., 2017); K transporter gene in Maize (Zhang et al., 2012) and Barley (Santa-María et al., 1997); and Zn transporter gene in Arabidopsis (Grotz et al., 1998), rice (Ishimaru et al., 2005;Lee et al., 2010;Neeraja et al., 2018) and Maize (Li et al., 2013).

CONCLUSION
Enhancing the level of grain Fe and Zn content among elite cultivars of pearl millet become the priority area in answering the nutritional imbalance among poor inhabitance of the semi-arid tropics, especially women and children. Unravelling the genetic principles involved in inheritance of complex traits like grain micronutrient content in pearl millet via, trait-specific mapping enables rapid genetic gains. Huge genetic variation among RIL population observed under the present study can be further used to develop high yielding micronutrient rich cultivar through heterosis breeding. The two co-localized QTLs on LG 2 and LG 3, for both the minerals found to promising as it could simultaneously transferred to the parental lines of elite pearl millet hybrids through marker-assisted breeding programs.

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

AUTHOR CONTRIBUTIONS
CS, SPS, and AK designed and supervised the overall research and contributed in the preparation of the manuscript. CS and SPS provided technical guidance for creating liaison among multi-environment locations, and edited the manuscript for final submission. TS, SMS, JB, and MM executed the field experiments. TS and SMS performed the phenotyping, carried out statistical analysis, and prepared the manuscript draft. CB and NA reviewed the manuscript. All authors contributed to the article and approved the submitted version.