Identification of Genomic Regions Associated with High Grain Zn Content in Polished Rice Using Genotyping-by-Sequencing (GBS)

Globally, micronutrient (iron and zinc) enriched rice has been a sustainable and cost-effective solution to overcome malnutrition or hidden hunger. Understanding the genetic basis and identifying the genomic regions for grain zinc (Zn) across diverse genetic backgrounds is an important step to develop biofortified rice varieties. In this case, an RIL population (306 RILs) obtained from a cross between the high-yielding rice variety MTU1010 and the high-zinc rice variety Ranbir Basmati was utilized to pinpoint the genomic region(s) and QTL(s) responsible for grain zinc (Zn) content. A total of 2746 SNP markers spanning a genetic distance of 2445 cM were employed for quantitative trait loci (QTL) analysis, which resulted in the identification of 47 QTLs for mineral (Zn and Fe) and agronomic traits with 3.5–36.0% phenotypic variance explained (PVE) over the seasons. On Chr02, consistent QTLs for grain Zn polished (qZnPR.2.1) and Zn brown (qZnBR.2.2) were identified. On Chr09, two additional reliable QTLs for grain Zn brown (qZnBR.9.1 and qZnBR.9.2) were identified. The major-effect QTLs identified in this study were associated with few key genes related to Zn and Fe transporter activity. The genomic regions, candidate genes, and molecular markers associated with these major QTLs will be useful for genomic-assisted breeding for developing Zn-biofortified varieties.


Introduction
Micronutrient malnutrition or hidden hunger has become a major issue for people in developing countries, especially in Asia and Africa [1,2]. Marginal intake of micronutrients has been shown to contribute to increasing mortality rates, affect livelihoods, and result in adverse effects on millions of school-going children and pregnant women [3]. About three billion people worldwide are estimated to be affected and more than 24,000 people die daily owing to micronutrient malnutrition or hidden hunger [4]. Recently, the United Nations (UN) declared that tackling micronutrient deficiencies is one of the sustainable development goals (SDGs) and set a goal of SDG 2 to be achieved by 2035 (https://sustainabledevelopment.un.org/sdgs (20 September 2022)). Among different micronutrients, zinc (Zn) and iron (Fe) are the most essential and are associated with the development of physical and mental health in humans and reducing the risk of micronutrient malnutrition [5].
Microelements, such as iron (Fe), zinc (Zn), copper (Cu), and manganese (Mn), are found in many enzymes and have great importance for maintaining normal metabolic pathways in humans. Zn serves as a major co-factor for more than 300 enzymes and 2000 transcription factors and is involved in the metabolism of carbohydrates, lipids, proteins, and nucleic acids in humans [6,7]. It is the only metal present in all six enzyme classes (oxidoreductase, transferase, hydrolases, lyases, isomerases, and ligases) [8]. One third of the human population, particularly children and women, suffer from Zn deficiencyassociated health complications such as stunting, diarrhea, reduced immunity, poor cognitive development, and skin problems [9]. According to World Health Organization (WHO), Zn deficiency is sixth among the top 10 major causes of illness in developing nations, affecting 27-30% of the global population. A large portion of the global population suffers from mineral malnutrition, as they depend on plant-based diets that have a low mineral/calorific density [10].
Rice is a predominant and major staple food crop that provides the caloric needs of half of the global population [11,12]. A total of 483.3 million tons of milled rice is produced globally, according to a recent survey (ICRISAT., 2018). To meet the need, however, production must rise to 800 to 900 million tons by 2025 [13]. Mega varieties under cultivation are a poor source of micronutrients, particularly Zn in polished form (12 ± 14 ppm), and provide only one fifth of daily recommended Zn requirements (15 mg/day) [14,15]. In other words, 11 to 15 mg of Zn must be present in 220 gm of rice, which is the per capita amount of rice identified for India. The present international target as well as that from the AICRIP biofortification trial is 28 ppm, which can only meet half of the RDA. Through the AICRIP rice biofortification trial, medium slender, long slender, and short bold grain type rice varieties were released to some of the states in the country. Therefore, there is a need to develop biofortified varieties suitable for rice consumers under aromatic grain and other grain types. To ameliorate this situation, researchers are working to enhance the grain Zn content in rice by different approaches. Biofortification has emerged as a costeffective and sustainable strategy to tackle mineral deficiencies by enhancing micronutrient content without compromising the yield achieved through agronomic, conventional, and biotechnological approaches [16,17]. The availability of genetic variability in grain mineral nutrients in the rice germplasm can be employed to develop biofortified high-yielding rice varieties [18,19].
Understanding the genetic potential of the genotypes as well as interaction between genotype and environment is essential for developing successful biofortified rice varieties. Developing Zn-biofortified rice is challenging, as it involves complex genetics of the trait, genetic interactions such as epistasis, and environmental factors such as soil and water [14,[20][21][22]. Several genetic studies have also been carried out to identify markers linked with quantitative trait loci (QTLs) for high Zn in grains, which would help in the development of biofortified rice varieties through marker-assisted selection (MAS) [23]. QTL mapping provides information for identification of genomic regions associated with targeted traits of interest (Zn) by combining genetic information with phenotypic data [24]. Several QTLs with moderate to high phenotypic variance were reported for grain Zn on all 12 chromosomes of rice [25]. Recent studies on QTL mapping for mineral elements in rice using bi-parental and multi-parental mapping populations have identified multiple loci and demonstrated the genetic complexity of grain micronutrient traits [26].
Compared with other molecular markers (RFLP, RAPD, SSRs, etc.) single nucleotide polymorphisms (SNPs) are preferred because of their uniform distribution and wide occurrence across the genome (one marker every 100-500 bp, specific to species), thus making them the ideal choice for constructing high-density linkage maps and identification of markers closely associated with the trait of interest [27,28]. Next generation sequencing (NGS) technology with reduced cost provides abundant sequence information along with significant improvements in genome coverage and time [29]. Genotyping by sequencing (GBS), an advancement in NGS technology based on genome target reduction and restriction enzyme use has gained popularity as a cost-effective method in the development of genome-wide markers for genetic studies [30].
The present study aims to construct a high-density linkage map and identify QTLs for high grain Zn in polished rice and related agronomic traits using a recombinant inbred line (RIL) population derived from a MTU1010 X Ranbir Basmati cross. Putative candidate genes present in the QTLs identified for high grain Zn in polished rice were retrieved. The QTLs and putative candidate genes identified in the study can be deployed to advance genetic research and breeding applications for enhancing grain Zn content in rice.

Phenotypic Variation in MTU1010 X Ranbir Basmati RIL Population
Wide variations were observed among the RILs across the four seasons for the nine agronomic, yield-, and mineral-related traits, zinc (Zn) in brown rice (ZnBR; 8. In WS 17, G3 RIL showed Zn content above the threshold value of 24 ppm as well as SPY above 26 gm of yield check (IR64) (Figure 1a). In WS 20, G5 RIL showed Zn content above the threshold value (24 ppm) and G3 along with another seven RILs noted SPY. In the DS 18, G1 noted Zn content above the threshold value (24 ppm) and nine lines noted SPY above the yield check (26 gm). In DS 20, G6 noted Zn content above the threshold value and G3 line noted on-par yield with yield check. Based on mean Zn content and SPY of the four seasons, none of the genotypes showed Zn content above the threshold value, whereas seven RILs showed SPY ≥ the yield check (26 gm) ( Figure 1b).

Correlation Analysis
The relation between the mean values of individual traits varied from wet season to dry season in the correlation analysis. In the wet season of 2017-2020, positive correlations were observed between SPY and DFF at 0.05 or 0.01 probability. SPY showed significant negative correlation at 0.05 or 0.01 probability. Panicle length showed significant negative correlation at 0.05 probability with FeBR. NT showed significant negative correlation at 0.05 probability with ZnPR. In the dry season of 2018-2020, DFF showed significant negative correlation with FeBR and FePR at 0.001 and 0.05 probability, respectively. Significant positive correlations were observed among the four mineral traits, namely ZnBR, ZnPR, FeBR, and FePR, except for FeBR and FePR in the dry season of 2018-2020 (see Table 2). Significant positive correlations were observed among the four mineral traits ZnBR, ZnPR, FeBR, and FePR, except for FeBR and FePR in the dry season of 2018-2020.

Correlation Analysis
The relation between the mean values of individual traits varied from wet season dry season in the correlation analysis. In the wet season of 2017-2020, positive correlatio were observed between SPY and DFF at 0.05 or 0.01 probability. SPY showed signific negative correlation at 0.05 or 0.01 probability. Panicle length showed significant negat correlation at 0.05 probability with FeBR. NT showed significant negative correlation 0.05 probability with ZnPR. In the dry season of 2018-2020, DFF showed signific negative correlation with FeBR and FePR at 0.001 and 0.05 probability, respective Significant positive correlations were observed among the four mineral traits, nam ZnBR, ZnPR, FeBR, and FePR, except for FeBR and FePR in the dry season of 2018-20 (see Table 2). Significant positive correlations were observed among the four mineral tra ZnBR, ZnPR, FeBR, and FePR, except for FeBR and FePR in the dry season of 2018-202   The combined ANOVA across four seasons or environments indicating significant variance at probability < 0.001 was observed among the genotypes (Supplementary Table S5).

PCA
To explain the relationships among agronomic, yield-related, and mineral traits, PCA (Principal Component Analysis) with 92 RILs was performed using the mean values of all the four seasons.
The nine traits were resolved into nine principal components (PCs). The first four PCs noted an eigenvalue ≥ 1 ( Figure 2) with a cumulative proportion of 75.56% (see Table 2). ZnBR, ZnPR, and FeBR were significant in PC1. PH and PL were negatively significant in PC2; DFF, NT, and SPY were significant in PC3. DFF and NT were significant in PC4 ( Figure 2 and Table 3). The combined ANOVA across four seasons or environments indicating significant variance at probability < 0.001 was observed among the genotypes (Supplementary Table  S5).

PCA
To explain the relationships among agronomic, yield-related, and mineral traits, PCA (Principal Component Analysis) with 92 RILs was performed using the mean values of all the four seasons.
The nine traits were resolved into nine principal components (PCs). The first four PCs noted an eigenvalue ≥1 ( Figure 2) with a cumulative proportion of 75.56% (see Table  2). ZnBR, ZnPR, and FeBR were significant in PC1. PH and PL were negatively significant in PC2; DFF, NT, and SPY were significant in PC3. DFF and NT were significant in PC4 ( Figure 2 and Table 3).

High-Density Genetic Map
The genome-wide polymorphic SNPs were evaluated for expected segregation ratio using chi-square analysis in MTU1010 X Ranbir Basmati RIL population. Out of 32,24 SNPs between the parents, a total of 23,183 SNPs with contrasting alleles between th parents were found to be homozygous polymorphic and displayed segregation within th mapping population. A genetic map was constructed using 2746 SNPs spanning a geneti distance of 2445 cM ( Table 4). The total number of polymorphic SNPs mapped across 1 chromosomes ranged from 28 (Chr04) to 418 (Chr06), and the average genetic distanc ranged from 122 cM (Chr04) to 371 cM (Chr02) ( Table 4). A large variation in marke density was seen across the chromosomes, with the highest marker density on Chr12 (2.0 SNPs/cM) and lowest on Chr04 (0.23 SNPs/cM) ( Table 4).

High-Density Genetic Map
The genome-wide polymorphic SNPs were evaluated for expected segregation ratio using chi-square analysis in MTU1010 X Ranbir Basmati RIL population. Out of 32,245 SNPs between the parents, a total of 23,183 SNPs with contrasting alleles between the parents were found to be homozygous polymorphic and displayed segregation within the mapping population. A genetic map was constructed using 2746 SNPs spanning a genetic distance of 2445 cM ( Table 4). The total number of polymorphic SNPs mapped across 12 chromosomes ranged from 28 (Chr04) to 418 (Chr06), and the average genetic distance ranged from 122 cM (Chr04) to 371 cM (Chr02) ( Table 4). A large variation in marker density was seen across the chromosomes, with the highest marker density on Chr12 (2.05 SNPs/cM) and lowest on Chr04 (0.23 SNPs/cM) ( Table 4).

Analysis of Epistatic Interactions
A total of 170 epistatic interactions were detected for agronomic and yield (83) along with grain mineral (87) traits across four seasons (Supplementary Table S1). Eight epistatic interactions for five traits (PH, NT, FeBR, FePR, and ZnPR) were considered as the main effect QTLs, with PVE ranging from 11.20 to 26.83%. For two QTLs of ZnPR, one on Chr05 interacted with Chr10 (PVE 12.15%) and another QTL on Chr07 interacted with Chr11 (PVE 17.24%). Epistatic interaction for FeBR trait was observed on Chr04, which interacted with the locus on Chr05 (PVE 16.16%). Similarly, another epistatic interaction was observed for FePR on Chr02, which interacted with locus on Chr10 (PVE 26.83%). Two interactions were observed for NT on Chr01 with two loci of Chr02 and Chr08, and another QTL for NT on Chr02 interacted with loci on Chr05. Two digenic epistatic QTL interactions within the two loci of Chr02 were observed for PH, with PVE 8.74% and 24.49%, respectively ( Figure 5).

Candidate Genes Underlying QTLs for Zn QTLs
A total of 1637 genes present within the identified QTL region for Zn brown (ZnBR) and polished (ZnPR) on Chr02, Chr06, and Chr09 were retrieved (Supplementary Table  S2). They were annotated and categorized broadly on cellular, molecular, and biological processes using WEGO (Web Gene Ontology Plot) software (Supplementary Figure S3). Among several functional categories, 36 putative candidate genes (seven on Chr02, 18 on Chr06, and 11 on Chr09) belonging to transporter activity roles were shortlisted for the targeted QTLs (Supplementary Table S3). Out of 36 genes, two candidate genes on Chr06 showed a putative function of Zn and Fe transport; thus, the role of these candidate genes in Zn transport and metabolism was evaluated using knetminer software.
Os06g0705700 present in the QTL region of qZnPR.6.2 encoded by peptide transporter was found to be linked with 17 genes and three QTLs. Knetminer network analysis showed that this gene has a role in Zn and Fe transport (Figure 6a), and another gene Os06g0706100 present in the same QTL region on Chr06 is encoded by NRAMP4 and has a role in metal ion transport, particularly Zn transport. The two putative candidate genes (Os06g0705700 and Os06g0706100) within the identified QTL region are being

Candidate Genes Underlying QTLs for Zn QTLs
A total of 1637 genes present within the identified QTL region for Zn brown (ZnBR) and polished (ZnPR) on Chr02, Chr06, and Chr09 were retrieved (Supplementary Table S2). They were annotated and categorized broadly on cellular, molecular, and biological processes using WEGO (Web Gene Ontology Plot) software (Supplementary Figure S3). Among several functional categories, 36 putative candidate genes (seven on Chr02, 18 on Chr06, and 11 on Chr09) belonging to transporter activity roles were shortlisted for the targeted QTLs (Supplementary Table S3). Out of 36 genes, two candidate genes on Chr06 showed a putative function of Zn and Fe transport; thus, the role of these candidate genes in Zn transport and metabolism was evaluated using knetminer software.
Os06g0705700 present in the QTL region of qZnPR.6.2 encoded by peptide transporter was found to be linked with 17 genes and three QTLs. Knetminer network analysis showed that this gene has a role in Zn and Fe transport (Figure 6a), and another gene Os06g0706100 present in the same QTL region on Chr06 is encoded by NRAMP4 and has a role in metal ion transport, particularly Zn transport. The two putative candidate genes (Os06g0705700 and Os06g0706100) within the identified QTL region are being further investigated (Figure 6b).

Expression Analysis of Identified Candidate Genes
To predict the expression pattern of selected Zn-and Fe-related candidate g (Os06g0705700 and Os06g0706100), expression profiling was carried out using the parental lines MTU1010 and Ranbir Basmati, their derived inbred lines (I-42, I-44, an Figure 6. (a) Network analysis of candidate gene (Os06g0705700) in a major QTL qZnPR.6.2 using Knetminer analysis; (b) Network analysis of candidate gene (Os06g0706100) in a major QTL qZnPR.6.2 using Knetminer analysis.

Expression Analysis of Identified Candidate Genes
To predict the expression pattern of selected Zn-and Fe-related candidate genes (Os06g0705700 and Os06g0706100), expression profiling was carried out using the two parental lines MTU1010 and Ranbir Basmati, their derived inbred lines (I-42, I-44, and I-83), two mega rice varieties (IR64, BPT5204), a set of Zn-biofortified varieties (ZINCORICE, PAUSTIC1, PAUSTIC9) and a high-Zn landrace karuppunel. The recipient parent MTU1010, which has low Zn content, showed low expression of both the genes, whereas Ranbir Basmati showed higher expression levels (Supplementary Figure S4).

Discussion
One third of the global population is reported to be suffering from a lack of sufficient Zn nutrition [39]. Biofortification is considered as long-term sustainable strategy to address micronutrient malnutrition by enhancing the grain micronutrient density. Over the last decade, using conventional breeding approaches, a few biofortified rice varieties with high grain Zn have been developed and released globally, as Zn plays a critical role in several cellular and metabolic activities [40,41]. Though wide variability for grain Fe is available in brown rice, because of polishing, around 70 to 80% of grain Fe is lost [42].
Marker-assisted breeding has a vast potential to achieve desirable phenotypic variations in less time through the deployment of molecular markers linked to QTLs for desirable traits [43]. Large QTL regions identified with low-density SSR markers may have undesirable linkages, resulting in unsuccessful introgression [44]. In some of the recent studies, undesirable linkage groups were successfully eliminated by employing NGS-generated markers such as SNPs by identifying the recombinants. Next generation sequencing (NGS) technologies have become potential tools for the discovery of millions of SNPs in a cost-effective manner. The genotyping-by-sequencing (GBS) technique has facilitated the identification of key genomic regions for both complex and simple traits and has accelerated the breeding process [45][46][47][48]. ddRAD sequencing facilitates the identification of SNPs by reducing genome complexity in rice and other crops [49].
In the present study, we developed a high-density genetic linkage map from an RIL population MTU1010 X Ranbir Basmati using a GBS approach to identify QTLs for grain Zn and Fe using multi-seasonal phenotype data of the mapping population. Detailed analysis of phenotype data revealed a wide genetic variability within the RIL population for grain Zn and Fe along with agronomic and yield-related traits and showed normal distribution across four seasons. Previous studies indicated that complex inherited traits and micronutrient values exhibit wide variation, and these can be attributed to various genetic and environmental factors [50,51]. We examined data from the two wet and dry seasons in the current study and confirmed our results based on the previous variability studies in rice grain mineral content across the seasons (Kharif and Rabi) [52,53].

Correlation
Highly significant positive correlations were obtained among ZnBR, ZnPR, FeBR, and FePR, irrespective of season and location; however, there was a negative correlation with single plant yield (SPY). Several studies have reported positive correlation between Zn and Fe in rice [54][55][56][57]. Swamy et al. (2018b) [58] attributed the positive correlation between Zn and Fe contents to commonalities in their pathways and genetic networks of Zn and Fe uptake, translocation, and loading in rice. Negative correlations between yield and Zn in rice were earlier reported by [20,59] (G3: Mean Zn and Mean SPY). A QTL (qGZn9a) identified in Australian wild rice strain O.meridionalis was reported to be associated with an increase in grain Zn levels also found to be concomitant with fertility reduction [60]. The negative linkages between the yield and Zn must be eliminated for developing a successful biofortified variety by generating a greater number of recombinants per cross for identifying positive combinations for grain Zn and yield.

Identification of Promising RILs
Six stable RILs were identified for four mineral traits in the present study with >24 ppm (Zn) based on stability and G × E interaction analysis across four environments. From the study, we could also identify promising RILs with high grain Zn in polished rice (35.6 ppm) and promising RILs with single plant yield (28.4 g to 38.4 g), which could be further deployed in the breeding of biofortified rice varieties. For rice grain Zn and Fe, stability and G × E analyses are used for identification of stable donors [47,[61][62][63].
Considering the wide variability observed for the breeding lines with high grain Zn and Fe, stability and G × E analyses are being applied for selecting promising RILs in cereals. The contribution of environmental variation to grain Zn and Fe along with other agronomic traits in a RIL population of sorghum was demonstrated through genotype × environment interaction, correlation, and GGE biplot analyses [64]. Stable RILs with higher grain Zn and Fe content were also identified in the RILs of pearl millet using AMMI and GGE biplot analyses [65]. Different stable breeding lines were identified for different environments among eight Zn biofortified lines through stability and G × E analyses [51].
Several researchers have mapped QTLs for grain Zn and Fe traits in rice using various mapping populations, such as RILs, ILs, F2, DH, and MAGIC populations, and have identified multiple loci and clearly demonstrated the genetic complexity of the grain micronutrient traits. The genomic region on Chr02 (SNP_24248265-SNP_25682947) has shown co-localization in the qZPR.2.1 region [42].
Co-localization of QTLs for different element contents in grain is common in rice [70]. In the study, co-localization was observed: the genomic regions of Chr02, 05, and 06 of the DFF trait was co-localized with mineral trait QTLs FeBR, FePR, ZnBR, and ZnPR on Chr02. The Fe QTLs on Chr05 were co-located with ZnBR QTL, suggesting that there may be a possibility of selecting high grain Zn and Fe lines using molecular markers; Ref. [71] reported co-location of Zn and Fe QTLs on Chr02, 03, 08, and 12 in RIL populations.
In the present study, eight significant epistatic interactions were observed for PH, NT, FeBR, FePR, and ZnPR on Chr01, 02, 04, 05, and 07 and accounted for 11.20-26.83% PVE, suggesting the complex genetic variation of the traits. Two digenic interactions were detected for PH on Chr02, with PVE 8.74% and 24.49%, respectively. None of the identified digenic interactions in the present study were found to be involved with the main effect QTLs. Similar observations were earlier reported for epistatic interactions for grain Zn in rice [52,59]. Recent QTL mapping studies suggest that epistasis is considered to be a crucial genetic component underlying complex quantitative traits. Epistasis should be underscored in studying complex traits because it can account for hidden quantitative genetic variations in natural populations [72].
Integration of genomics with conventional breeding efforts can decipher the molecular mechanisms underlying traits of interest. Using the rice genome sequence, it is now possible to identify the putative candidate genes for grain Zn within QTL regions. In the present study, major QTLs involved in grain mineral element content that show consistency across seasons merit further examination by fine mapping and candidate gene analysis for ultimate use in MAS. The identified five candidate genes Os02t0629200 on Chr02 with Zn binding and transport activity, two genes Os09t0268300 and Os09t0297300 on Chr09 with transporter activity, and two genes Os06t0705700 and Os06t0706100 on Chr06 with Zn and Fe transporter activity and metal ion transport function could be promising genes for further validation. The enhanced expression of Os06t0705700 and Os06t0706100 in the flag leaf samples of the donor parent (Ranbir Basmati) and some of the donors/RILs with grain Zn further supported the identified QTLs/candidate genes in the present study.
Out of 92 RILs evaluated, one promising RIL G3 showed the highest SPY and Zn content across 92 RILs. The G3 RIL can be used as donor in biofortification breeding programs. Likewise, three other RILs, G1, G4 and G6, with high Zn content donors, and RILs G73, G49, G70, G82, and G52, with promising yield, can be used for future breeding programs.

Phenotypic Evaluation of Agronomic and Mineral Traits
Across four environments (seasons), five uniform plants were tagged from the center of each row, and phenotype data was collected for the traits. Days to 50% flowering (DFF), plant height (PH) (cm), number of tillers per plant (NT), panicle length (PL) (cm) at heading/early filling stage, single plant yield (SPY) (g), grain zinc (ZnBR, ZnPR), and iron (FeBR, FePR) in brown and polished rice were estimated at the post-harvest stage.

Estimation of Grain Zn and Fe Content
For estimation of grain Zn and Fe in brown rice (ZnBR, FeBR) and polished rice (ZnPR, FePR), seed samples were dehusked using a JLGJ4.5 rice husker (Jingjian Huayuan International Trade Co., Ltd., Hangzhou, Zhejiang, China). The brown rice was polished using a polisher with non-ferrous and non-zinc components (Krishi International India Ltd., New Delhi, India). After a thorough cleaning, samples of brown and polished rice (5 g) were subjected to energy-dispersive XRF (ED-XRF) (OXFORD Instruments X-Supreme 8000, Abingdon, UK) as per the standardized protocol [73].

Statistical Analysis
For each season, descriptive statistics such as mean, standard error of the mean (SEM), skewness, kurtosis, and coefficient of variations (CV%) were calculated for five plants using MS-Excel (2010). Trait-wise frequency distribution histograms were generated using R software [74]. Pearson correlation analysis and principal component analysis were computed using R script among nine traits. Different R packages (R core Team, 2018, Vienna, Austria), namely ggplot2, gge, FactoMinerR, and factoextra, were used to generate frequency distribution plots and individual comparison of RILs; the 'Metan' package was used for ranking of RILs [75].

Genotyping by Sequencing (ddRAD-Seq)
Genomic DNA was extracted from pooled young leaves (20-25 days) of both parents and RILs using the DNeasy plant kit (Qiagen, Hilden, Germany) following the manufacturer's protocol. DNA quality of each sample was checked on 0.8% agarose gel. Furthermore, DNA quantification was done using a Qubit 2.0 fluorometer (Themo Fisher Scientific Inc., Waltham, MA, USA). ddRAD-Seq protocol (modified GBS) was followed in the present study.
To perform GBS, genomic DNA was double digested with SphI and MlucI restriction enzymes (NEB, UK) and fractionated in 2% agarose gel to check the product size. The digested fragments were cleaned (Agencourt AMpure XP beads, Invitrogen, Waltham, MA, USA) using standard protocols. The ligation enzyme, T4 DNA ligase (NEB, England) was used to ligate the unique barcode adapters (4-8 nt sequence) at 16 • C for 30 min and heat-inactivated at 65 • C for 10 min.
This was followed by indexing with the addition of index-1 and index-2 (6-8 nt long) for multiplexing sequencing libraries in NGS Illumina. These libraries were PCR amplified (8-12 cycles) using a Phusion TM polymerase kit (Fisher Scientific, Loughborough, UK) followed by AMPure bead cleanup for purification to remove excess adapters and were sequenced on the Illumina HiSeqX platform (Illumina Inc, San Diego, CA, USA).

Genotyping and Filtration
The sequence reads for the parents and RILs were obtained as FASTQ files, which are used for SNP discovery. Raw reads were de-multiplexed according to their barcodes, and the adapters/barcode sequences were removed using standard software [76]. Highquality reads were aligned onto the rice reference genome of Nipponbare (MSU7) using bowtie2-2.2.6. The mapped reads were exported in the form of a Sequence Alignment Map (SAM) file by SAM tools, version 0.1.19 [77]. The alignment file was then processed for SNP calling using a Bayesian approach at specific site. Potential SNPs were filtered using the following criteria: loci with >70% missing data and those that showed distorted segregation of the two parental genotypes were excluded. SNPs with minimum allele frequency (MAF) > 0.05 and 90% call rate were considered for further analysis. The variant annotation was performed based on rice gene models, using in-house pipeline software developed by Agri Genome Labs Pvt., Ltd., Kochi, India.

Genetic Linkage Map Construction
High-quality SNPs obtained after filtering were used for map construction using the linkage mapping function using IciMapping v4.2. [78]. The grouping and ordering of 2746 SNP events between adjacent markers was performed at a minimum logarithm of odds (LOD) value of 2.5. 'Rippling' was done for fine-tuning of the ordered markers on their respective chromosomes. The Kosambi map function was used for the construction of the genetic map and to convert the recombination frequencies into map distances in centiMorgans (cM) [79].

QTL Analysis
The main effect QTL analysis was carried out using IciMapping. Deploying the BIP function in IciMapping, additive QTLs were identified by inclusive composite interval mapping (ICIM) based on a 1000-permutation test at a 95% confidence level. The QTLs with >3.0 LOD and phenotypic variance explained (PVE) > 10% were considered as major effect QTLs, those with PVE <10% were considered as minor effect QTLs for a particular trait. Epistatic interactions with the logarithm of odds (LOD) threshold value at 5.0 were considered as significant epistatic QTLs. The QTLs were visualized using MapChart V2.3 [80].

Confirmation of Identified QTLs and Markers
Out of mined candidates from QTLs, primers were designed for two candidate genes with known function of Zn transport, namely Os06g0705700 (CS528-F primer: GTGACG-GCTTCGGATGAG and R primer: CCCGGTGTAGAAGAAGGTAATG), Os06g0706100 (CS533-F primer: CCAATGCCTGGCCTACTT and R primer: CCACGTGGACACGTTCTT) using the Primer Quest tool (https://eu.idtdna.com/pages/tools/primerquest (2 September 2022)) with qPCR parameters (amplicon size < 150 bp; GC content > 50%; melting temp 60 • C; no secondary structures) and were synthesized at Integrated DNA Technologies (USA). Based on the RiceXpro data, flag leaf samples during the anthesis stage (after 10 days of flowering) from 11 genotypes (with differential grain Zn) were collected in RNALater™ (Sigma Aldrich, St. Louis, MO, USA). The published protocols for RNA isolation, cDNA synthesis, and checking of integrity of RNA and cDNA and their concentration were followed [83]. To investigate the expression of Os06g0705700 and Os06g0706100 genes, qPCR was performed with three biological replicates. Each qPCR reaction was performed in three technical replicates containing 5 µL SYBR ® qPCR Master Mix (Promega Corporation, Madison, WI, USA), 0.6 µM forward and reverse gene-specific primers, and 100 ng template cDNA. The qPCR was performed using QuantStudio5 Applied Biosystems (Thermo Fisher Scientific, Waltham, MA, USA) Real-Time PCR Detection System with the following conditions: 40 PCR cycles of denaturation at 95 • C for 15 s, amplification at 58 • C for 15 s, followed by extension at 72 • C for 10 s. The melting curve (melting-95 • C for 10 s, 65 • C for 1 min, 97 • C for 1 s) was analyzed for checking the amplicon specificity. Two internal control genes, Membrane Protein (Memp) (LOC_Os12g32950.1-F Primer: GAGCGC AAAGTTCCAGAAGAA and R Primer: CGCCACTAGTTGCCGTCCTGAT) and Tumor Protein Homolog (TPH) (LOC_Os11g43900.1-F Primer: CATTGGTGCCAACCCATC and R Primer: AAGGAGGTTGCTCCTGAAGA) were used for the normalization [84]. Relative fold change was calculated using the 2(−∆∆C(T)) method as described by [85].

Conclusions
The developed genetic map in the present study using MTU1010 X Ranbir Basmati cloud be a foundation for fine mapping of grain zinc (Zn) QTLs, and the identified genomic regions could be utilized in future breeding programs. The genes present in the identified QTL regions in this study reportedly play a key role in transporter activity. Further fine mapping, sequencing, and functional validation is required to identify the candidate genes in these QTL regions.