Molecular dissection of connected rice populations revealed important genomic regions for agronomic and biofortification traits

Breeding staple crops with increased micronutrient concentration is a sustainable approach to address micronutrient malnutrition. We carried out Multi-Cross QTL analysis and Inclusive Composite Interval Mapping for 11 agronomic, yield and biofortification traits using four connected RILs populations of rice. Overall, MC-156 QTLs were detected for agronomic (115) and biofortification (41) traits, which were higher in number but smaller in effects compared to single population analysis. The MC-QTL analysis was able to detect important QTLs viz: qZn5.2, qFe7.1, qGY10.1, qDF7.1, qPH1.1, qNT4.1, qPT4.1, qPL1.2, qTGW5.1, qGL3.1 , and qGW6.1 , which can be used in rice genomics assisted breeding. A major QTL (qZn5.2 ) for grain Zn concentration has been detected on chromosome 5 that accounted for 13% of R2. In all, 26 QTL clusters were identified on different chromosomes. qPH6.1 epistatically interacted with qZn5.1 and qGY6.2 . Most of QTLs were co-located with functionally related candidate genes indicating the accuracy of QTL mapping. The genomic region of qZn5.2 was co-located with putative genes such as OsZIP5, OsZIP9, and LOC_OS05G40490 that are involved in Zn uptake. These genes included polymorphic functional SNPs, and their promoter regions were enriched with cis-regulatory elements involved in plant growth and development, and biotic and abiotic stress tolerance. Major effect QTL identified for biofortification and agronomic traits can be utilized in breeding for Zn biofortified rice varieties.


Introduction
Globally, more than 3 billion people depend on rice for their daily caloric intake and nutritional needs (Yadaw et al., 2006;Elert, 2014;Bouis and Saltzman, 2017). However, milled rice is a poor source of micronutrients; hence, majority of the resourcepoor rice consumers without access to adequate nutrition suffer mineral deficiencies (Dipti et al., 2012;Goudia and Hash, 2015;Garcia-Oliveira et al., 2018;Goloran et al., 2019;Ludwig and Slamet-Loedin, 2019). Iron (Fe) and Zinc (Zn) malnutrition is common across all age groups, especially among children and women in the developing world (Bouis et al., 2002;Ahmed et al., 2016;Palanog et al., 2019;Gupta et al., 2020;Calayugan et al., 2021). Increasing the mineral density in the edible portion of the staple crops is a food-based approach to tackle mineral deficiencies. It is reported to be the most economical and sustainable solution to address malnutrition (Dipti et al., 2012;Bouis and Saltzman, 2017;Goloran et al., 2019;Ludwig and Slamet-Loedin, 2019). However, for the successful adoption and consumption of biofortified rice varieties, they should be highyielding, agronomically-superior, and desirable in terms of grain quality (Swamy et al., 2016;Swamy et al., 2021b).
Rice has a vast amount of genetic diversity available in cultivars, landraces, and wild relatives (Khush, 1997;Jackson and Lettington, 2003;Huggins et al., 2019). In particular, aus is a genetically-distinct group of rice accessions mainly originated from Bangladesh and India (Garis et al., 2005;Norton et al., 2018;Swamy et al., 2018). They are well-known for their wider adaptability, tolerance to biotic and abiotic stresses, and nutritional value (Xu et al., 2006;Ali et al., 2011;Henry et al., 2011;Gamuyao et al., 2012;Zhu et al., 2016). Currently, aus accessions are being explored for Zn biofortification due to their high grain Zn concentration Rakotondramana et al., 2022).
Quantitative trait loci (QTL) mapping using bi-parental populations is widely popular in rice (Al-Shugeairy et al., 2014;Travis et al., 2015). However, it has limited recombination and lesser allelic diversity, which results in larger confidence intervals of the QTLs, and captures only a portion of the total genetic variation (Xu et al., 2016;Verdeprado et al., 2018). Further, QTLs have to be fine-mapped and validated in multiple genetic backgrounds for their use in Marker Assisted Breeding (MAB). An alternative approach is QTL mapping in multi-parent populations that can increase the power of detecting major QTLs, and helps to understand QTL interactions and QTL by genetic background effects Descalsota et al., 2018). Multi-parent advanced generation intercross (MAGIC) and nested-association mapping (NAM) populations increase the probability of identifying precise QTL (Bandillo et al., 2013).
In this study, we developed four connected recombinant inbred lines (RILs) populations derived from Kaliboro (IRGC 77201-1) -a high grain Zn germplasm crossed with four elite Zn breeding lines. These populations were analyzed using multi-cross QTL (MC-QTL) method to identify QTLs for agronomic, yield and biofortification traits, candidate genes were shortlisted for major effect QTLs, and detailed analyses of major effect Zn QTLs was conducted to characterize candidate genes and cisregulatory elements.

Field establishment and phenotyping
All the mapping populations were grown in an alpha lattice design with two replications at Zeigler Experimental Station at IRRI during the dry and wet seasons of 2017 (DS2017 and WS2017) and in the dry season of 2018 (DS2018). Seeds were sown in the seedbed, and 21-day old seedlings were transplanted in the field at 20 cm x 20 cm planting distance. Each plot consisted of 2 rows with 20 hills per plot. Inorganic nitrogen (N), phosphorus (P), and potassium (K) fertilizers were applied at the rate of 120:30:30 NPK kg ha -1 during the DS and 90:30:30 NPK during the WS. Standard crop management practices were employed to ensure good crop growth.
We collected data on days to 50% flowering (DF), plant height (PH), number of tillers (NT), number of productive tillers (PT), panicle length (PL), thousand-grain weight (TGW), grain length (GL), grain width (GW) and grain yield per hectare (GY). We followed the standard evaluation system to gather data (IRRI, 2014). The Fe and Zn concentrations in the rice grains was measured using 3 grams of milled rice per sample by using XRF-Bruker S2 Ranger. Samples were analyzed twice, and mean values were considered for statistical analysis. and descriptive statistics of traits were performed using R v3.5.2 in R Studio v1.0.153 (RStudio Team, 2015). Broad-sense heritability (H 2 ) of all the traits was estimated using following formula: Where s 2 g is the genotypic variance, s 2 p is the phenotypic variance,

DNA extraction and genotyping
Fresh leaf samples were collected from RILs and parental lines during the early seedling stage. Samples were ground after freezing in liquid nitrogen and genomic DNA was extracted from each sample using the modified CTAB method (Murray and Thompson, 1980), and the quality of DNA was estimated by using 1% agarose gel electrophoresis. High-quality DNA samples with desirable concentration (~50 ng) were submitted to IRRI-Genotyping Services Laboratory (GSL) for SNP genotyping using the 7K Infinium Chip (Morales et al., 2020).

SNP genotypic data analysis
A set of 7,086 high-quality SNP markers were generated after thorough filtering using Illumina's protocol based on >80% call rate, homozygosity, and polymorphism between respective parents. SNP markers were further filtered based on homozygosity of the parents and with high heterozygosity on the segregating families. There were 495 high quality SNPs selected based on the distribution across the genome with < 20% missing SNPs among the parental lines. The selected SNPs were used to generate individual genetic maps and merged to construct a consensus physical map using BioMercator v4.2.1 considering a conversion of 250 kb = 1 cM. The scoring of alleles for the combined populations was performed by comparing parental alleles for each SNP. Four recipient parents were scored as "A", donor allele was scored as "B", and ambiguous/missing data were scored as "X".

Multi-cross QTL (MC-QTL) analysis
Combined phenotypic (BLUP) and genotypic data of four populations were used to perform MC-QTL analysis using MC-QTL v1.0 (Jourjon et al., 2005;De Oliviera et al., 2016). A linear regression model and iterative QTL mapping method (iQTLm) was used to detect QTL (Haley and Knott, 1992;Charcosset et al., 2000). A logarithm of odds (LOD) threshold of ≥5 was used to declare putative QTLs (Churchill and Doerge, 1994). The Phenotypic variance (R 2 ) for each QTL and global R 2 for each trait were generated. Epistatic analysis was also performed.

QTL detection using inclusive composite interval mapping (ICIM)
QTL detection for each individual bi-parental mapping population was conducted using inclusive interval mapping (ICIM) with the aid of IciMapping v.4.1 (Wang, 2009). Critical threshold value for QTL detection was calculated by 1000 random permutations of the phenotypic data to establish an experimentwise significance value at 0.05 (Churchill and Doerge, 1994). Estimated phenotypic variance explained by QTL for each trait and corresponding additive effect was also generated.

Candidate gene analysis
Consistent major effect QTLs with overlapping regions and with consensus boundaries of ≤500 Kb, found across different types of QTL analysis were subjected to candidate gene analysis. Predicted candidate genes were searched within or near ( ± 500 Kb) the QTL for each trait using its flanking markers. Physical locations of the annotated genes were determined with the aid of the RAP-DB In silico gene prediction and cis-regulatory elements analysis of grain Zn QTL SNPs within the candidate genes related to Zn homeostasis were downloaded from the rice SNP-Seek 18 (http://snp-seek.irri.org) (accessed on March 22, 2019). Coordinates of the genomic regions of interest were identified using the reference genome (Nipponbare) by extracting flanking sequences and aligning these to the target reference genome assemblies of IR64 and Kaliboro 600 (representing the genomes of the parents). A region in the target genome with best flank hits was then extracted and characterized. Alternative SNPs between indica-IR64 (representing the recipient parents) and aus-Kaliboro 600 (representing the donor parent) were identified using Rice SNP-Seek 18 (Mansueto et al., 2017). SNP genotyping data were sourced from 3K Rice Genome Project (The 3,000 rice genome project, 2014; Wang et al., 2018). The 1 Kb upstream region of each gene was used to shortlist the SNPs in their promoter regions. All polymorphic SNPs in the promoter and coding regions were listed.
Predicted genes showing SNP polymorphisms between the query genomes were further subjected to cis-regulatory analysis. The 1Kb region from the coding sequences in the 5' upstream strand of the Zn homeostasis gene was scanned for putative cis-regulatory elements using PlantPAN (http://PlantPAN.mbc.nctu.edu.tw) (accessed on March 23, 2019).

Phenotypic analysis of connected populations
The data of the various traits from four populations evaluated over three seasons are presented in Table 1. All the traits had large phenotypic variations and showed typical normal distribution across seasons and populations (Table 1 and Figures S1, S2). The highest mean values were recorded for GY, TGW, and GW during the DS; while DF, PH, and PL were higher during the WS. Grain Zn and Fe concentration were high during the DS. A low coefficient of variation (CV, <10%) was observed for DF and GL, whereas a moderate to high CV (10-45%) was observed for other traits. High broad sense heritability (H 2 ) was observed for Zn, DF, PH, and TGW (>70%) while it was low to moderate (10 to 60%) for all other traits (Table 1). Zn exhibited moderate to strong negative correlation with GY across populations ( Figures 1A-D). Likewise, it has moderate negative correlation with TGW and GL in Pop2 ( Figure 1B), and strong negative correlation with DF in Pop3 ( Figure 1C). Similarly, Fe was negatively correlated with several agronomic traits: PH, PL, DF, TGW, and GL in Pop1, Pop2 and Pop3 ( Figures 1A-C), while showed positive correlation with GW, NT, PT and GY in Pop4 ( Figure 1D). Meanwhile, a strong positive correlation between NT, PT, TGW, and GL; with PH and PL were

Multi-cross QTL analyses of connected populations
Overall, 495 SNPs across four populations were used to construct a consensus genetic map and used for MC-QTL analyses. The number of SNPs per chromosome ranged from 11 to 75 with lowest and highest number of SNPs on chromosomes 9 and 1, respectively ( Figure S5). The total map length was 1627 cM with an average distance of 6.0 cM between SNPs. MC-QTL analyses identified MC-156 QTLs for 11 traits ( Table 2). The global R 2 captured by the QTLs varied from 12 to 58% for various traits, while the R 2 explained by individual QTL ranged from 0.1 to 22%. The two major QTLs with highest R 2 (22% and 13%, respectively) were qPH 1.1 and qZn 5.2 . The tall parent (P2) allele       (Table 2). The qDF 7.1 located between id7005418 and c7p27670416 on chromosome 7 explained 7.3% R 2 , P1 allele increased maturity by 1.70 days while P5 allele reduced maturity by 0.80 days. Whereas, qNT 4.1 and qPT 4.1 were detected between 3647133 and 3653113 on chromosome 4, each explained by a R2 >3%. Meanwhile, qPL 1.2 detected on chromosome 1 explained 7% R 2 , and qTGW 5.1 identified between 5020568 and 5056493 increased TGW by 0.52 g. Among the 21 QTLs identified for GL, qGL 3.1 explained 8% R 2 , P1 allele increased GL by 0.09 mm. For grain width, qGW 6.1 contributed 4.2% R 2 and P3 allele increased grain width 0.03 mm. The qGY 10.1 explained 4% R 2 with highest positive allelic contribution from P2 (85.5 kg/ha). A total of 27 QTL was detected for Zn that explained R 2 of 0.3 to 13.0% with a global R 2 of 58%. As expected, the donor parent (P2) had the highest combined additive effect (2.6 ppm). Interestingly, the recipient parent (P5) total allele contribution was 2.5 ppm (Table 2). A major QTL (qZn 5.2 ) with an R 2 of 13% was detected on chromosome 5 at 22.03 Mb flanked by markers 5711540 and 5726844. It has a narrow genomic region of 600 Kb among the Zn QTLs. There were three QTLs, qZn 3.2 , qZn 11.1a and qZn 12.1b each contributed >3% R 2 . Similarly for Fe, 14 QTLs were detected with R 2 ranging 0.7-3.2% and a global R 2 of 25%. These QTL were located on all the chromosomes except chromosomes 5 and 10. The qFe 7.2 , located on chromosome 7 at 21 Mb, accounted for the highest R 2 (3.2%). However, the highest positive allele contribution with 0.16 ppm increase in Fe concentration observed for qFe 2.1 .

Epistasis
Epistatic interactions were detected between qZn 5.1 and qPH 6.1 , and qPH 6.1 and qGY 6.2 ; each explained a R 2 of 5% (Table S2). qPH 6.1 is a minor QTL flanked by id6007220 and 6228191 with genetic interval of 8.08 Mb. The positive allele (P3) increased PH by 1.25 cm while P4 allele reduced PH by 1.70 cm. On the other hand, qZn 5.1 is also a minor allele associated with Zn flanked by 514403 and SNP-5.9394949 with positive allele contributed by P2 (0.60 ppm). Minor QTL, qPH 6.1 is also interacted with qGY 6.2 which was flanked by markers 6907224 and id6016547 with a narrow interval of 0.07 Mb. The positive alleles were contributed by P3 and P4 with an additive effect of 100 kg/ha (Table 2).

Candidate genes for agronomic, yield and biofortification traits
Largest-effect QTL for grain yield, agronomic traits, Fe, and all the QTLs detected for grain Zn were subjected to candidate gene analysis (Table 2 and Figure S3). Genes conferring biological functions related to the trait were shortlisted (Table 3). The OsMADS18 is co-located with qFe 7.1 , it helps in Fe transport, cellular and inter-cellular responses under Fe deficiency, it is also known to regulate seed maturation and days to maturity. Similarly, HUA2 was linked to flowering locus qDF 7.1 and regulates flowering time and reproductive development. A pyruvate dehydrogenase kinase gene (OsPdk1) is found nearest to qPL 1.2 which is involved in root hair length and panicle threshability. Meanwhile, Glycogen Synthase Kinase (GSK) is found to influence days to maturity, seed maturation, brown rice shape, grain width and grain weight. GSK2 was co-located with major effect QTL, qTGW 5.1 , while GSK4 was co-located with qGW 6.1 . Similarly, GS3 was within the confidence interval of qGL 3.1 . Candidate genes OS01G0899425, OS04G0350700, and OS10G0326900 were co-located with qPH 1.1 , qPT 4.1 , and qGY 10.1 , respectively (Table 3).
Candidate gene prediction and cisregulatory analysis of qZn 5.2 The details of 20 candidate genes shortlisted from qZn 5.2 are presented in Table 4. Most of the candidate genes were either Zn finger or metal cation transporters. Four putative candidate genes identified within or nearest the QTL Viz: LOC_OS05G39540 (OsZIP9), LOC_OS05G39560 (OsZIP5), LOC_OS05G40490 and LOC_OS05G41790. These genes are differentially expressed in roots, leaves, stems, flowers, and meristems ( Figure S4). OsZIP5 was highly expressed in stem, internode and seeds, and moderately expressed in roots; OsZIP9 highly expressed in roots differentiation zone, and in inflorescence. LOC_OS05G40490 had minimum expression in leaves during drought stress, but LOC_OS05G41790 has high expression in seeds and had moderate expression in leaves under well-watered and drought conditions ( Figure S4). The in-silico gene prediction analyses of all these genes using the 3k genome identified 75 non-synonymous SNPs, 23 transversions, 34 transitions, and 18 deletions. Out of 20, 15 genes have non-synonymous SNPs found in their promoter and coding sequences (Table S3). The highest number of transitions (A/C, A/T, or G/C) was identified for OsZIP9 and LOC_OS05g40490, while OsZIP9 had the highest number of transversions (A/G or T/G). LOC_OS05g40490 had the highest number of deletions. The 15 predicted genes that show nonsynonymous SNPs were further investigated for cis-regulatory elements in their promoter regions using PlantPAN. A large number of cis-regulatory elements were detected and they were mainly associated with important physiological processes involved in submergence tolerance, light-regulation, meristematic tissue activities, and disease resistance (Figure 3).

Discussion
Breeding for improved nutrition has been prioritized in rice and other major staple food crops (Bouis and Saltzman, 2017;Gaikwad et al., 2020). There have been significant efforts over the last one decade to breed for high Fe, Zn and Vitamin A enriched crop varieties, and to deploy them on a large scale to create a health impact (Mwanga et al., 2021;Swamy et al., 2021a;Biswas et al., 2021). In rice, Zn biofortification breeding has been taken up on a large scale by IRRI and its partners leading to release of several high Zn rice varieties in Asia and Africa (Swamy et al., 2016;Tsakirpaloglou et al., 2019;Calayugan et al., 2021). Recently, Zn mainstreaming breeding has been initiated to incorporate Zn as a must trait in all the future rice varieties (Stangoulis and Knez, 2022). Identification of diverse donors, QTLs, genes, and a better understanding of molecular basis of grain Zn concentration, agronomic and yield traits are essential for the efficient Zn biofortification of rice (Swamy et al., 2016Calayugan et al., 2020;Rakotondramana et al., 2022) We used high Zn aus accession Kaliboro (IRGC 77201-1) to develop four RILs mapping populations. It has acceptable yield potential (~4900 kg/ha) and twice the amount of grain Zn (~40 ppm) in comparison to popular rice varieties (14-16ppm). In addition to that aus accessions are also genetically diverse and adaptable to a wide range of environmental conditions (Redoña and Mackill, 1996;Xu et al., 2006;Jagadish et al., 2008;Norton et al., 2010;Henry et al., 2011;Gamuyao et al., 2012;Al-Shugeairy et al., 2014). All the four recipient parents were high Zn breeding lines (12-19 ppm) developed at IRRI.
A wide range of variation for GY, Zn, Fe, and other major agronomic traits was observed in all the populations. Most of the traits exhibited normal distribution indicating the complex genetic basis of these traits (Table 1 and Figures S1, S2). Broad-sense heritability (H 2 ) was high for DF, PH, TGW, and Zn across populations (Table 1), which permits effective phenotypic selection for their improvement. These results are in consonance with earlier reports on traits distributions and heritability (Widodo et al., 2010;Du et al., 2013;Zhang et al., 2014;Swamy et al., 2016;Chang et al., 2018;Swamy et al., 2018). Consistent positive  (Jan et al., 2006) qTGW 5.1 5 6.59-7.50 GSK2 OS05G0207500 6.67 Average grain weight, protein concentration, total amylase activity trait, days to maturity, potassium uptake, seed maturation, seed weight, grain width  qGY 10.1 10 9.06-9.10 OS10G0325400 9.00 grain number, days to maturity, seed maturation, brown shape rice qGL 3.1 3 16.73-16.95

GS3 OS03G0407400
16.73 Grain length (Fan et al., 2006) qGW 6.1 6 20.14-20.98 GSK4 OS06G0547900 20.72 Grain width, days to maturity, seed maturation, brown rice shape, average grain weight (Yoo et al., 2006) Palanog et al. 10.3389/fpls.2023 Frontiers in Plant Science frontiersin.org correlations were observed between Fe and Zn ( Figure 1) indicating the possibility of simultaneous improvement (Swamy et al., 2016). However, concomitant breeding for high Zn and Fe with high GY will be challenging due to their strong negative correlations. Hence, appropriate breeding strategies would be essential for the successful Zn biofortification (Welch, 2004;Chang et al., 2018). One approach is to use high Zn donors with acceptable GY potential coupled with rapid generation advancement, speed breeding and genomic selection (Swamy et al., 2016;Calayugan et al., 2021;Swamy et al., 2021a;Yang et al., 2021). Moreover, breeding materials or populations that exhibit weak or no correlation between GY and Zn and transgressive seggregants that possess both high GY and Zn increases the chances of successful biofortification breeding (Figures 1, S1).  Cis-regulatory elements in the upstream region (500 Kb) of predicted genes for grain Zn. Various elements are associated with physiological functions such as submergence tolerance (blue), merismatic (orange), disease resistance (yellow), light regulation (green), amylase synthesis (red) were indicated using different colors.

MC-QTL analysis detected major QTLs and epistasis
Genetic analyses using multi-parent derived populations or combined analysis of multiple biparental populations that share common parentage help to dissect major effect and stable QTLs, which can work across genetic backgrounds without any epistatic or genetic background effects (Specht et al., 2001;Jannink and Jansen, 2001;Blanc et al., 2006;Doust et al., 2014). We carried out MC-QTL analyses for 11 traits using a set of 495 SNPs distributed across 12 chromosomes ( Figure S5). A high rate of SNPs polymorphism among parents was expected as Kaliboro (aus) and IR14M141, IR14M110, IR14M125, and IR95044:8-B-5-22-19-GBS (indica) are from genetically-distinct subgroups (Garis et al., 2005;Chen et al., 2014;Tang et al., 2016). The MC-QTL analysis detected more QTLs but effects were smaller compared to QTLs discovered by Inclusive Composite Interval Mapping (ICIM) (Tables 2 and S4). For instance, MC-QTL was able to detect 14 QTLs for Fe whereas only 2 QTLs were identified by ICIM. Some of the important QTLs detected using MC-QTL analysis viz: qZn 5.2 , qFe 7.1 , qGY 10.1 , qDF 7.1 , qPH 1.1 , qNT 4.1 , qPT 4.1 , qPL 1.2 , qTGW 5.1 , qGL 3.1 , and qGW 6.1 , which can be used in rice MAB and Genomic Selection (GS). Interestingly, most of the positive alleles for high GY and early DF were contributed by high yielding (IR14M110) and early maturing parents (IR95044:8-B-5-22-19-GBS) respectively. High-grain yield and early-maturity are two most desirable traits for rice varietal development (Zhao et al., 2011;MacKill and Khush, 2018;Li et al., 2019). Similarly, QTLs for PL and GW were contributed by Kaliboro, while TGW and GL were derived from IR14M141 (Table 2). Thus, it is evident that Kaliboro has contributed many positive alleles for yield and yield related traits and for improved grain Zn concentration. These favorable alleles can be pooled through genomics assisted breeding to develop rice varieties (Blanc et al., 2006;Grenier et al., 2015;Varshney et al., 2021;Xu et al., 2021;Reyes et al., 2022).
The QTLs of correlated traits (positive or negative) tend to cluster together on a chromosomal region Yu et al., 2021). The GY, TGW, GL, and GW were highly correlated, and their respective QTLs co-located on chromosome 5. All the QTLs for PT and NT, and Fe and Zn were co-located (Figures 1, 2). However, these co-locations may be due to pleiotropy or linkage or shared regulatory genes (Collard and MacKill, 2007;Du et al., 2013;Swamy et al., 2016;Descalsota et al., 2018). Positively correlated traits and their co-located QTLs can be readily mobilized into the breeding programs, while negative linkages must be removed by pre breeding, phenotypic selections, and fine mapping for their deployment in the rice breeding (Du et al., 2013;Swamy et al., 2016;Descalsota et al., 2018;Descalsota-Empleo et al., 2019a;Pradhan et al., 2020).
A notable result is that a major QTL for grain Zn (qZn 5.2 ) is consistently identified in all the populations and also in the combined population analysis. It had a R 2 of 13% in MC-QTL analysis and went up to 19.75% in ICIM (Table S4), with an additive effect ranging from 0.44 to 2.10 ppm (Tables 2 and S4), and had a narrow genetic interval of 600Kb. It is also devoid of any epistatic interaction or background effects, making it one of the potential candidate loci for MAB or GS to improve grain Zn concentration. This locus has been frequently reported from aus derived populations. Several studies have also identified similar genomic regions (Lu et al., 2008;Garcia-Oliveira et al., 2009;Zhang et al., 2014;Descalsota-Empleo et al., 2019b;Islam et al., 2020).
We observed very minimal epistasis among the QTLs. The PH (qPH 6.1 ) interacted with GY (qGY 6.1 ) and Zn (qZn 6.1 ), but their R 2 (<5%) and additive effects (0.12-1.25 cm, 9-282 kg/ha, 0.13-1.0 ppm, respectively) were low to make any significant changes in their phenotypic expression (Table S2). There were reports on genetic interactions between Zn and PH that affected their phenotypic performance, and also genetic background altering the epistasis effects leading to variable trait expression (Jiang et al., 2008;Lekklar et al., 2019). Therefore, epistatic QTLs and the influence of genetic backgrounds on major QTLs should be verified before their use in the breeding programs Islam et al., 2019;Lekklar et al., 2019).
The qZn 5.2 was also enriched with Zn/Metal homeostasis genes, especially Zn-finger proteins (Table 3). They play an important role in Zn uptake, transport and loading into the grains, and their expression provides abiotic stress tolerance under diverse climatic conditions (Goel et al., 2011;Zhang et al., 2014;Hara et al., 2017;Papierniak et al., 2018;Lekklar et al., 2019). It is noteworthy that these genes were also enriched with cis-regulatory elements in their promoter regions, and actively regulates physiological processes involved in metal homeostasis, submergence tolerance, meristematic tissue activities, light-regulation and disease resistance (Pandey et al., 2015;Biłas et al., 2020) (Table 4 and Figure 3).
Functional polymorphisms within genes play an important role in phenotypic diversity (Srivastava et al., 2014). We report that 15 of the 20 genes shortlisted from qZn 5.2 had functional polymorphic SNPs; notable ones are OsZIP9, OsZIP5, and LOC_OS05G40490. These transcription factors regulate expression of Zn homeostasis genes such as OsNAS1, OsNAS2, and OsNAS3, therefore play a role in metal transport, partitioning and loading into grains (Guerinot, 2000;Chandel et al., 2010). There are previous reports showing upregulation of OsZIP5 in roots and flag leaves, and OsZIP9 in roots of rice (Banerjee and Chandel, 2011;Chadha-Mohanty et al., 2015), but they are homologs and complement in phytosiderophoremediated Zn uptake (Lee et al., 2010). Recent haplotype analysis of OsZIP5 and OsZIP9 using a 3K rice panel revealed that superior haplotypes for these two genes were in high frequency in aus and indica sub groups (Gaurav et al., in press). On the other hand, there is less information on LOC_OS05G40490, a Zn knuckle domaincontaining protein gene involved in metal binding. Some of these genes can be pyramided using gene specific markers, explored for superior haplotypes, and can be used to develop high Zn transgenics or genome edited rice lines.

Conclusion
We carried out MC-QTL and ICIM analysis for agronomic, yield and grain micronutrient traits using connected rice populations. In all 156 MC-QTLs, 26 QTL co-locations and two epistatic interactions were detected. The donor parent, Kaliboro, contributed many positive alleles. A major QTL (qZn 5.2 ) for grain Zn concentration has been detected on chromosome 5 that accounted for 13% of R 2 . There were 20 genes within qZn 5.2 that are involved in metal binding, uptake, transport, partitioning and loading in different tissues. Three primary putative genes, OsZIP5, OsZIP9, and LOC_OS05G40490, along with 12 other genes showed functional SNP polymorphisms indicating their role in metal homeostasis. These genes were also enriched with cis-regulatory elements that regulate mineral uptake, growth and development, and tolerance to biotic and abiotic stresses etc. All the major effect QTLs can be readily mobilized into breeding programs through genomics assisted breeding. Promising candidate genes can be characterized to understand their functional roles. Over all our results provide insights into molecular bases of Zn homeostasis, and identified major QTLs/genes useful for Zn biofotification of rice.

Data availability statement
The phenotypic and genotypic data used in the analyses are provided as supplementary data in this article.

Author contributions
AP and BS designed the paper. AP led the conduct of the experiment and wrote the whole manuscript. CN, GD-E, MC, and ZS analyzed and gathered some genotypic and phenotypic data. AA and MI-A assisted in the establishment of field experiments and collection of phenotypic data. JH, PS, TB, AL, RM, KM, and BS thoroughly edited and contributed significantly to the improvement of the manuscript. All authors contributed to the article and approved the submitted version.

Funding
The work was supported by HarvestPlus (5402), Bill and Melinda Gates Foundation (OPP1194925), Department of Science and Technology (DOST) -Philippines. The funding entities did not influence the design, data collection, and interpretation, and writeup of the paper.