The study utilized 174 genotypes, comprising a subset of the 3K genome panel (Wang et al., 2018) of 170 genotypes and four check varieties. Initiated with 192 genotypes (188 germplasm lines and four checks), the study could continue only with 174 genotypes, as 14 genotypes suffered establishment issues rendering them unusable in the downstream analyses. The checks were contrasting excess Fe-responsive genotypes, Shahsarang (tolerant), Megha SA2 (tolerant), Pusa 44 (sensitive), and IR64 (sensitive). Shahsarang is a medium-tall (100–110 cm), non-lodging, and medium-duration (120–130 days) indica rice variety with moderate resistance to blast disease and high tolerance to Fe toxicity. Shahsarang performs well in acidic and lowland soils, and its grains are notable for their high Fe and Zn content. Megha SA2, another variety tolerant to excess Fe (Sonu et al., 2024), is a late-maturing cultivar well suited to mid- to low-altitude lowland ecosystems. It produces medium-aromatic, long slender grains with awns and red-colored kernels. Among the sensitive checks, IR64 is an early-maturing, high-yielding variety with comparatively lower grain micronutrient levels, while Pusa 44 is a long-duration, semi-dwarf, high-yielding cultivar known for its sturdy stem. The subset assembly had predominantly Indian origin genotypes comprising cultivars, landraces, and breeding lines ( Supplementary Table S1 ). All materials were sourced from the ICAR-Indian Agricultural Research Institute (IARI), New Delhi.

Field evaluations were conducted at three diverse agro-ecological locations in India. These sites were New Delhi (DEL), Aduthurai (ADT), and Barapani (BPN) having four sites with varying soil Fe levels. At DEL, the experiment was conducted during Kharif 2022–2023 at the experimental plots of the Division of Genetics, ICAR-IARI (28°38′ N; 77°10′ E; 246 m), while at ADT, the field at the IARI Rice Breeding and Genetics Research Centre (RBGRC, Tamil Nadu) was used (11°00′ N; 79°28′ E; 41 m) in Rabi 2022–2023. Located at the experimental farm of the ICAR-Research Complex for North Eastern Region, Umaim, Meghalaya (25°41′ N; 91°54′ E; 980 m), at BPN, there were two different lowland field conditions. Site 1 was with normal Fe (24.6 ppm) (BPN-N) and site 2 was with excess Fe-stressed (72.0 ppm) (BPN-S) conditions. DEL soil was relatively low in Fe (16 ppm) (Rani et al., 2021) (site 3), while ADT had moderate soil Fe content (24.6 ppm) (site 4) (Jeyasingh et al., 2023). The soil Zn content at all of the sites was moderate, with DEL having a DTPA-extractable Zn content of 2.5 ppm, BPN with 1.96 ppm, and ADT with 1.2 ppm. The soil types and pH varied across sites, with DEL having sandy loam soil (pH 8), ADT with alluvial soil (pH 6.28), and BPN having red loam acid soil (pH 5.63 for BPN-N and pH 5.2 for BPN-S).

The field experiments have been laid out in augmented randomized complete block design at all three sites under irrigated transplanted conditions. The nursery was raised on elevated beds for 21 days and transplanted into well-puddled field with a spacing of 20 × 15 cm. The field was divided into four blocks, with 47 genotypes in each block to accommodate the initial panel of 188 genotypes. The check genotypes were replicated four times over four blocks. Each entry was transplanted in 0.9-m² plots of three rows having 10 plants each. The entire trial was managed with recommended agronomic practices at all of the locations. Experimental layout and randomization were made using PBTools v1.4 (Sales et al., 2013). As previously mentioned, only 174 genotypes were used further for analyses.

A QTL validation population was developed by crossing Shahsarang/IR64. These two parents were selected based on their contrasting responses to Fe toxicity rather than being among the top 12 high-iron and zinc accessions identified in the study. Cross was made in 2020, and F₁ seeds were grown and tested using SSR markers to confirm hybridity. 192F₂ plants were raised along with the parental lines at ICAR-IARI, New Delhi, during the Kharif season of 2021–2022. The plants were selfed and harvested separately. The F_2:3 seeds were evaluated in the lowland area at ICAR-NEH region, Umiam, Meghalaya, during Kharif 2022–2023, and phenotyping was done for grain Fe and Zn.

At maturity, the grains were harvested from each genotype separately from five uniform-looking plants per plot from each of the four experiments. The harvested grains were pooled per genotype, cleaned by removing any discolored or ill-filled spikelets and sun-dried for 3 days to bring to uniform moisture content. The grains were then dehusked using Satake^® THU35C-T testing husker. Grain Fe (GFe) and Zn (GZn) levels were directly determined using Hitachi^® X-Supreme 8000 (Oxford Instruments Plc., UK) energy-dispersive X-ray fluorescence (ED-XRF) spectrometer calibrated as per Paltridge et al. (2012) for high-throughput screening of brown rice (Chandu et al., 2020). To this, 5 g of dehusked and cleaned rice kernels was placed in 30-mm aluminum sample cups and sealed with Poly-4 XRF film. The samples were uniformly distributed by gentle shaking before analysis. The GFe and GZn contents were recorded in mg/kg.

Soil Fe content (SFe) was estimated using the diethylenetriamine pentaacetate (DTPA) method (Lindsay and Norvell, 1978). Soil samples were collected from the experimental plots and dried sufficiently. Moreover, 20 g of air-dried sample was thoroughly mixed with 40 mL of DTPA solution, allowed to stand for 1 to 2 h, and centrifuged. The supernatant was extracted and filtered with Whatman 41 filter paper, and the filtrate was analyzed for Fe concentration using atomic absorption spectrophotometry.

For the purpose of understanding the influence of soil Fe content on grain Fe concentration, we propose an index, grain iron assimilation efficiency (GIAE), following similar indices calculated for other nutrients (Xu et al., 2012; Saito et al., 2021).

Site-wise analysis of variance (ANOVA) was used to compute coefficients of variation for phenotypic (PCV) and genotypic (GCV) components and broad sense heritability (H ²). R package augmentedRCBD (Aravind et al., 2023) was used in the computations. Using the pooled data, best linear unbiased predictors (BLUPs) were generated using a restricted maximum likelihood (REML) approach using the lme4 package integrated with the software PBTools v1.4 (IRRI 2014). Jitter box plots were drawn in R using ggplot2 package. Correlation analysis was done and presented graphically in R using corrplot package. Graphical analysis of genotype (G), site (E), and genotype-by-site interaction (GSI) was done using GGE biplot analysis using GGEbiplotGUI package and factorial regression analysis (FRA) was done using GEA-R software (https://data.cimmyt.org). FRA particularly emphasized how different co-factors of different environments influence the genotype performance (Denis, 1988; Van Eeuwijk et al., 1996; Vargas et al., 1999). A factorial regression (FR) model that incorporates environmental covariates into the genotype × environment interaction (GEI) can be expressed as:

where Y _ij is the response variable of the i ^th genotype (i=1,…,I) in the j ^th environment (j=1,…,J). µ represents the grand mean, g _i and e _j are the genotype and environment deviations respectively from the grand mean; Z _ih are the environmental covariates; ζ _jh are the genotype factor; H (H< J) is the number of environmental covariates, and ϵ _ij is the error term.

To estimate the population structure of the test germplasm, genome-wide data of 3,341,271 SNPs were downloaded from the Rice-SNP Seek database (Alexandrov et al., 2015). Filtered across taxa for a maximum of 10% missing call rate and across the sites for minor allele frequency of 0.05 and heterozygosity >0.3, a total of 140,487 markers were used for analysis. Chromosome-wise SNP distribution was visualized using an SNP density plot in SRplot software (Tang et al., 2023). The number of subgroups in the association mapping panel was estimated using the LEA-R package (Frichot et al., 2014) along with vcfR and PCA. Analysis was performed with assumed subgroups (K) ranging from 1 to 10, replicated 10 times. The total number of ancestral populations was determined using cross-entropy and the elbow method (Bholowalia and Kumar, 2014). PCA was conducted with GAPIT (Lipka et al., 2012), and the significant number of PCs was determined using a scree plot. Pairwise LD was calculated using a 50-bp sliding window in TASSEL 5.0, and LD (r²) was plotted against marker distance (Bradbury et al., 2007). Only r² values with p <0.05 were considered for LD decay analysis, and the LD decay plot was constructed in R 4.2.3 as per Remington et al. (2001).

Bayesian information and linkage disequilibrium iteratively nested keyway (BLINK) model adopted in GAPIT package was used to identify marker–trait associations (MTAs) because of its high computing efficiency and statistical power. Significant MTAs were identified after a Bonferroni multiple test correction calculated from the reciprocal of the total number of markers used for analysis [p < 4.78E-07; -log10(p) > 6.32]. Circular Manhattan plots and symphysic Q-Q plots were used for the graphical presentation of significant MTAs. The percentage of phenotypic variance explained (PVE) by individual SNP was calculated through the single-marker analysis. These MTAs were named following QTL naming conventions.

The marker interval which harbors significant marker–trait association across the locations was studied further to identify putative candidate genes for grain Fe and Zn content. The probable expressed genes present between the marker positions were downloaded from the Rice Annotation Project Database (RAP-DB) and compared as to the physical positions with those of the previously reported quantitative trait loci (QTLs). In addition to the literature survey, the Gramene QTL database (https://archive.gramene.org/qtl/) and KnetMiner (Hassani-Pak et al., 2021; https://knetminer.com) were searched to identify the physical locations of the previously reported QTLs.

The parental genotypes of the validation population were tested for parental polymorphism for the SNP-linked SSR markers using microsatellite (SSR) markers, the segregation ratio of polymorphic markers was checked using chi-square test, and distorted markers were removed (Zhang et al., 2010). Phenotypic data was tested for ANOVA, and best linear unbiased predictors (BLUPs) were generated through the REML approach using the lme4 package integrated with the software PBTools v1.4 (Sales et al., 2013) with genotypes as random variables in the model. The F₂ genotypic data was regressed on the phenotypic BLUPs of F_2:3 and the corresponding individual F₂s to perform single-marker analysis. The markers which showed significant variation among the parental classes for the target trait were considered validated for the corresponding linked QTL.

Pooled analysis of variance showed significant variations for genotype (G), site (S), and GSI components for GFe and GZn ( Table 1 ). Variation due to checks and genotype at BPN-N was not significant for GFe, while grain GZn exhibited significant variation. Boxplots across all of the sites showed the distribution for both traits, confirming the quantitative inheritance of the trait ( Figure 1A ). The average GFe and GZn content across the sites is presented in Supplementary Table S2 . GFe showed the highest mean of 22.3 mg/kg at BPN-S, with a range of 17.3 to 30.2 mg/kg, followed by DEL with an average of 17.9 mg/kg and a range of 10.5 to 33.3 mg/kg. No significant difference was observed for GFe content at ADT and BPN-N with average GFe of 14.9 and 15.9 mg/kg, respectively. Similarly, GZn content with the highest average was observed at BPN-S condition with a value of 37.2 mg/kg and a range from 34.9 to 41.9 mg/kg, followed by DEL, ADT, and BPN-N. Genetic variability analysis revealed that the highest PCV and GCV for GFe were observed to be highest at ADT (22.0 and 20.8, respectively), followed by DEL (19.6 and 18.1), BPN-S (19.2 and 17.1), and BPN-N (15.6 and 13.0). GZn displayed the highest PCV and GCV at BPN-N (19.0 and 18.3, respectively), followed by BPN- S (16.4 and 16.3), DEL (16.2 and 15.3), and ADT (14.6 and 14.0). Highest broad sense heritability was observed at ADT for GFe, while the lowest was observed for BPN-N, and the GZn maximum heritability was observed at BPN-S and the lowest at DEL ( Table 1 ).

According to the critical grain micronutrient thresholds for biofortification set by HarvestPlus—grain Fe >18 mg/kg and grain Zn >28 mg/kg—several genotypes were found to exceed these levels when grown at different sites, particularly at sites BPN-S and DEL. For grain Zn, in addition, numerous genotypes exceeded the critical level at ADT also. However, the genotypes did not consistently maintain high micronutrient levels across all sites. For grain Fe, only one genotype, IRG173, consistently showed over the threshold GFe across all sites. In the case of GZn, three genotypes—IRG189, IRG192, and IRG72—consistently exceeded the threshold across all locations ( Figure 1B ).

Pearson’s correlation coefficient between the grain micronutrient contents showed a significant association between sites ( Figure 1C ). The GFe and GZn contents between ADT and BPN-S showed a significant positive association. The GFe content at DEL showed a positive trend with that of other sites as well as with the GZn content except for BPN-N. Except for the positive association between GFe and GZn at BPN-N, GFe at this site did not show any relation with the micronutrient status at other sites, whereas the GFe recorded at BPN-S has had a significant association with GZn content at ADT, DEL, and BPN-N. Other than these, the GZn content at different sites showed a low association with the grain micronutrient contents of other sites. Exception on this were the associations observed between GZn contents BPN-S and those of ADT and DEL. Interestingly, no associations in the grain micronutrient contents could be observed between normal and stressed conditions at BPN, signifying the effect of Fe toxicity in grain micronutrient accumulation.

The average Fe assimilation efficiency (IAE) of genotypes varied significantly between the sites, and so were the variances. IAE was highest in DEL (1.20) and lowest at BPN-S (0.31). The IAE at ADT (0.62) and BPN-N (0.65) conditions were close and were intermediate to that at the other two sites. The sample variance also varied between 0.04 (BPN-S) and 0.25 (DEL) ( Figure 2B ).

Site-wise soil analysis revealed that the DEL site had a pH of 8.0, EC of 1.2 dS/m, and 16 mg/kg soil Fe content. ADT soil was alluvial, having a pH of 6.3, EC of 1.7 dS/m, and 24.0 mg/kg Fe content. The BPN-S site was highly acidic with a pH of 4.5–5.2, having a Fe content of 72.0 mg/kg and EC of 0.11 dS/m. The BPN-N site soil had a Fe content of 24.6 mg/kg with a pH range of 5 to 6 and EC 0.18 dS/m (BPN-N). The FRA demonstrated the varied levels of significance in the influences of G, S, and GSI components on both GFe and GZn. Partitioning GSI further, all three soil variables, pH, EC, and soil Fe were found to have a significant influence ( Table 2 ). Interaction of genotype × soil pH had a profound share of 47.9% in genotype × site interaction for GFe content, followed by SFe (36.1%) and EC (32.4%). A similar pattern of influence was observed on GZn also, in which genotype × soil pH interaction had 58.2% share in the GSI, followed by EC (30.9%) and SFe (21.8%).

The biplots of GGE analysis indicated the differential response of the test genotypes under different sites. In both cases, DEL was the closest to the average environment axis (AEA), and BPN-N remained the farthest. Among the genotypes, IRG169 and IRG204 showed high GFe content under BPN-S, while IRG260 showed higher GFe under DEL location. There was a distinct difference among the sites for Fe content, while there was no significant variation among the sites for GZn content. IRG193, IRG247, and IRG197 were the top rankers in GZn but were found unstable across environments. However, genotypes such as IRG142, IRG205, IRG163, IRG207, etc., were found to be stable as well as having higher GZn. In the case of GFe, IRG218 and IRG134 were found stable ( Figure 2A ).

The chromosome-wise distribution of 140,487 SNPs among the panel of 170 genotypes revealed that the markers were distributed across with an average distance of 2.66 kb ( Table 3 ). The highest marker density was found on chromosome 1 with an average marker distance of 1.99 kb and the lowest was on chromosome 9 with an average marker distance of 3.73 kb ( Figure 3A ). Linkage disequilibrium (LD) calculated based on the r ² values denoted a critical LD of 158.9 kb based on the 95th percentile of the cumulative pairwise distances between markers on each linkage group. The markers present between a physical distance of less than 159 kb tend to inherit together as a single haplotype unit ( Figure 3B ). The cross entropy from the estimated admixture coefficients plotted against the expected number of subpopulations revealed a significant fall in the slope at K = 3, revealing the elbow point ( Figure 3C ). The three sub-populations were designated as POP1, POP2, and POP3. A PCA analysis performed on the genotypic data also indicated a significant presence of three groups in the population structure, explaining about 40% of the total variation ( Figure 3D ). With a maximum cutoff value of 0.95 for admixture (co-ancestry) coefficients (Q), 49.4% of admixtures were identified in the whole population. When plotted in the bar chart, POP1 had 30 genotypes, of which 18 genotypes were with less than 5% admixing, and the remaining 12 had 5% to 41% admixing. POP2 contained 14 genotypes, nine of which were least admixed and the remaining five were admixed. Among 126 genotypes in POP3, 59 were least admixed and 67 were admixtures ( Figure 3E ). The details of the genotypes in the three subpopulations are provided in Supplementary Table S3 .

The Q-Q plots for all the traits showed the expected distribution of log₁₀(p) that implied perfect control of false discovery ( Figure 4 ). Among the six MTAs identified, five were mapped for GFe and one for GZn content. Three MTAs detected for GFe at ADT were distributed on chromosomes 1(qGFe1.1 ^ADT), 8 (qGFe8.1 ^ADT), and 12 (qGFe12.1 ^ADT) and positioned respectively at 2.6, 19.4, and 18.5 Mbp. These MTAs explained a phenotypic variance of 10.8%, 8.3%, and 9.7%, respectively. At BPN -S, the only MTA detected was on chromosome 2 (qGFe2.1 ^BPN-S) located at 0.48 Mbp, which explained a PVE of 42.3% for GFe. There were two MTAs identified under BPN-N both located on chromosome 12. One of these MTAs was associated with GFe (qFe12.2 ^BPN-N) and was located at the physical position of 21.2 Mbp, explaining 41.7% variation. The other MTA was associated with GZn (qGZn12.1 ^BPN-N) and was located at 14.2 Mbp, accounting for 32.2% of the total phenotypic variance ( Table 4 ).

A total of 92 annotated gene models were identified within the haplotype surrounding the identified MTAs. These included 42 candidate genes belonging to different protein family classes such as auxin transporters, cation efflux proteins, heavy-metal-associated domains, nicotinamine synthase, 6-phosphoribosyl transferase domains, iron-regulated metal transporters, and zinc ion binding proteins. qGFe1.1 ^ADT was found to be associated with genes like OsMT2D (Os01g0149200) and OSMT2A (Os01g0149800) encoding metallothionein (MT) and controlling iron and zinc homeostasis. qGFe2.1 ^BPN-S was found to be associated with Cytochrome p450, CYP1, and HPL1 gene encoding heme-binding protein, while qGFe8.1 ^ADT was found proximal to genes OsHMP39, HMP39, OsHIPP46, OsaHIP46, HIP46, OsHMP40, HMP40, OsHPP7, OsaHPP07, and HPP07 (Os08g0403300) encoding heavy-metal-associated protein 39 and heavy-metal-associated isoprenylated plant protein 46. Candidate genes around qGFe12.1 ^ADT included genes encoding OSM protein (oxidative stress management genes), while candidate genes around qGFe12.2 ^BPN-N were in close proximity with the genes OSGRX28 and OSGRX29 and genes encoding metal binding protein. The MTA was also found to be associated with genes such as VIL1 (Os12g0533500) and VIP5 (Os12g0535900) associated with grain yield and biomass. The genes associated with qGZn12.1 ^BPN-N include OsFUSED, CCD1, and LAP1, controlling traits like pollen development, lipid transport, and zinc ion binding protein ( Supplementary Table S4 ).

The ANOVA of the F_2:3 population derived from Shahsarang/IR64 showed significant differences for GFe and GZn contents among the segregating progenies. The frequency distribution histogram of F_2:3 demonstrated that both of the traits were normally distributed in the population, confirming the quantitative inheritance pattern ( Figure 5A ). The phenotypic evaluation revealed GFe ranging from 8.1 to 31.1 mg/kg, with a mean value of 15.2 mg/kg and high PCV and GCV, while GZn ranged from 22.1 to 39.0 mg/kg, with a mean value of 28.0 mg/kg and low GCV and medium PCV. Both of the traits exhibited high heritability ( Supplementary Table S6 ).

Using SSR markers, two QTLs out of six, qFe12.2 ^BPN-N and qZn12.1 ^BPN-N, were validated. There were 30 SSR markers present in the vicinity of the MTA region on chromosome 12. A polymorphism survey between the two parents identified six polymorphic markers between them ( Supplementary Table S7 ). The marker, RM1986, showed a significant association with the GFe under both conditions (BPN-N and BPN-S) and RM179 with GZn at BPN-S. The genotypes which carry the homozygous allele of IR 64 (RM1986) showed an average GFe of 12.5 mg/kg at BPN-S and 10.6 mg/kg at BPN-N, whereas genotypes which carry homozygous Shahsarang alleles showed 17.6 mg/kg at BPN-S and 12.6 mg/kg at BPN-N. The PVE of the 12:21260572 SNP was found to be 55.3% and 16.4% with additive effects of 0.9 and 2.5 mg/kg at BPN-S and BPN-N, respectively, indicating that the positive QTL allele was contributed by Shahsarang. Genotypes with homozygous RM179 alleles of IR 64 showed an average GZn content of 27 mg/kg, whereas genotypes with homozygous Shahsarang alleles showed 31.5 mg/kg. The PVE of the marker was 42.7%, and the additive effect was found to be 2.2 mg/kg, indicating that the QTL was contributed by Shahsarang ( Table 5 ).

The qFe12.2 ^BPN-N was found networked with two thioredoxin fold domain containing protein genes, Glutaredoxin 28 (OsGRX28) and Glutaredoxin 29 (OsGRX29), located within the 21.4–21.5-Mbp region. Another prominent network associated with qZn12.1 ^BPN-N was linked to the OsFUSED gene (Os12g0433500) ( Figure 5B ).