Genetic Analysis of Agronomic Traits and Grain Iron and Zinc Concentrations in a Doubled Haploid Population of Rice (Oryza sativa L.)

The development of micronutrient dense rice varieties with good agronomic traits is one of the sustainable and cost-effective approaches for reducing malnutrition. Identification of QTLs for high grain Fe and Zn, yield and yield components helps in precise and faster development of high Fe and Zn rice. We carried out a three-season evaluation using IR05F102 x IR69428 derived doubled-haploid population at IRRI. Inclusive composite interval mapping was carried out using SNP markers and Best Linear Unbiased Estimates of the phenotypic traits. A total of 23 QTLs were identified for eight agronomic traits and grain Fe and Zn concentration that explained 7.2 to 22.0% PV. A QTL by environment interaction analysis confirmed the stability of nine QTLs, including two QTLs for Zn on chromosomes 5 and 12. One epistatic interaction for plant height was significant with 28.4% PVE. Moreover, five QTLs were identified for Fe and Zn that harbor several candidate genes, e.g. OsZIP6 on QTL qZn5.1. A number of QTLs were associated with a combination of greater yield and increased grain Zn levels. These results are useful for development of new rice varieties with good agronomic traits and high grain Zn using MAS, and identification of genetic resources with the novel QTLs for grain Zn.

). For agronomic traits, DF ranged from 71.7 to 106.2 days with a mean of 86.9 days. PH ranged from 70.1 to 127.7 cm with a mean of 95.5 cm while YLD ranged from 2669.9 to 7409.1 kg ha −1 with a mean of 5446.1 kg ha −1 . GL ranged from 6.4 to 10.3 mm with a mean of 9.02 mm and GW ranged from 1.8 to 2.3 mm with a mean of 1.97 mm. In grain Fe and Zn concentrations, Fe ranged from 2.8 to 5.8 ppm with mean of 4.1 ppm while Zn ranged from 8.7 to 19.7 ppm with a mean of 12.6 ppm ( Table 1). The parent IR05F102 exhibited higher values for agronomic traits NT, NP, and YLD in all three seasons while IR69428 exhibited higher values for DF, PH, TGW, GW, Fe and Zn across three seasons and BLUEs. Coefficient of variation for agronomic traits ranged from 2.88% to 19.20%. Low CVs (<10%) were recorded for DF, PH, TGW, GL, and GW. Genotypic effects were highly significant and heritability values were high (0.52-0.98) for all agronomic traits and grain Zn concentration while lower heritabilities were observed for Fe, viz. 0.22 in S3 and 0.49 in S1 (Table 1). Pearson's correlation co-efficient between BLUEs and the three season's data showed significant positive correlations ranging between 0.52 and 0.97. Of the 45 possible correlations; 27 were significant in S1, with 12 positively correlated and 15 negatively correlated; while 26 were significant in S2, with 11 positively correlated and 15 negatively correlated; and 20 were significant in S3, of which 11 were positively correlated and nine were negatively correlated. A total of 29 correlations were significant using BLUEs, 13 of them were positive and 16 were negative  Table 2). In PC2, TGW and GW contributed most of the variation with loadings of 0.37 and 0.40, respectively.

Identification of QTLs.
A total of 23 QTLs for agronomic traits and grain Fe and Zn concentrations were identified in the DH population based on BLUEs (Table 2, Fig. 3). They were distributed on chromosomes 1, 2, 3, 5, 6, 7, 9, 10, and 12. The phenotypic variance explained (PVE) by the QTLs ranged from 7.18% to 22.03%. Two QTLs for DF (qDF 1.1 and qDF 3.1 ) were identified on chromosomes 1 and 3. These QTLs explained 14.82% and 12.63% PV. QTL qDF 3.1 was contributed by the high-yielding parent IR05F102. One QTL each was identified on chromosome 3 for PH, NT and NP, each explained more than 12% PV. Four QTLs for YLD (qYLD 1.1 , qYLD 3.1, qYLD 7.1 and qDF 12.1 ) were identified on chromosomes 1, 3, 7, and 12 contributed by IR05F102. Three of them QTL co-locations. Three QTL co-locations consisting of QTLs for two traits each were observed on chromosomes 3 and 9, two for agronomic traits and one for grain mineral elements (Fig. 3). Interestingly, these traits

QTL x environment interaction analysis.
To understand the stability and interaction of the identified QTLs with environment, QEI analysis was conducted. Estimated effects of 10 consistent QEI QTLs can be found in Table 3. One consistent QEI QTL was detected for DF on chromosome 1; two for YLD on chromosomes 3 and 7; two for GL on chromosomes 1 and 10; three for GW on chromosomes 3, 5 and 9; and two for Zn on chromosomes 9 and 12 (Fig. 4). All identified QEI QTLs were stable as indicated by the larger LOD A (1.55-9.67) than LOD AE (0.12-3.83) values and larger PVE A (1.93-13.91%) than PVE AE (0.05-3.21%) values (Table 3). epistatic interaction analysis. One epistatic interaction for PH was identified in the DH population (Supplementary Table 4; Fig. S2). The major di-genic epistatic interaction was identified between loci on chromosomes 3 with SNP marker interval id3495083 -3500757 Mb and on chromosome 7 with SNP marker interval id7862147 -7892971 Mb with a PVE value of 28.40%.
Candidate genes underlying QTLs for agronomic traits, Fe and Zn. Seven agronomic trait QTLs and five Fe and Zn QTLs were harbored in or near candidate genes (Table 4). For agronomic traits, two candidate genes (OsLFL1 and OsHD6) were identified for DF located on chromosomes 1 and 3. OsLFL1 with locus id Os01g0713600 was identified as a gene for flower development while OsHD6 with locus id Os03g0762000  www.nature.com/scientificreports www.nature.com/scientificreports/ was for days to heading. One candidate gene, OsLTS1 was identified for PH on chromosome 3 with locus id Os03g0837300. Four candidate genes were identified for YLD such as, OsDST, OsIPT4, OsMED5_3, and OsNAC139 on chromosomes 3, 7, and 12, respectively. Two candidate genes for GW, OsSRS3 and OsGS5 were identified on chromosome 5. Among QTLs for grain mineral element concentrations, for Fe, OsLysM-RLK10 and OsSWEET13 were candidate genes identified on chromosome 9 and 12, while candidate genes OsGATA8, OsSar1b, OsZIP6, and Os09g0511500 on chromosomes 1, 5, and 9 were identified for Zn (Supplementary Figs. S3 and S4). There were significant sequence differences in the candidate genic regions of the major QTLs (Supplementary Table S5).

QTL classes for YLD and Zn.
To understand the effects of combining the identified QTLs for YLD and Zn in this study, QTL classes were calculated. The mean YLD of DH lines without QTLs was 4430.9 kg ha −1 . On the other hand, the mean Zn of the DH lines without QTLs was 11.18 ppm. The genotype classes for YLD and Zn showed a wide range of 4129.4-6373.0 kg ha −1 and 11.1-17.7 ppm, respectively. An increase in YLD and Zn was observed with an increase in the number of QTLs (Tables 5 and S3). The IR05F102 allele at all four loci increased YLD. DH lines with three to four QTLs (genotype class 15 and 16) showed the highest YLD, 6373.0 kg ha −1 and 6285.7 kg ha −1 , respectively, among all YLD genotype classes. The IR69428 alleles at all four loci increased Zn. DH lines with four QTLs (qZn 1.1 + qZn 5.1 + qZn 9.1 + qZn 12.1 ) showed the highest Zn, 17.7 ppm, among all Zn genotype classes.

Identification of DH lines with good agronomic traits and high grain Zn. The top ten high Zn lines
were identified using BLUEs of the 148 DH lines (Table 6). IR 91153-AC 113-1 was early for DF and had the highest YLD (6697.23 kg ha −1 ) among the top lines. All four YLD QTLs (qYLD 1.1 , qYLD 3.1 , qYLD 7.1 , qYLD 12.1 ) and one Zn QTL (qZn 12.1 ) were detected in this line. IR 91153-AC726-1 had the highest value for Zn and Fe of 16.21 and 4.82 ppm, respectively. Two Zn QTLs were identified in this line (qZn 5.1 , qZn 12.1 ) and one YLD QTL (qYLD 1.1 ). These DH lines can be used as donors in breeding programs or can be directly tested in multi-location trials to further evaluate their performance and potential to be released as high Zn rice varieties.

Discussion
Compared to other cereals, rice is a poor source of essential micronutrients to fulfill daily human nutritional requirements 33 . Zn biofortification is considered to be a major solution that appears to be the most sustainable and cost-effective approach for addressing micronutrient malnutrition or hidden hunger 34 . Knowledge of genetic variation in agronomic traits, grain mineral elements concentration such as Fe and Zn and genes underlying allelic variation is vital for rice biofortification breeding to fast track the development of rice varieties that can make a positive contribution to human health. In this study, a DH population of rice was used to map QTLs for www.nature.com/scientificreports www.nature.com/scientificreports/ agronomic traits, Fe and Zn concentrations with SNP markers and BLUEs. Epistatic interactions, QTL by environment interactions, QTL pyramiding effects and haplotype analysis of major QTLs for YLD and Zn were also examined.
The DH population showed wide variation in all traits studied which specifies the traits complex and polygenic nature. BLUEs for grain Fe showed relatively little variation while grain Zn ranged from 8.7 to 19.7 ppm. The variation observed in the DH population for agronomic traits and grain Zn is important and a prerequisite for QTL mapping 35 . The variation indicates that a proportion of the phenotypic variance can be attributed to genotypic variance and that the agronomic traits and grain Zn can be exploited through selective breeding 36 . A wide range of variation also displayed the role of genotype as well as the effects of environment on the expression of the studied traits 37 . Several prior studies have also reported significant genetic variation specifically for grain Zn concentration [38][39][40] which supports that selection for high levels of Zn coupled with advanced molecular marker technologies is a feasible approach.
High H 2 values (>0.50) for all agronomic traits and grain Zn concentration were observed in the DH population which is consistent with previous studies that also work with rice agronomic traits and grain micronutrient concentrations [41][42][43] . However, the heritability was much lower for grain Fe concentration. Reliability of early generation selection of highly heritable agronomic traits (DF, PH, TGW, GW, and GL) is highly feasible, which may result in significant response to selection 44 . The low heritability observed for Fe indicates significant influence of environment and early selection for Fe is not feasible.
The first two PCs explained 55.94% of the total variation observed in the DH population. In general, principal components with vector coefficients more than 0.3 irrespective of their direction of influence are considered important 45 . Component loadings for each principal component revealed agronomic traits, such as YLD, TGW, www.nature.com/scientificreports www.nature.com/scientificreports/ GL and GW were among the phenotypic traits contributing positive projections while NT, NP, Fe and Zn contributing negative projections in two principal components which explained 30.34% and 25.60% of variation respectively. Through PCA we identified a number of traits which are responsible for the observed genotypic variation in the DH population, thereby identifying the traits with greatest impact on the phenotype of DH lines, which is informative for the selection of lines in breeding programs.
Analysis of relationships among quantitative traits is important for assessing the feasibility of joint selection of two or more traits and the effect of selection for secondary traits in genetic gain for the primary trait 46 . A positive genetic correlation between two desirable traits makes the job of the plant breeder easy in improving both traits simultaneously 47 . Clear insight on genetic correlations of agronomic traits and grain Fe and Zn concentrations can help plant breeder devise a suitable breeding strategy to enhance micronutrient density in rice. In the present study, DF was negatively correlated with TGW in BLUEs and across seasons, while PH was positively correlated with YLD. The trait NT was positively correlated with NP and YLD, but negatively correlated with GW. Similarly, GL was negatively correlated with Fe, while Fe and Zn were positively correlated. Previous studies also reported a positive correlation between grain Fe and Zn concentrations in rice 15,16,39,40 which infers common mechanisms for uptake, translocation and loading of Fe and Zn, and suggests that grain Zn can also be used for grain Fe selection. Negative correlations between YLD and Zn have also been observed 39,40 , though it was not consistent, and these   www.nature.com/scientificreports www.nature.com/scientificreports/ highlight the need for biofortification breeding programs to give significant weight to YLD while selecting for high grain Zn genotypes.
Considering the agronomic traits such as DF, PH, and YLD, Fe and Zn, the present study identified a DH line for which YLD was 6697.23 kg ha −1 with a Zn content of 15.11 ppm (IR 91153-AC 113-1). The high Zn lines identified in this study could be used in biofortification breeding programs to improve the grain Zn levels in rice resulting in improved rice varieties that could directly impact human nutrition, especially in populations that heavily rely on rice-based diets.
In this study, seventeen QTLs were mapped for agronomic traits and six QTLs were identified for Fe and Zn. It was notable that all the QTLs for grain Fe and Zn concentrations as well as QTLs for GW were derived from IR69428. Meanwhile, IR05F102 contributed all the QTLs for agronomic traits PH, NP, NT, YLD, TGW, and GL. These QTLs are good candidates for further studies such as gene fine-mapping and cloning. Our QTL results corroborated with recent studies on QTL mapping for agronomic traits and grain Fe and Zn concentrations in rice using diverse populations such as RILs 26 , BILs 28 , BC 2 F 3 derived lines 19 , DH lines 15 , and MAGIC lines 40 that revealed multiple loci located on all rice chromosomes. Consequently, these QTLs were detected in different environments and genetic backgrounds. The above results clearly showed the complexity of agronomic traits and grain Zn concentration 14,15,42 . Since the present study used BLUE values which excluded environmental effects and magnify genetic effects, the QTLs identified in this study may well prove more useful in rice biofortification breeding programs than the QTLs identified using single environment phenotypic values.
Haplotypes for two stable and major effect QTLs for grain Zn were further examined (Fig. 5). For qZn 5   www.nature.com/scientificreports www.nature.com/scientificreports/ two loci. The CA haplotype is associated with high Zn. Similarly, the two haplotypes in qZn 12.1 were GC and AA. The latter was associated with high Zn. The stability of the markers in these loci has clearly shown its potential for use in biofortification breeding programs although further validation is still necessary.
Co-localized QTLs are important for simultaneous improvement of traits of interest such as grain Fe and Zn concentrations in rice from plant breeder point of view 48 . QTLs controlling correlated traits were usually mapped in the same or close chromosome regions 49 . Some previous studies have reported the phenomenon of QTL co-locations for rice micronutrient concentrations 15,19,23 . In this study, several co-located QTLs linked with correlated traits were identified on chromosomes 3 and 9. For instance in agronomic traits, qDF 3.1 and qPH 3.1 were co-localized on chromosome 3 in the same SNP marker interval 3500757-3501392. For grain Fe and Zn, qFe 9.1 and qZn 9.1 were co-localized in the same SNP marker interval 9809545-9819278. These co-located QTLs may be the result of pleiotropy or multiple tight linkage of genes controlling the traits 50 . These regions might be useful for development of high grain Zn rice with highly productive agronomic traits. The co-location of QTLs for Fe and Zn will be beneficial for their simultaneous improvement with Marker Assisted Breeding 15,19,23 .
Identification of environment-specific and stable QTLs having consistent genetic effects across a wide range of environments is of great importance in rice biofortification breeding 37,40 . Several QTLs for agronomic traits and grain Fe and Zn concentrations have been identified, but the positions, effect sizes and genetic effect directions of these QTLs were mostly genotype and environment specific 51 . Studies on QEI contribute to the effective use of marker-assisted selection (MAS) in biofortification breeding, better understanding of the genetic architecture of grain Zn and agronomic traits, and the interactions of the genotype and environment 52 . The QTLs that were only detected in one environment mostly showed poor stability. All the consistent QTLs identified in this study: one for DF; two for YLD, GL, and Zn; and three for GW based on QEI analysis confirms stability as indicated by higher LOD A and PVE A values as compared to LOD AE and PVE AE values. The environment plays only a minor role in the additive effects of these QTLs. For instance, the additive effects of Zn QTL qZn 12.1 were 0.38 in S1, whereas in S2 the additive effects were only 0.11.
Expression of the phenotype is a result of polygenes interacting to control complex traits 53 . Understanding epistatic interactions can be useful to target genes through marker-assisted selection strategies for the improvement of complex traits such as YLD and grain Zn. In this study, we detected one epistatic interaction on chromosomes 3 and 7 for PH. The total variation explained by epistatic loci for PH were significant suggesting that epistasis in the form of additive by additive interactions, plays a very important role in controlling PH.
Several QTLs harbor genes for agronomic traits and grain Fe and Zn concentration. In this study, genes for flowering development on the chromosome 1 and days to heading on chromosome 3 were co-located with QTLs identified for DF 54,55 . OsDST, OsIPT4, OsMED5_3, and OsNAC139 genes for grain number, 1000-seed weight, and grain width 56 co-located with QTLs identified for YLD on chromosomes 3, 7, and 12. The candidate genes for grain Fe and Zn concentrations are involved in uptake, translocation, and metal homeostasis within rice plants. A network analysis of genes Sar1a and SWEET13 were conducted to understand their interaction with other traits and genes through the Knetminer program ( Supplementary Figs. S3 and S4). Sar1a was positioned within 12.71-13.79 Mb on chromosome 1 for Zn while SWEET13 was positioned within 16.98-17.57 Mb on chromosome 12 for Fe. Six QTLs and 35 traits were linked to Sar1a, while 12 QTLs and 42 traits were linked to SWEET13. A QTL for Zn on chromosome 5 was linked to candidate gene OsZIP6, which translocates Fe and Zn and may be responsible for high concentrations of these micronutrients in the grains 57,58 These candidate genes which were tightly linked with QTLs identified are worthy of further investigation.
QTL pyramiding is successfully utilized in rice breeding to develop lines resistance/tolerance to biotic and abiotic stresses 59,60 . To evaluate the pyramiding effect of each QTL in the present study, we compared the means of YLD and Zn between genotype classes in all combinations. Although, the arithmetic mean of YLD of all the genotype classes were higher than that of class 1 (no QTLs), no significant differences were detected. However, the means of class 15 (6373.00 kg ha −1 ) and class 16 (6285.74 kg ha −1 ) were significantly higher than that of class 4 (4129.43 kg ha −1 ). On the other hand, the mean Zn concentration of classes 3 to 11 ranged from 11.84 to 15.73  www.nature.com/scientificreports www.nature.com/scientificreports/ ppm, which was higher than that of class 0, though there were no significant differences detected among these groups. A significant difference was detected between class 12 (17.70 ppm) and class 1 (11.18 ppm). Hence, it will be necessary to pyramid favorable alleles of Zn and YLD QTLs in popular variety background to accumulate a large number of desirable alleles.

conclusion
The DH population showed wide variation for agronomic traits, grain Fe and Zn. High heritability was obtained for DF, PH, and GL and moderate for TGW and Zn. Significant positive correlations between PH & YLD, YLD & GL, TGW & GW, and Fe & Zn were observed. Most of the QTL alleles associated with improved grain Zn and Fe were contributed by the donor parent IR69428. Six QTLs contributed more than 15% of the phenotypic variance; qNT 3.1, qNP 3.1, qYLD 3.1, qTGW6 .1, qGL 6.1, and qZn12 .1. Based on QEI analysis, all the consistent QTLs were relatively stable. The epistatic QTL analysis detected one significant di-genic interaction for PH with 28.40% PVE. QTLs for grain Fe and Zn concentrations were associated with the candidate genes invovlved in metal homeostasis. A positive correlation between YLD and Zn levels was found in several combinations of three to four QTLs. The high Zn lines identified in this study could be used in biofortification breeding program to improve the grain Zn levels in rice resulting in improved rice varieties that could directly impact human nutrition.  www.nature.com/scientificreports www.nature.com/scientificreports/ Phenotyping. The population was evaluated for eight agronomic traits and two grain mineral elements. The agronomic traits were measured following the standard evaluation system 61 , including days to 50% flowering (DF), plant height (PH), number of tillers (NT), number of panicles (NP), grain yield (YLD), thousand grain weight (TGW), grain length (GL), and grain width (GW). For grain Fe and Zn analysis, 50 g paddy samples were dehulled using a Satake dehuller and milled for one minute using a K-710 mini-lab rice polisher. Milled rice samples weighing at least 3 g representing each plot were analyzed using X-ray Fluorescence Spectrometry (XRF) (Oxford) 62 . Measurements were done twice per sample and was expressed in parts per million (ppm). The average reading per plot was considered for further statistical analysis.

Statistical analysis.
A summary of basic statistical parameters was generated using STAR v.2.0.1. Analysis of Variance (ANOVA) was performed on single-season data and combined season data using one-stage multi-environment analysis implemented in PBTools v1.4. Best linear unbiased estimates (BLUEs) were generated by setting genotype effects as fixed and season effects as random. The BLUEs were used to perform QTL analysis. Boxplots and Pearson ,s correlation coefficients between pairs of traits were estimated using R Core team 63 .
The model used for ANOVA was: where μ is the overall mean, α i is the effect of the i th genotype; r j is the effect of the j th season, l kj the effect of the k th replicate within the j th season, b ojk was the effect of the o th block at the j th season of the k th replicate and ε ijk the error. The genotypes were considered fixed while replicates and blocks within replicates were random. Broad-sense heritability (H 2 ) for each trait in each season was calculated as: Genotyping. Plant genomic DNAs were extracted from leaf tissue using the cetyl trimethylammonium bromide (CTAB) method 64 . Quality check of DNA samples was carried out in 1% agarose gel. We submitted ~50 ng of DNA from the complete set of the doubled-haploid mapping population, along with the parents for genotyping with 7 K SNP array technology at the Genotyping Services Laboratory (GSL) at IRRI. Scanned image calls and automatic allele calling were loaded in the Illumina Genome Studio data analysis version V2011.1.
Linkage mapping and QTL analysis. The linkage map of the DH population was constructed using 379 high quality SNP markers using IciMapping v4.1 65 . The distribution of SNP markers varied across chromosomes. The number of SNPs per chromosome ranged from 15 SNPs (chromosomes 8 and 9) to 57 SNPs (chromosome 3). The total length of the linkage map was 1,629.6 cM while the average interval length was 4.3 cM (Supplementary Fig. S5). The linkage map was created using Kosambi function 66 . For analysis of QTLs, the BLUEs of each line in the DH population were used. The permutation method was used to obtain an empirical threshold for claiming QTLs based on 1000 runs of randomly shuffling the trait values at 95% confidence level using the BIP function, whereas epistatic interactions were identified by setting the logarithm of odds (LOD) threshold value at 6.0 utilizing inclusive composite interval mapping (ICIM) model in QTL IciMapping ver.4.1. The MET function in ICIMapping was used to study QTL x Environment Interaction. QTLs were visualized using MapChart v.2.3 67 . candidate gene analysis. The physical position of each QTL was determined by the position of the flanking SNP markers and genes physically located within or near QTLs for agronomic and grain Fe and Zn traits were considered candidate genes. Annotated genes with functions related to agronomic traits, metal transport and homeostasis were compiled and the physical positions of annotated genes were determined using the RAP DB Genome Browser 68 (http://rapdb.dna.affrc.go.jp/viewer/gbrowse/ irgsp1). Annotation and functions attributed to different candidate genes were downloaded from Oryzabase 69 (https://shigen.nig.ac.jp/rice/oryzabase/gene). A network of the top hit candidate genes and QTLs was created using K-netminer program 70 (http://knetminer. rothamsted.ac.uk/Oryza_sativa/).