Meta‐QTL s and haplotypes for efficient zinc biofortification of rice

Biofortification of rice with improved grain zinc (Zn) content is the most sustainable and cost‐effective approach to address Zn malnutrition in Asia. Genomics‐assisted breeding using precise and consistent Zn quantitative trait loci (QTLs), genes, and haplotypes can fast‐track the development of Zn biofortified rice varieties. We conducted the meta‐analysis of 155 Zn QTLs reported from 26 different studies. Results revealed 57 meta‐QTLs with a significant reduction of 63.2% and 80% in the number and confidence interval of the Zn QTLs, respectively. Meta‐quantitative trait loci (MQTLs) regions were found to be enriched with diverse metal homeostasis genes; at least 11 MQTLs were colocated with 20 known major genes involved in the production of root exudates, metal uptake, transport, partitioning, and loading into grains in rice. These genes were differentially expressed in vegetative and reproductive tissues, and a complex web of interactions were observed among them. We identified superior haplotypes and their combinations for nine candidate genes (CGs), and the frequency and allelic effects of superior haplotypes varied in different subgroups. The precise MQTLs with high phenotypic variance, CGs, and superior haplotypes identified in our study are useful for an efficient Zn biofortification of rice and to ensure Zn as an essential component of all the future rice varieties through mainstreaming of Zn breeding.


INTRODUCTION
Zinc (Zn) is one of the most important mineral elements required for the successful completion of life cycle in both plants and animals (Yashona et al., 2018).It provides structural and functional stabilities to hundreds of enzymes, and proteins thereby play a vital role in cell division, immune function, cognitive development, sensory function, reproductive health, and neurobehavioral development in human beings (Cakmak & Hoffland, 2012;Gibson, 2012;Stein, 2010).Therefore, an adequate daily intake of Zn is recommended for the proper growth, development, and to maintain health.The prolonged Zn deficiency leads to stunting, diarrhea, vulnerability to diseases, poor cognitive development, infant and maternal mortalities, and so on (Brooks et al., 2004;Calayugan et al., 2021;Sazawal et al., 2007;Tielsch et al., 2007;Young et al., 2014).These health issues are widespread globally but highly visible in the populations that are entirely dependent on staples for their nutritional needs (Hotz & Brown, 2004;Stein, 2010;Cakmak & Hoffland, 2012;Gibson, 2012).Zn malnutrition affects over two billion people worldwide; especially children, women, and elderly are highly vulnerable.Thus, addressing malnutrition has been accorded a major global health priority (UN General Assembly, 2015).Biofortification of food crops with enhanced nutritional value is the most cost-effective and a sustainable food-based solution to combat malnutrition (Bouis et al., 2011;Harvest-Plus, 2021;Siwela et al., 2020;White & Broadley, 2009).As a result, several biofortified crop varieties with higher Zn, Fe, and vitamin A contents have been commercialized, and subsequent bio-efficacy studies using biofortified crops were found to be effective in improving the nutritional status (De Moura et al., 2014;Sazawal et al., 2018).In rice, several Zn biofortified varieties with an additional 8-10 ppm of grain Zn content have been released for commercial cultivation in the target countries of Asia, and efforts are being made to deploy them widely (Calayugan et al., 2021;Inabangan-Asilo et al., 2019).However, there is an urgent need to expedite the development of Zn biofortified rice and also to mainstream the grain Zn as an essential component of the rice breeding product profiles, so that all the future rice varieties will have enhanced levels of grain Zn content.A better knowledge of the molecular basis of grain Zn accumulation, factors that influence the stable performance of high Zn lines, and combining a suite of agronomic and other desirable traits along with high grain Zn content are essential for the success of rice Zn biofortification program.Efforts are being made to mainstream Zn in rice breeding by using advanced breeding techniques, such as rapid generation advance, speed breeding, recurrent selection, and multi-trait genomic selection (Calayugan et al., 2021;CGIAR, 2018).

Core Ideas
• Meta-analyses of 155 Zn QTLs identified 57 MQTLs with reduced confidence intervals.• Major metal homeostasis genes colocated with MQTLs were prioritized.• In silico expression and co-expression analyses revealed the complex network of genes involved in metal homeostasis.• Identified superior haplotypes are useful for efficient rice Zn biofortification.
The enormous wealth of genetic resources and molecular information available for rice has enabled the genomic dissection of complex traits, including grain Zn and Fe contents, to identify major quantitative trait loci (QTLs) or genes (Jiang et al., 2012;Swamy et al., 2016;Yang et al., 2013).Although several grain Zn QTLs have already been identified, very few of them have been successfully used in the breeding programs, mainly because of their small effect size, unavailability of diagnostic markers, large confidence intervals (CIs), and many of them have not been either validated or fine mapped (Gong et al., 2018;Mackill, 2006).In addition, QTL identification is highly influenced by the choice of parents, markers, size and type of population, and environmental factors leading to less accuracy or precision of QTLs location or effects, which obstructs their transferability across different breeding programs (Bradshaw et al., 1995(Bradshaw et al., , 1998;;Garrin et al., 2021;Li et al., 2003;Yu et al., 2011).Thus, an identification of genetically and environmentally stable major effect Zn QTLs and their use in genomics-assisted breeding (GAB) is necessary to fast-track the development of Zn biofortified rice varieties.
Meta-QTL (MQTL) analysis consolidates QTLs mapped for a trait in multiple independent genetic backgrounds and environments to provide precise, stable, and narrow genomic regions for a cluster of QTLs (Arcade et al., 2004;Goffinet & Gerber, 2000).Recently, many studies have reported MQTLs for grain yield, biotic and abiotic stress tolerances, grain quality, and nutrition-related traits in different crops (Gudi et al., 2022;Guo et al., 2018;Saini et al., 2022;Sheoran et al., 2022;Soriano et al., 2021).MQTLs with improved genetic resolution have been identified for grain Zn and Fe contents in rice, wheat, maize, and common beans (Izquierdo et al., 2018;Jin et al., 2013Jin et al., , 2015;;Raza et al., 2019;Shariatipour et al., 2021;Singh et al., 2022).All these studies have clearly shown that MQTL analysis is a powerful tool to identify stable QTLs with narrow CI.Particularly in rice, there are already 2 studies reporting 22 and 46 MQTLs for grain Zn content, and 7 of these MQTLs were common; there was a 75% reduction in the total number of QTLs, and also, their CIs were significantly reduced (Jin et al., 2015;Raza et al., 2019).But recently, single-nucleotide polymorphisms (SNPs) have been widely used in QTL mapping and genome-wide association studies (GWAS) for grain Fe and Zn contents in rice.Thus, there is a need for an updated meta-analysis by consolidating Zn QTLs identified using high-density maps.
Major effect MQTLs are more useful for GAB and to mine candidate genes (CGs) for designing functional markers and for their functional characterization (Collard & Mackill, 2008;Islam et al., 2019).With the availability of multi-omics platforms, now it is easier to understand the functions of CGs, to assess their diversity, and identify superior haplotypes in the germplasm panels such as 3K-RGP (3K rice genome panel) (Abbai et al., 2019;Bevan et al., 2017;Li et al., 2014;Sato et al., 2011;Varshney et al., 2018).
Therefore, the major objectives of our study were to conduct the review of literature to compile QTLs identified for grain Zn content, to identify consensus and precise Zn-related QTL regions through MQTL analysis, to mine CGs associated with Zn metabolism and understand their expression pattern, co-expression network, and to identify superior haplotypes of selected CGs.The results from this study will be useful for the efficient Zn biofortification of rice through GAB.

Literature search to compile Zn QTL information
We carried out Google Scholar search to compile QTLs identified for grain Zn in rice; preferably those identified using simple sequence repeats (SSRs) and SNP markers.We retrieved Zn QTL information from 2007 to 2020 using various keywords, such as QTL mapping, association mapping, Zn, minerals, and rice.The unpublished Zn QTL results from MS and Ph.D. thesis conducted at the International Rice Research Institute (IRRI) were also used for the analysis.In all 26 QTL mapping studies were utilized for the analyses.The details of the parental lines used in mapping population development, population size, number and type of markers used, range of Zn content in the population, and number of QTLs identified are presented in Table 1.The QTL information compiled are QTL names, linkage group, peak position, the logarithm of odds (LOD), phenotypic variance explained (R 2 ), and CI.Wherever peak QTL position information was lacking, the mid-point between the two flanking markers was used.Similarly, a LOD threshold of 3 was used for the missing LODs.The following formulas were used to estimate the missing CIs of QTLs: where N represents the population size, and R 2 is the phenotypic variance percentage explained by individual QTL.Equation (I) was used to estimate the CI of recombinant inbred lines, whereas Equation (II) was used for backcrossderived populations and F 2 populations (Darvasi & Soller, 1997;Weller & Soller, 2004).Equation (III) was used to calculate missing R 2 (Nagelkerke, 1991).The CI of the QTLs from the GWAS was considered 0.5 Mb each upstream and downstream from the position of maker trait associations.
After gathering all the information, a separate QTL input file was prepared for each chromosome.

Consensus map construction and QTL projection
We constructed a consensus map using 2180 SSRs, 1K-RiCA SNPs, and 7K Infinium SNPs (Arbelaez et al., 2019;IRGSP, 2005;Morales et al., 2020).The missing flanking markers information of the QTLs was retrieved from Gramene and added to the consensus map using their physical positions as a reference (http://www.gramene.org).All the markers used in this study were arranged on their respective chromosomes based on their physical positions.Then the physical positions were converted to cM using the conversion factor 1 cM = 250 kb for building the consensus map (Raghavan et al., 2017).A separate map file was prepared for each chromosome.The QTL and map input files were used for the analysis.The QTL file was projected onto the consensus map to consolidate all the QTL information.

2.3
Meta-QTL analysis MQTL analysis was performed on the QTL cluster using BioMercator v4.2 (Arcade et al., 2004).We followed the Veyrieras approach for chromosomes with >10 QTLs, and the Goffinet and Garber method for chromosomes with ≤10 QTLs (Veyrieras et al., 2007).The Akaike information  criterion (AIC) was used to select the best model for defining the number of MQTLs.Model with the lowest AIC value is regarded as the best model because it corresponds to the least amount of information lost (Akaike, 1978;Sosnowski et al., 2012).The average R 2 of the initial QTLs was used to calculate the R 2 of the MQTLs.

Identification of CGs underlying MQTL regions
All the MQTLs were subjected to CG analysis.The physical positions of the markers that flank each MQTL were used to extract the underlying genes within the MQTLs.The batch retrieval option was used to extract all the genes within each MQTL region from the Rice Annotation Project Database (RAP-DB) (Sakai et al., 2013).Thereafter, genes were filtered based on gene ontology (GO) term or description related to Zn ion binding/transport, Fe ion binding/transport, heme binding, cation ion binding/transport, and metal ion binding, transport, and homeostasis.GenFam was used to determine the enriched gene families among the filtered genes using the options of the Fisher exact test for statistical enrichment, and the Benjamini-Hochberg method to control false discovery rate (Bedre & Madnani, 2019).The higher the -log10(p-value) value, the greater the confidence in the enrichment of the gene family.Finally, 20 CGs involved in Zn or Fe homeostasis were selected for the in silico expression, co-expression, and haplotype analysis.

Spatiotemporal and co-expression analyses
RiceXPro v3.0 (https://ricexpro.dna.affrc.go.jp/) was used for the spatiotemporal expression analysis of the selected CGs using the RXP_0001 dataset (Sato et al., 2011).Then the Rice Functionally Related gene Expression Network Database (RiceFREND v2.0) (https://ricefrend.dna.affrc.go.jp/) was used for the construction of co-expressed gene network for each gene using single guide search (Sato et al., 2013).For the co-expression analysis, the mutual rank (MR) was set to 10 and the hierarchy to 1.A lower MR value represents stronger co-expression.

Haplotype analysis of selected CGs across the 3K-RGP
The selected CGs were subjected to haplotype analysis using the inbuilt tool of the SNP seek database (Mansueto et al., 2017).The "3kfiltered" SNP set representing all the 12 subpopulations was utilized for the analysis (3K-RGP).Only nonsynonymous or functional SNPs were considered for the analysis.Heterozygotes were considered missing, and then, genotypes with more than 20% missing SNPs were removed before carrying out the analysis.Flapjack 1.22.04.21 was used to visualize the haplotypes of selected genes (Milne et al., 2010).

Identification of superior haplotypes
We used previously generated Zn phenotypic data of aus panel at IRRI.These trials were conducted using an alpha lattice design with two replications in 2014WS (Wet Season) and three replications each in 2015WS and 2015DS (Dry Season), with a varying number of blocks in each season.The grain Zn content was determined in polished rice using X-ray fluorescence (XRF)-Bruker S2 RANGER.The best linear unbiased estimates (BLUEs) of Zn data for each line were estimated using PBTools v1.4.A total of 134 of these accessions were found to be common between the seasons and were part of the 3K-RGP.The BLUEs of grain Zn content were used for the phenotypic validation of haplotypes to identify the superior haplotypes (Table S1).The haplotypes found in at least two lines were considered for the analysis.Further, the relationship between Zn content and haplotypes was investigated to identify superior haplotypes.There were lines with more than one superior haplotypes, so haplotype combination effects on grain Zn content were also estimated.The F-statistic and Tukey honestly signficant difference (HSD) tests were used to investigate significant differences among haplotypes and haplotype combinations at p < 0.05 (Tukey, 1949).

Literature search for grain Zn QTLs
Bibliographic search for the published papers and thesis on Zn QTL mapping identified 26 studies.These included both association mapping and bi-parental QTL mapping, but majority of the studies were carried out on bi-parental populations, namely, RILs, BILs DH, and F 2S , which were developed using elite rice varieties/breeding lines crossed with high Zn donor lines, such as landraces, wild rice, aus accessions, and flood-tolerant rice varieties.The association mapping studies were conducted on multiparent advanced generation intercross (MAGIC) populations and germplasm panels.The Zn content was measured either in brown rice or milled rice or both by using advanced instruments, such as XRF, ICP, and atomic absorption.Predominantly SSRs and SNPs were used for mapping, the number of markers varied from 45 to 31k.Furthermore, the population size ranged from 50 to 378 lines or accessions.The Zn content among the populations ranged from 0.4 to 157.4 ppm.In all 220 Zn QTLs reported, the number of QTLs per study varied from 2 to 17, with 5 or more QTLs identified in 21 studies, and their R 2 values varied from 0.01% to 47.6% (Table 1).The QTLs lacking physical position information, very low R 2 values (<1%), and those without overlapping CIs with other QTLs at the same chromosomal region were excluded from the analyses.After careful assessment, a total of 155 Zn QTLs were finalized for further analyses, which were distributed throughout the genome (Figure 1).

Consensus map and MQTL analysis
We constructed a high-density consensus map using 10,211 markers (2230 SSRs and 7981 SNPs) that included all the flanking makers of 155 Zn QTLs.The consensus map had a length of 372.24 Mb or 1488.9 cM with a marker density of one marker per every 36 kb.As expected, chromosome 1 with 1294 markers had the highest length of 172.9 cM, whereas chromosome 9 with 625 markers had the shortest length of 90.72 cM.On other chromosomes, the number of markers varied from 633 to 1027.Similarly, the map length on other chromosomes also varied from 92.3 to 145.4 cM.Because we used the physical position of the markers, their grouping and order on the chromosomes was accurate, which helped us to precisely project 155 Zn QTLs from different studies on to the consensus map.The number of projected QTLs varied from 5 on chromosome 11 to 20 QTLs on chromosome 5, and in 10 of the 12 chromosomes ≥10 QTLs were distributed (Figure 1).We performed meta-analyses on QTL clusters chromosome-wise and used least AIC to declare the number of MQTLs at each of the QTL clusters.A total of 57 MQTLs were identified across all the 12 chromosomes.The MQTLs identified on different chromosomes ranged from 2 to 7. The highest number of seven MQTLs was detected on chromosomes 3-5, whereas lowest two MQTLs were detected on chromosome 1.MQTL 11.1 and MQTL 11.2 were found to be overlapping.The number of initial QTLs involved in MQTLs ranged from 2 to 8, but in majority of the MQTLs (49), the initial QTLs involved were 2-4.Four MQTLs (MQTL 1.2 , MQTL 5.4 , MQTL 6.3 , and MQTL 7.2 ) were detected from a cluster of five QTLs each on chromosomes 1, 5, 6, and 7, whereas three MQTLs (MQTL 5.1 , MQTL 7.5 , and MQTL 7.6 ) involved six initial QTLs on chromosomes 5 and 7, and highest eight QTLs cluster resulted in only one MQTL (MQTL 12.2 ) on chromosome 12.The CIs ranged from 0.05 Mb (MQTL 7.4 and MQTL 10.3 ) to 3.7 Mb (MQTL 3.3 ).In 26 MQTLs (46%), the CI was ≤300 kb, whereas in 8 MQTLs, it was 300-500 kb, and in the rest of the MQTLs, CI was >500 kb.There was an average 80% reduction in the CI of MQTLs compared to their mean CI of initial QTLs.The R 2 of MQTLs varied from 2.5% to 45.5%, and 30 (52.6%) MQTLs had R 2 > 10%.The identified MQTLs are shown in Figure 2, and details of each MQTL are presented in Table 2.

Mining of CGs underlying MQTL regions
A total of 4196 genes within MQTL regions were obtained from RAP-DB.Following filtering based on the functional role and GOs, 318 genes were selected within 51 MQTLs (Table S2).Based on GO, 35.8% of genes were associated   with Zn ion binding and Zn ion transport, 14.2% with Fe ion binding and transport, 9.7% with cation binding and transport, 9.7% with metal ion binding and transport, 6.3% with heme binding, and 18.2% genes were enriched with oxidation-reduction process, cell redox homeostasis, electron carrier activity, and transmembrane transport.A new computational approach, GenFam allows annotation, classification, and enrichment of genes based on their gene family.Among 318 genes identified, 190 genes were annotated and classified into gene families by GenFam, resulting in 60% interaction/coverage with the GenFam database.A total of 49 gene families were identified with 23 of them being overrepresented (p-value <0.05).The most enriched gene family was the RING finger domain gene family followed by Cytochrome P450.The other enriched families include the Zn finger gene family, copper chaperone gene family, glycoside hydrolase gene family, and 2OG-Fe(II) oxygenase gene family (Table S3 and Figure S1).

Expression profile and co-expression network of selected CGs
For all the genes except OsMTP7 and OsbHLH133, expression information was available in the expression database (RXP_0001).The heat map shows the time of sampling, the amount of tissue sampled for expression analysis, and the level of expression of different genes in 12 different vegetative and reproductive tissues at different stages of their development (Figure S2).Most of the genes except OsMTP1 and OsNRAMP7 were differentially expressed with moderate-tohigh level of transcripts in both vegetative and reproductive tissues.However, they were predominantly expressed in leaf blades, leaf sheaths, and roots compared to reproductive tissues.Some of the genes were highly expressed in stem (OsFRDL1, OsZIP5, and OsMT1a), ovary (OsFRDL1, OsD-MAS1, OsNAS3, and OsTOM2), inflorescence (OsFRO1), T A B L E 3 List of selected candidate genes and their chromosomal distribution.anther (OsRab6B2, OsDMAS1, OsFRO1, OsMT1a, and OsMT1g), and pistil (OsZIP5).Similarly, OsZIP5 is expressed in roots, leaf blade, leaf sheath, stem, inflorescence, pistil, lemma, palea, and ovary.However, OsTOM1 was expressed in roots and leaf sheath, OsTOM2 expressed in ovary, and in almost all the stages of embryo and endosperm development, including roots and leaf sheath.High expression was also observed in endosperm for OsRab6B2, OsDMAS1, OsMT1f, and OsMT1g.

RAP ID
Results of the co-expression network analyses are presented in Figure S3.OsNAS3 was found to co-express with OsNAAT2 and other genes involved in oxidoreductase activity and abiotic and biotic stress tolerances.OsDMAS1 co-expressed with OsMTK2, a gene involved in methionine biosynthesis.OsSAMS1 co-expressed with LOC_Os12g42876, a gene having Zn binding function, as well as other genes that control shoot branching, nitrogen use efficiency, and the synthesis of secondary metabolites.The OsNRAMP1, OsIRO2, OsIMA1, and OsIMA2 involved in Fe homeostasis were co-expressed with OsFRO1.OsZIP5 displayed co-expression with trehalose biosynthetic gene.The genes associated with signal transduction and gluconeogenesis co-expressed with OsZIP7.Zn transporters OsZIFL9 and OsTOM1 co-expressed.Similarly, OsTOM2 and OsTOM3 co-expressed with Zn transporters OsZIFL10 and OsZIFL12, respectively, as well as with other genes involved in cellular metabolic processes and oxidoreductase activity.The co-expressed genes of OsMT1f were associated with copper stress, salt stress, and pathogen infection.In the case of OsMT1d, the co-expressed genes like OsSa1T, OsHAC4, OsEXPB2, and OsMT1Ld were related to salt tolerance, arsenic tolerance, root system architecture, and metallothionein, respectively, whereas OsMT1a also showed co-expression with OsSa1T.Regarding the co-expression network of OsRab6B2, genes regulating biotic stress tolerance, carbohydrate metabolism, and transmembrane transport were present.

Identification of superior haplotype of selected CGs
Haplotypes were identified for the selected 20 CGs using the rice SNP seek database representing 3K-RGP.The OsMTP1 did not reveal any haplotype.For the remaining genes, the number of haplotypes varied from 2 to 16 (Figure 3a,b and Figure S4).Two haplotypes for 5 genes, 3 haplotypes for 6 genes, 4 haplotypes for 4 genes, 5 haplotypes for 2 genes, and 11 and 16 haplotypes for 1 gene each were identified from 3K-RGP.The frequency of haplotypes of the selected genes is presented in Figure 3e.
We carried out the validation of haplotypes and haplotype combinations using the Zn data of 134 aus accessions evaluated for three seasons at IRRI.The range of Zn content of various haplotypes across three seasons is presented in Table S4.However, five genes had only one haplotype and were removed before conducting the analysis (Figure 4a).Finally, 14 genes were selected for the analysis, out of which 5 genes did not show any significant haplotype effects (Figure 4b).The haplotypes of nine genes were further tested by mean comparison using the Tukey HSD post hoc test, and superior haplotypes were identified (Table 4).The superior haplotypes with highest Zn content were Hap1 for OsZIP5 (17.66 mg/kg), OsNRAMP7 (21.17 mg/kg), and OsMT1f (17.63 mg/kg); Hap2 for OsRab6B2 (18.13 mg/kg), OsDMAS1 (19.77 mg/kg), and OsZIP7 (18.35 mg/kg); Hap3 for OsZIP9 (17.77 mg/kg) and OsFRO1 (18.35 mg/kg), and Hap4 for OsMT1g (21.92 mg/kg), respectively (Figure 4c,d and Figure S5).The Zn content of the lines increased with an increase in the number of superior haplotypes, among the six superior haplotype combinations for nine genes, lines with five haplotypes (20.6 mg/kg) were significantly different from all the other groups, whereas there was a nonsignificant decrease in the Zn content of lines with six superior haplotypes over lines with five superior haplotypes (Table S5).
The distribution of superior haplotypes of the abovediscussed nine genes varied in different subpopulations of 3K-RGP, and some of them were predominantly present in certain groups (Figure 3c,d and Figure S6).The superior haplotypes for OsZIP7, OsFRO1, and OsRab6B2 were found in high frequency in aus subgroup, OsZIP9 and OsDMAS1 in indica, whereas superior haplotypes for OsZIP5 and OsMT1f were found in aus and indica, whereas for OsMT1g in japonica, subtropical, temperate, tropical, and aromatic subgroups, whereas superior haplotype for OsNRAMP7 was distributed across the subgroups with low frequency in aus.

DISCUSSION
Biofortification of rice with high grain Zn content is one of the feasible options to address Zn malnutrition in rice-consuming Asian countries (Bouis & Saltzman, 2017;Kassebaum et al., 2014;Nestel et al., 2006;Wessells et al., 2012).Several conventional bred high Zn rice varieties have been released for commercial cultivation but the efforts are targeted, isolated, and the whole process has been slow (Calayugan et al., 2021;Inabangan-Asilo et al., 2019).Mainstreaming of grain Zn trait as an important component of all the future rice varieties helps to address Zn malnutrition on a larger scale (Palanog et al., 2019).In our study, we carried out MQTL analyses of Zn QTLs, shortlisted candidate metal homeostasis genes, and identified superior haplotypes to support mainstreaming the Zn biofortification of rice.
It is important to note that >200 QTLs for grain Zn and Fe have been reported in rice (Raza et al., 2019).If this valuable public genomics resource is processed and used in Zn breeding, it has the potential to enhance the efficiency of rice Zn biofortification (Swamy et al., 2016).MQTL analyses is a popular, user-friendly, and in silico approach for the consolidation of all the QTLs for a trait, and identification of MQTLs with stable expression and narrower CIs enables their use in MAB and GS (Arcade et al., 2004;Goffinet & Gerber, 2000;Gong et al., 2018;Norton et al., 2008;Swamy et al., 2016).As a result, a number of MQTL studies were carried out for various traits in multiple crops reporting significant reduction in redundant QTLs and their CIs (Guo et al., 2018;Islam et al., 2019;Khahani et al., 2020;Soriano et al., 2021;Swamy et al., 2011;Selamat & Nadarajah, 2021;Saini et al., 2022;Sheoran et al., 2022).But it is worth noting that most of these studies used sparse linkage maps with different types of markers, resulting in a lack of common markers among studies and with that of reference map leading to difficulties in QTL consolidation, to accurately locate MQTLs and to precisely narrow down their CIs (Cai et al., 2019;Jin et al., 2015;Raza et al., 2019;Zhang et al., 2017).In two previous MQTL studies on grain Zn content in rice as well sparse linkage maps were used (Jin et al., 2015;Raza et al., 2019).Therefore, the use of highdensity linkage maps is useful in identifying stable MQTLs with much narrower CIs for use in breeding and in prioritizing CGs for functional characterization.
We used genetic studies conducted on grain Zn content in rice utilizing mainly SNPs and SSRs.We also constructed a high-density consensus map that consisted of 10,211 markers (SNPs and SSRs) with an average marker interval of 0.14 cM.The meta-analyses of 155 Zn QTLs summarized them to 57 MQTLs with a significant reduction of 63.2% in Zn QTLs, which suggests that most of the previously mapped Zn loci were shared by multiple populations.The physical intervals of 15 and 19 MQTLs identified in the current study matched with the MQTLs reported in the prior 2 studies, respectively, and 5 of these MQTLs were shared by all 3 studies (Jin et al., 2015;Raza et al., 2019) (Table S6).The CIs were also reduced on average by 80%, which is higher than achieved in previous studies, thus improving the mapping resolution of the genomic regions.The CIs of MQTLs ranged from 0.05 to 3.66 Mb, whereas they ranged from 0.07 to 11.64 Mb in previous reports (Jin et al., 2015;Raza et al., 2019).In our study, 46% of the MQTLs had CI < 300 kb, 14% had 300-500 kb, 29% had 500 kb-1 Mb, and 10.5% had >1 Mb.In previous study by Raza et al. (2019), only 8.6% MQTLs had <300 kb, 8.7% with 300-500 kb, 21.7% with 500 kb-1 Mb, and 61% had >1 Mb, whereas all the MQTLs except one locus had CI > 1 Mb in another study (Jin et al., 2015).All the common MQTLs except MQTL3.3 had narrower CIs in our study.The R 2 of MQTLs varied from 2.5% to 45.5% with 52.6% of the MQTLs having >10% R 2 .Further, five MQTLs (MQTL 1.2 , MQTL 5.4 , MQTL 6.3 , MQTL 7.5 , and MQTL 12.2 ) were found to be some of the interesting MQTLs based on the following criteria: average phenotypic variance >10% and involvement of a large number of initial QTLs within the MQTL.These MQTLs have the potential to be the most effective candidate targets for use in marker-assisted selection (MAS) after being validated in various genetic backgrounds and environments utilizing linked flanking markers.The MQTLs identified in this study provide more closely linked markers for use in GAB to improve grain Zn content in rice.
However, a deeper insight into the MQTLs revealed that at least 20 major metal homeostasis genes were associated with 11 MQTLs (Table 3).Most of the genes were differentially expressed in both vegetative and reproductive tissues.Overexpression of OsNAS3 in roots and leaf blades in rice has been previously reported (Inoue et al., 2003).Likewise, the expression of OsTOM2 in endosperm and OsTOM1 in roots is comparable to that previously described (Ricachenevsky et al., 2011).The majority of these genes co-expressed with other metal homeostasis genes, genes involved in growth and development, and biotic and abiotic stress tolerances.OsNAAT2 was found to be strongly coexpressed with OsNAS3 suggesting its role in Fe/Zn uptake.In similar manner, the co-expressions of OsFRO1 with OsNRAMP1, OsIRO2, OsIMA1, and OsIMA2 suggest their coordinated involvement in Fe/Zn homeostasis.OsMTK2, a 5-methylthioribose (MTR) kinase gene, has role in Met (Methionine) cycle (Sauter et al., 2004).According to the study by Aung et al. (2018), OsDMAS1 and OsMTK2 were both found to be suppressed in roots under Fe toxicity, and it was also suggested that the Met cycle might play a role in nicotianamine synthesis under excess Fe.Similarly, OsTOM1 was found to be co-expressed with OsZIFL9, OsTOM2 with OsZIFL10, and OsTOM3 with OsZIFL12, probably due to their high similarity (Ricachenevsky et al., 2011).However, more confirmatory research is needed to validate and comprehend the co-expression network, which could lead to new insights into plant responses to nutrient deficiency.There have been several earlier studies on physiological, molecular, and transgenic validations of these genes proving their role in the production of root exudates, chelation, and minerals uptake, transport, partitioning, and loading into grains in rice (Table 3).It clearly shows that micronutrient movement from soil to grain is a complex process, which involves a coordinated function of the multiple metal homeostasis genes at different stages of crop growth, especially during reproductive stage to determine the grain mineral content of rice varieties (Sperotto et al., 2014;Swamy et al., 2016;Waters et al., 2012;Xia et al., 2020Xia et al., , 2021)).
Identification of haplotypes for different genes and their use in crop improvement through allele mining, pyramiding, or genomic selection is found to be efficient in achieving better genetic gains (Abbai et al., 2019;Addison et al., 2020;Anandan et al., 2022;Bhat et al., 2021;Du et al., 2022).Using genome editing approach, two yield-related genes Gn1a and DEP1 were altered and genome-edited lines revealed that one mutant allele of Gn1a and three of DEP1 provides higher yield as compared to other natural high-yield alleles (Huang et al., 2018).We identified 2-16 haplotypes for 19 metal homeostasis genes in 3K-RGP, and their distribution varied across the subgroups.We were able to phenotypically validate haplotypes of 14 genes in a subset of aus lines and identified superior haplotypes for 9 genes (Figure 4b, Table 4).The frequencies of superior haplotypes varied in different subpopulations but were predominantly present in certain groups (Figure S6).The superior haplotypes for OsZIP7, OsFRO1, OsRab6B2, OsZIP5, and OsMT1f were in high frequency in aus and indica subgroups, OsZIP9 and OsDMAS1 in indica, whereas OsMT1g in japonica and aromatic subgroups.Previous research has also found that haplotype distribution is determined by evolutionary and population genetic factors such as mutation and recombination rates, as well as selection (Magwa et al., 2016;Zaitlen et al., 2005).It provides an idea to the breeders for the selection of best donors for their trait of interest from the gene bank.For example, aus accessions were found to be a rich source of micronutrients, and several high Zn donor accessions, such as Kaliboro (IRGC 77201-1), Jamir (IRGC117765), Lalsaita (IRGC 43915-1), and UCP122 (IRGC 8794-1), have been identified and extensively used in the Zn biofortification program at IRRI (Calayugan et al., 2021;Palanog et al., 2019;Swamy et al., 2016).
It is worth noting that in this study, the Zn content of the lines increased with an increase in number of superior haplotypes; the group of lines with five superior haplotypes was significantly different from all the other groups and had a Zn increment of 2.5 mg/kg (Table S5).Previous studies showed that the pyramiding of QTL superior alleles improved the seed storability (Wu et al., 2021).In the same way, haplotype analysis of 21 yield-and grain-quality-related genes in rice showed that lines with multiple superior haplotypes performed better (Abbai et al., 2019).Overall, haplotype-based breeding has a huge potential for developing tailor made rice varieties with improved yield, grain quality, and nutrition.The superior haplotype-specific markers would be utilized in order to combine them into single background as well as to choose the ideal combination of haplotypes controlling the particular phenotype.However, the robust validation of haplotypes, selection of the right combination of haplotypes, understanding the environmental impact, haplotype and environment interactions, pleiotropic effects, and interactions among haplotypes are prerequisite for the success of a haplotype-based breeding program (Abbai et al., 2019;Mayer et al., 2020).

CONCLUSION
MQTL analysis is a popular in silico approach to identify precise and consistent genomic regions for use in breeding and to prioritize CGs for functional validation.In the present study, MQTL analyses of 155 Zn QTLs using a high-density genetic map identified 57 MQTLs, with a significant reduction in QTLs (63.2%) and their CIs (80%).MQTLs were found to be enriched with metal homeostasis genes/gene families with diverse functions indicating their importance in metal uptake, transport, partitioning, and loading into the grains.
In silico differential expression and co-expression analyses revealed a complex network of genes involved in metal homeostasis.The CG haplotypes analyses using 3K-RGP showed variations with their predominance in certain subgroups, and superior haplotypes were identified using Zn data of a subset of aus lines.On an average superior haplotype contributed 2 mg/kg more Zn and a group of lines with five superior haplotypes performed better over the others.The MQTLs, CGs, and superior haplotypes identified in our study are useful for the efficient Zn biofortification of rice.

A C K N O W L E D G M E N T S
Financial support for this study was provided by IRRI-ICAR collaborative project.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interests.

D A T A AVA I L A B I L I T Y S T A T E M E N T
All data supporting the findings of this study are available within the paper and within its Supporting Information published online.

R E F E R E N C E S
population: DH, doubled haploid; RILs, recombinant inbred lines; BILs, backcross inbred lines; MAGIC, multiparent advanced generation intercross.b Population size.c Number of markers.d Type of markers: RFLP, restriction fragment length polymorphism; SSRs, simple sequence repeats; SNP, single-nucleotide polymorphism.e Range of phenotypic variance of the QTLs.f Type of rice: BR, brown rice; MR, milled rice.g Method of assessment of grain Zn: ICP-OES, inductively coupled plasma optical emission spectrometry; AAS, atomic absorption spectrometry; ICP-MS, inductively coupled plasma mass spectrometry; EDXRF, energy-dispersive X-ray fluorescence; ICAP, inductively coupled argon plasma spectrometer; ICPS, inductively coupled plasma analyzer; XRF, X-ray fluorescence.

F
I G U R E 1 Number of Zn quantitative trait loci (QTLs) projected and meta-quantitative trait loci (MQTLs) identified on each rice chromosome.
Details of identified meta-quantitative trait loci (MQTLs) for grain Zn.
Abbreviations: AIC, Akaike information criterion; CI, confidence interval.a Flanking markers of the MQTL.

F
I G U R E 3 Parts (a) and (b) represent the haplotypes identified across 3K rice genome panel (3K-RGP) for OsZIP7 and OsZIP5, respectively; parts (c) and (d) show the distribution of haplotypes in different subpopulations within 3K RG panel for the respective genes; and (e) circular bar chart showing frequency of haplotypes within various selected genes.

F
I G U R E 4 (a) Haplotypes distribution of various genes in aus subpopulation of 3K rice genome panel (3K-RGP); (b) haplotype distribution in 134 lines of aus subpopulation for the genes used for phenotypic validation; parts (c) and (d) represent the range of Zn content in 134 lines of aus subpopulation across three seasons for haplotypes of OsZIP7 and OsZIP5, respectively.Different alphabets denote statistically significant differences and vice versa from Tukey honestly significant difference test.T A B L E 4Haplotype-Phenotype relationship for the genes associated with grain Zn. number of accessions; -: absent in 134 selected lines; NV: not validated.Different alphabet denotes statistically significant differences and vice versa.