Association Mapping of Ferrous, Zinc, and Aluminum Tolerance at the Seedling Stage in Indica Rice using MAGIC Populations

Excessive amounts of metal are toxic and severely affect plant growth and development. Understanding the genetic control of metal tolerance is crucial to improve rice resistance to Fe, Zn, and Al toxicity. The multi-parent advanced generation inter-cross (MAGIC) populations were genotyped using a 55 K rice SNP array and screened at the seedling stage for Fe, Zn, and Al toxicity using a hydroponics system. Association analysis was conducted by implementing a mixed linear model (MLM) for each of the five MAGIC populations double cross DC1 (founders were SAGC-08, HHZ5-SAL9-Y3-Y1, BP1976B-2-3-7-TB-1-1, PR33282-B-8-1-1-1-1-1), double cross DC2 (founders of double cross were FFZ1, CT 16658-5-2-2SR-2-3-6MP, IR 68, IR 02A127), eight parents population 8way (founders were SAGC-08, HHZ5-SAL9-Y3-Y1, BP1976B-2-3-7-TB-1-1, PR33282-B-8-1-1-1-1-1, FFZ1, CT 16658-5-2-2SR-2-3-6MP, IR 68, IR 02A127), DC12 (DC1+DC2) and rice multi-parent recombinant inbred line population RMPRIL (DC1+DC2+8way). A total of 21, 30, and 21 QTL were identified for Fe, Zn, and Al toxicity tolerance, respectively. For multi tolerance (MT) as Fe, Zn, and Al tolerance-related traits, three genomic regions, MT1.1 (chr.1: 35.4–36.3 Mb), MT1.2 (chr.1: 35.4–36.3 Mb), and MT3.2 (chr.3: 35.4-36.2 Mb) harbored QTL. The chromosomal regions MT2.1 (chr.2: 2.4–2.8 Mb), MT2.2 (chr.2: 24.5–25.8 Mb), MT4 (chr.4: 1.2 Mb Mb), MT8.1 (chr.8: 0.7–0.9 Mb), and MT8.2 (chr.8: 2.2–2.4 Mb) harbored QTL for Fe and Zn tolerance, while MT2.3 (chr.2: 30.5–31.6 Mb), MT3.1 (chr.3: 12.5–12.8 Mb), and MT6 (chr.6: 2.0–3.0 Mb) possessed QTL for Al and Zn tolerance. The chromosomal region MT9.1 (chr.9: 14.2–14.7 Mb) possessed QTL for Fe and Al tolerance. A total of 11 QTL were detected across different MAGIC populations and 12 clustered regions were detected under different metal conditions, suggesting that these genomic regions might constitute valuable regions for further marker-assisted selection (MAS) in breeding programs.


INTRODUCTION
Trace amounts of metal ions are essential for plant growth and development; however high concentrations result in perturbations in physiological processes, and ultimately productivity. Rice is among the cereal crop able to accumulate high levels of metal, including Ferrous (Fe), Aluminum (Al), and Zinc (Zn). These metals are naturally present at very low levels in paddy soils, but long-term use of chemical fertilizers result to acidity and increased concentration of phytotoxic ion form. Low soil pH (<5.0) favors the production of soluble phytotoxic Al 3+ ion easily taken up by the root system inhibiting cell division. These results to poor root growth and therefore low ion and water uptake (Panda et al., 2009;Tanaka and Navasero, 2012). Moreover, the anaerobic nature and low redox potential of paddy soils results to the reduction of Fe to Fe 2+ , another soluble ion easily taken up by plants. Excessive amount of Fe 2+ catalyze the formation of reactive oxygen species (ROS) causing irreversible damage in the cells (Becker and Asch, 2005). Zn is taken up by plant in the form of Zn 2+ ion during early stages of the plant, which is highly phytotoxic. It was reported that Zn 2+ has a key role in photosynthetic system. Specifically, it interferes in the photochemical reaction of chloroplast as proven by chlorophyll degradation in lichens, thereby decreasing its photosynthetic activity (Rout and Das, 2003).
Remediation strategies for metal contaminated paddy fields via chemical, physical, or biological means are necessary. However, the available methods are not effective or practical to use due to high input and running costs as well as low efficiency. Breeding new varieties with low metal accumulation in the grain constitutes a cost effective and efficient method to reduce the risk of low rice productivity and improve food safety. Understanding the genetics of metal tolerance is crucial to developing metal tolerance rice varieties.
Quantitative Trait Loci (QTL) mapping is an effective means of dissecting the genetic factors underlying agronomic traits such as metal tolerance. A number of studies have reported QTL for tolerance to a variety of metal. With regards to Fe toxicity in rice, a total of 197 QTL have been reported (Wu et al., 1997(Wu et al., , 1998(Wu et al., , 2014Wan J. L. et al., 2003a;Wan J. M. et al., 2003b;Wan et al., 2005;Shimizu et al., 2005;Ouyang et al., 2007;Dufey et al., 2009Dufey et al., , 2012Dufey et al., , 2015Fukuda et al., 2012;Zhang et al., 2013;Zhao et al., 2013;Matthus et al., 2015;Ruengphayak et al., 2015;Liu et al., 2016). A total of four chromosomal regions (CR), including CR1 on chromosome 1 between markers RM246 and RM443; CR2 on chromosome 2 between markers RM526 and R758; CR3 on chromosome 3 between markers C515 and C25; and CR4 on chromosome 7 between markers R1245 and RM429, have been found to be involved in the resistance of rice to Fe toxicity (Dufey et al., 2014). However, no genes have yet been cloned. A total of 148 QTL were identified for Al tolerance in rice via linkage mapping using biparental crosses (Wu et al., 2000;Nguyen et al., 2001Nguyen et al., , 2002Nguyen et al., , 2003Ma et al., 2002;Mao et al., 2004;Xue et al., 2006aXue et al., ,b, 2007 and association mapping using natural populations (Famoso et al., 2011;Zhang et al., 2016). The association mapping of 383 rice accessions and linkage mapping populations of recombinant inbred lines (RIL) derived from IR64/Azucena and backcross inbred lines (BIL) derived from Nipponbare/Kasalath//Nipponbare were reported by Famoso et al. (2011). Al-tolerant QTL (Alt TRG 12.1) encompassing the ART1 locus on chromosome 12 exhibited a large effect (LOD = 7.85, R 2 = 19.3%) in a RIL population. Moreover, three regions corresponding to induced Al-sensitive rice mutants (ART1, STAR2, and Nrat1) were identified through biparental QTL mapping (Famoso et al., 2011). Association mapping for relative root elongation (RRE) was performed using a core collection of 150 accessions of rice landraces with the highest phenotypic variation (R 2 ) explained by significant associations of 20.03% (for PSM365) at 21.4 Mb on chromosome 11 (Zhang et al., 2016). A transcription factor, ART1 (Al resistance transcription factor 1), has been identified for Al tolerance in rice. ART1 regulates the internal and external detoxification of Al by affecting at least 30 genes Ma et al., 2014). It is clear that the exposure of the roots to Al triggers both the induction and expression of many Al-resistance genes in rice including OsNrat1, OsSTAR1/2, OsALS1, and OsFRDL4 (Huang et al., , 2012Xia et al., 2010;Yokosho et al., 2011). Three studies have reported 53 QTL for Zn tolerance using RIL populations (Dong et al., 2006;Zhang et al., 2013;Liu et al., 2016). Among of them, the major QTL qZNT-1 at interval marker XNpb93-C3029C on chromosome 1 explained 21.9% of the variance (Dong et al., 2006). Four QTL were detected using two independent backgrounds in reciprocal introgression populations, with QSh2b, QSh7, and QSdw5 simultaneously identified in both Teqing-ILs and Lemont-ILs backgrounds , while qZRRDW3 was detected in both MH63-ILs and 02428-ILs backgrounds . The expression of these QTL in different genetic backgrounds suggests that they might be widespread. Two QTL, QSdw5 at 17.3-19.5 Mb on chromosome 5  and qFRSDW11 at C11S49-C11S60 on chromosome 11  were expressed under Fe and Zn stress, suggesting that there is a genetic overlap in Fe and Zn toxicity tolerance at the seedling stage.
The genetic architecture of tolerance to Fe, Zn, and Al toxicity in rice appears complex and not yet fully understood. No major locus has been identified, fine-mapped, or cloned thus far. Two limitations of previous studies include that biparental populations cover only a small genetic variability and a small number of genetic markers, resulting in the detection of only a few QTL confined to few genetic backgrounds (Matthus et al., 2015). The use of MAGIC populations with high allelic and phenotypic diversity, in combination with high-density genotyping is an effective means of increasing the genetic mapping resolution for metal tolerance.
Effective phenotyping techniques are a prerequisite in QTL mapping for metal tolerance. Artificial hydroponic systems provide an effective method for metal tolerance screening, as environmental factors such as temperature and acidity of the culture solution are highly controlled (Marmiroli et al., 2011). Moreover, growth rate, leaf color, and the extent of plant injury are common phenotypic traits used to measure the metal tolerance at the seedling stage in rice. The greatest indicator of plant sensitivity to Fe and Zn toxicity is related to leaf symptoms (Audebert and Fofana, 2009). Leaf discoloration and the leaf bronzing index (LBI) are both used to measure the extent of Zn toxicity (Dong et al., 2006), whereas the LBI is used for Fe toxicity (Dufey et al., 2009).
In this study, we aimed to illuminate the genetic basis of tolerance to Fe, Zn, and Al toxicity at the seedling stage by screening highly diverse MAGIC populations in a hydroponic system and employing a 50 K single nucleotide polymorphism (SNP) array. A genome-wide association study (GWAS) was conducted using mixed linear models (MLMs) to determine the loci associated with metal tolerance. The results will be valuable in gene cloning and marker assisted selection (MAS)-based breeding for metal tolerance.

Metals Tolerance Screening
The screening for Al experiment was conducted under 200 µmol AlCl 3 and control conditions from July 15 in 2015 in a greenhouse with the day/night temperature of 30/25 • C at CAAS, Beijing, China. Rice seeds of 873 MAGIC lines and 8 parents were surface sterilized with 5% sodium hypochlorite, rinsed with water, pre-germinated for 48 h under dark condition and temperature of 30 • C. Experiment was conducted using randomized complete block design (RCBD) with two replicates: for one replicate 24 germinated uniform seeds of each line were randomly sown in 96 wells PCR plate (8 × 12). One line was sown in three different plates (8 × 3) that plates have perforated wells at the bottom to facilitate the roots to fully contact with the standard rice nutrient solution (Yoshida et al., 1976). The screening for Fe and Zn experiments were conducted under 300 mg/L Fe 2+ (as FeSO 4 × 7H 2 O), 300 mg/L ZnSO 4 and control conditions from September 10 in 2014 in a greenhouse (NG04-02) with the day/night temperature being about 30/25 • C at IRRI (International Rice Research Institute), Laguna, Philippines. Rice seeds including 8 parents and 873 RILs were soaked in demineralized water and germinated at 30 • C in the dark for 48 h. Subsequently, 56 uniform seeds per line were directly sown into perforated Styrofoam sheets covered with nylon net at the bottom placed on standard rice nutrient solution. Each hole contained two seeds of the same line or parent. An Augmented RCBD was adopted with six sets of incomplete blocks 8 parents being replicated six times. The 10 days old plants were exposed to Fe, Zn, and Al stress for further 20 days. The pH of the control and all treated groups were adjusted to 4.5 using 1 N NaOH or 1 N HCl, pH were maintained every day and solutions were changed after every 5 days.
The metal tolerance scores (MTS) were measured at 20 days after treatment according to the modified standard evaluation system (Dufey et al., 2012). The MTS for Fe and Zn indicates the severity of toxicity: 1 (Highly tolerance: Normal growth with no leaf symptoms), 3 (Tolerance: Nearly normal growth, but leaf tips or few leaves whitish and rolled), 5 (Moderately tolerant: Growth severely retarded; most leaves rolled; only a few are elongating), 7 (Susceptible: Complete cessation of growth; most leaves dry; some plants dying) and 9 (Highly susceptible: Almost all plants dead or dying).
Three plants of each line from each replicate were harvested for measuring the shoot length (SL) and root length (RL). The roots and shoots were put in craft paper envelop and dried at 60 • C for 72 h in an oven. The shoot dry weight (SDW) and root dry weight (RDW) were determined. The indexes of toxicity tolerance, relative shoot length (RSL), relative root length (RRL), relative shoot dry weight (RSDW) and relative root dry weight (RRDW) were calculated according to the following formula: relative trait value (%) = (trait value in treatment)/(trait value in control) × 100.

Genotyping with the 55 K SNP Array
A total of 873 RILs from the MAGIC population, plus the eight parents were genotyped. Approximately 200 ng of DNA from each inbred line was used for genotyping using the Affymetrix R GeneTitan R platform with the Affymetrix R Axiom R Rice Genotyping Array, conducted by the CapitalBio Technology (Beijing, China) according to the manufacturer's protocol. The raw signal CEL files were processed using the Axiom R genotyping best practices. A total of 881 plates passed Dish Quality Control (QC). The probe QC was then determined using samples that passed the QC and were classified into six major categories including "Poly High Resolution, " "Mono High Resolution, " "No Minor Homozygote, " "Off Target Variant, " "Call Rate Below Threshold, " and "Other." Finally, 39,066 high-quality (Poly High Resolution) SNPs displaying genetic diversity were selected from the common SNPs. In addition, a three-step filtering strategy was applied to select high quality SNPs for QTL mapping. Firstly, markers that were non-polymorphic among the parents were removed. Secondly, all heterozygous genotypes were set as "missing" and markers with more than 10% missing values were removed. Finally,  Table S1).

Statistical Analysis
Adjusted trait value for each RIL was obtained using the PBTools developed by IRRI (http://bbi.irri.org/). Association analysis was conducted using MLM implemented in TASSEL v5.2 (Bradbury et al., 2007). The model uses PCA and kinship to control population structure and the familiar relationship. Bonferronicorrected threshold probability based on individual tests was calculated to correct for multiple comparisons, using 1/N, where N is the number of individual trait-SNP combinations tested. Significant marker trait associations (MTAs) were identified based on probability level of 1.0 × 10 −4 . Peaks exhibiting significance threshold level within a physical distance of 1.0 Mb were delineated into a single QTL. A QTL explaining more than 10% of the phenotypic variation was considered a major QTL. QTLs detected for different metal stresses with an overlapping confidence interval of 1.0 Mb were defined as a QTL cluster. The R-qqman package was used for creating the Manhattan plot (Turner, 2014). All other statistical analyses were conducted in R 3.3.1 (R Development Core Team, 2016).

Phenotypic Variation
All the parental lines differed significantly for MTS under Fe and Zn stress ( Table 2). For Fe, the parental lines A (1.0) and B (2.0) had the lowest MTS, while G (4.5) and H (5.0) had the highest MTS among the eight parental lines. For Zn, the parental lines C (3.0) and G (3.0) were the most tolerant, while the parental line B (6.5) was the most sensitive. Transgressive segregations were observed in both directions for all the populations (Figure 1), though the mean scores were not significantly different among populations (

Trait Correlations
The trait correlations were in the same direction for all MAGIC populations with respect to all the trait pairs (

QTL for Fe Toxicity Tolerance
In the five MAGIC populations, a total of 21 QTL were mapped for MTS, SL, RL, SDW, RDW, RSL, and RRDW on all the chromosomes, with the exception of chromosome 10. Among these, two, one, eight, two, and 11 were detected in the DC1, DC2, 8way, DC12, and RMPRIL populations, respectively (Table 3; Supplementary Figure S1). For MTS, three putative QTL were detected on chromosome 2. qFeMTS2.1 and qFeScore2.2 were identified in RMPRIL, and explained 2.2 and 2.4% of the phenotypic variation, respectively. qFeMTS2.3 was detected in the 8way population and explained 3.9% of the phenotypic variation. For SL, seven putative QTL were detected on chromosomes 1, 3, 7, and 9. qFeSL1.1, qFeSL3.3, qFeSL7, and qFeSL9 were detected in the RMPRIL population and accounted for 1.8-2.0% of the phenotypic variation, whereas qFeSL3.1 and qFeSL3.2 were detected in the 8way population and explained 3.5 and 4.0% of the phenotypic variation, respectively. qFeSL1.2 was detected in the DC1, 8way, DC12, and RMPRIL populations and explained 3.6-18.9% of the phenotypic variation. For RL, three QTL were detected on chromosomes 6, 7, and 8. qFeRL6 and qFeRL8 were detected in the 8way population and explained 3.4 and 4.3% of the phenotypic variation. qFeRL7 was detected in DC12 and explained 3.6% of the phenotypic variation. For SDW, two QTL were detected on chromosomes 2 and 5. qFeSDW2 was detected in the DC1 population and accounted for 7.5% of the phenotypic variation. qFeSDW5 was detected in the RMPRIL population and explained 1.9% of the phenotypic variation. For RDW, only one QTL, qFeRDW5, located on chromosome 5 was detected in the 8way population for RDW and explained 4.4% of the phenotypic variation. For RSL, four QTL were detected on chromosomes 6, 8, 11, and 12. qFeRSL6 and qFeRSL11 were detected in RMPRIL and explained 2.7 and 1.8% of the phenotypic variation, respectively.
qFeRSL8 and qFeRSL12, which explained 3.8 and 8.1% of the phenotypic variation, were detected in both the 8way and DC2 populations. For RRDW, only one QTL located on chromosome 4 was detected in the RMPRIL population, which accounted for 1.9% of the phenotypic variation. One major QTL, qFeSL1.2, located at 38.4 Mb on chromosome 1 was detected for SL in all populations except in DC2. The QTL explained phenotypic variation of about 18.9% in the DC1 population and 12.8% in the DC12 population. The allele from parental line C increased SL with an average effect of 6.4 and 7.6 cm.

QTL for Zn Toxicity Tolerance
In the five MAGIC populations, a total of 30 QTL were mapped for nine traits on all chromosomes except 5 and 9.
No SNP markers were detected to be significantly associated in the DC2 population, while six, two, eight, and 23 QTL were detected in the DC1, 8way, DC12, and RMPRIL populations, respectively (Table 4; Supplementary Figure S2). For MTS, one QTL, qZnMTS7, was detected on chromosome 7 and explained    7.7, 4.2, and 2.0% of the phenotypic variation in the DC1, DC12, and RMPRIL populations, respectively. For SL, 12 putative QTL were detected on chromosomes 1, 2, 3, 4, 6, 8, and 12, of which six QTL (qZnSL1.1, qZnSL2.1, qZnSL3.1, qZnSL3.2, qZnSL8.1, and qZnSL8.2) were detected in RMPRIL and explained 1.8-2.8% of the phenotypic variation. Three QTL (qZnSL4.1, qZnSL4.2, and qZnSL12) were detected in DC12 and explained 3.9-5.1% of the phenotypic variation. qZnSL6 was detected in DC1, explaining 8.4% of the phenotypic variation, while qZnSL2.2 was detected on chromosome 2 and explained 2.7-8.8% of the phenotypic variation in DC1, DC12, and RMPRIL, respectively. qZnSL1.2 was associated in the DC1, 8way, DC12, and RMPRIL populations, and explained 6.2-11.8% of the phenotypic variation. For RL, four QTL were detected on chromosomes 2, 6, and 11, the qZnRL2 was detected in the DC1 population and explained 9.0% of the phenotypic variation. qZnRL6.1, qZnRL6.2, and qZnRL11 were detected in the RMPRIL population and explained 1.7-3.7% of the phenotypic variation. For SDW, five QTL were detected on chromosomes 1, 2, 3, 4, and 8. qZnSDW1, qZnSDW2, qZnSDW3, qZnSDW4, and qZnSDW8 were detected in the RMPRIL population and explained 3.6% of the phenotypic variation. For RDW, two QTL qZnRDW2.1 and qZnRDW2.2 were detected on chromosome 2 in the DC12 and RMPRIL populations, respectively, accounting for 3.7 and 2.3% of the phenotypic variation. For RSL, two QTL were detected on chromosomes 2 and 4. qZnRSL2 was detected in both the DC1 and RMPRIL populations and explained 7.8 and 2.4% of the phenotypic variation, respectively. qZnRSL4 was detected in DC12 and RMPRIL and was associated with 3.7 and 1.8% of the phenotypic variation, respectively. For RSDW, two QTL (qZnRSDW8 and qZnRSDW10) were detected on chromosomes 8 and 10 in the 8way and RMPRIL populations, explaining 4.1 and 1.9% of the phenotypic variation, respectively. For RRL, a single QTL, qZnRRL3, was detected in RMPRIL and explained 1.8% of the phenotypic variation. For RRDW, one QTL qZnRRDW4 was detected on chromosome 4 in the RMPRIL population, which explained 2.1% of the phenotypic variation. One major QTL qZnSL1.2 (chr.1: 38.4 Mb) was detected for SL in all populations except DC2, which explained 11.8 and 11.3% of the phenotypic variation in DC1 and DC12, respectively. The allele from parental line C increased SL with an average effect of 6.7 and 9.1 cm.  Zhang et al., 2013 qZnRSL2

QTL for Al Toxicity Tolerance
A total of 21 QTL were mapped for the eight traits on all chromosomes except chromosome 4. Among them, seven, three, three, five, and 10 QTL were detected in the DC1, DC2, 8way, DC12, and RMPRIL populations, respectively (Table 5; Supplementary Figure S3). For SL, three putative QTL were detected on chromosomes 1 and 3. qAlSL1.1 was detected in DC2 and explained 4.0% of the phenotypic variation, while qAlSL3 was detected in the RMPRIL population and explained 2.8% of the phenotypic variation. qAlSL1.2 was detected in the DC1, 8way, DC12, and RMPRIL populations and explained 6.2-15.0% of the phenotypic variation. For RL, seven QTL were detected on chromosomes 1, 6, and 12. qAlRL1.1 and qAlRL1.2 were detected in the DC2 population and both explained 0.1% of the phenotypic variation. qAlRL6, qAlRL12.1, and qAlRL12.2 were detected in RMPRIL, explaining 1.9, 1.9, and 2.1% of the phenotypic variation, respectively. qAlRL1.3 was detected in the DC1 and DC12 populations and explained 9.4 and 5.1% of the phenotypic variation, respectively. qAlRL1.4 was detected in the DC2 and RMPRIL populations and accounted for 0.1 and 2.2% of the phenotypic variation, respectively. For SDW, five QTL were detected on chromosomes 1, 5,  Chen et al., 2012 qAlSL3 Yamaji et al., 2009 Chen et al., 2012 Krill et al., 2010 Krill 6, 7, and 9. qAlSDW1 was detected in 8way and explained 4.2% of the phenotypic variation. qAlSDW6 was detected in the DC1 population and explained 9.2% of the phenotypic variation. qAlSDW7 and qAlSDW9 were detected in the RMPRIL population and accounted for 2.0 and 2.2% of the phenotypic variation, respectively. qAlSDW5 was detected in the DC1 and DC12 populations and explained 10.2 and 4.9% of the phenotypic variation, respectively. For RDW, qAlRDW3 was detected in the DC1 and DC12 populations and explained 8.2 and 4.4% of the phenotypic variation, respectively. qAlRDW6 was detected in the 8way population and explained 4.4% of the phenotypic variation. For RSL, RRL, RSDW, and RRDW, only one QTL was detected in the single population. qAlRSL11, qAlRRL2, qAlRSDW10, and qAlRRDW12 were detected in the 8way, RMPRIL, DC1, and RMPRIL populations on chromosomes 11, 2, 10, and 12 and explained 2.0-9.6% of the phenotypic variation.
Two major QTL were detected for SL and SDW. The major QTL qAlSL1.2 (chr.1: 38.4 Mb) was detected in all populations except in DC2, with a phenotypic variation of 15.0% for the DC1 population and 14.6% for the 8way population. The allele from parental line C increased SL with an average effect of 3.7 and 8.0 cm. Another major QTL qAlSDW5 (chr.5: 5.4 Mb) was detected in the DC1 and DC12 populations and explained 10.2 and 4.9% of the phenotypic variation, respectively. The allele from parental line A increased SDW with an average effect of 0.0037 g.

Mapping Power and Resolution of the Magic Populations
We were able to overcome the limitations of bi-parental populations, in which only two alleles are analyzed and the genetic recombination is limited. In this study by using multiple parents to produce mapping populations with high allelic and phenotypic diversity, in combination with high levels of recombination events brought about by several cycles of intermating. The efficacy of MAGIC populations for high resolution mapping in rice has been previously reported (Bandillo et al., 2013;Li et al., 2013;Meng et al., 2016a,b). The eight parents used for the MAGIC population in this study displayed obvious differences in agronomic traits and biotic/abiotic stresses, such as resistance to drought, salinity, ferrous, zinc, aluminum, cadmium, and bacterial blight ( Table 1). These populations possess a slightly higher genetic diversity than a population of 248 of IRRI's breeding lines (Meng et al., 2016a). Due to the multiple hybridizations and selfing used in developing the populations, the DC1, DC2, and 8way populations exhibited no clear population structure and thus the effect of population structure on the mapping results is negligible. These populations were also used to identify QTL for 14 yield-related traits in the dry season (DS) and wet season (WS) during 2014 at the headquarters of the IRRI (Meng et al., 2016b). The mapping resolution and power of the DC12 and RMPRIL populations was higher than the DC1, DC2, and 8way populations as a result of the larger population size of the former. Greater metal toxicity-tolerant QTL were identified in the RMPRIL population (44 QTL) in comparison to the DC12 (15 QTL), 8way (13 QTL), DC1 (15 QTL), and DC2 (four QTL) populations. A total of 11 QTL were detected across the different MAGIC populations, particularly qFeSL1.2, qZnSL1.2, and qAlSL1.2, which were identified in all the populations except DC2. The QTL corresponding to known genes were also identified in the RMPRIL (13 QTL), DC12 (7 QTL), 8way (seven QTL), DC1 (six QTL), and DC2 (one QTL) populations. This suggests that the combined populations possess substantially higher QTL detection power and precision brought by the large population size and high level of recombination (Coles et al., 2010;Steinhoff et al., 2011). Although the MAGIC populations exhibited no clear population structure, approximately 20 lines showed exceptionally high kinship, and 32 lines were slightly different compared to the other lines in the DC2 population (Meng et al., 2016a). These results explain the observed low detection power in DC2 compared with the DC1 and 8way populations. The use of DC12 (DC1 + DC2) and RMPRIL (DC1 + DC2 + 8way) improved the mapping power and resolution. In addition, the mapping resolution in previously reported studies was limited by low marker density (several hundreds), which only covered a limited number of chromosomal recombination events occurring in biparental crosses (Huang and Han, 2014). In the present study, the high density 55 K rice SNP array genotyping platform was used to increase the mapping resolution.

Comparison of Identified and Reported QTL
In the present study, a total of 21, 30, and 21 putative QTL were detected under treatment with Fe, Zn, and Al, respectively. With respect to Fe tolerance, qFeSL1 and qFeSL3, located on chromosomes 1 and 3, were mapped closely to the leaf bronzing score-based QTL identified by Wan et al. (2005). qFeRDW5 and qFeSDW5 on chromosome 5 co-located with the QTL detected by Dufey et al. (2015) for RDW and MTS. qFeSDW2 colocalized with the QTL detected for SDW and MTS by Shimizu et al. (2005), and qFeSL3 associated with SL (Chr.3: 16.6 Mb) detected in the 8way population was located 1 Mb away from the previously reported OsIRO3 gene (helix-loop-helix family transcription factor) (Zheng et al., 2010). The qFeRL8 gene associated with RL detected in 8way (Chr.8: 0.7 Mb) was located 0.5 Mb away from the IDEF1 gene (Iron Deficiency-responsive Element-binding Factor 1) (Inoue et al., 2009). The QTL peaks of qZnRDW2, qZnRDW2, and qZnSL4 on chromosomes 2 and 5 were previously reported by Zhang et al. (2013) using reciprocal advanced backcross introgression lines (cross between Lemont and Teqing) for RDW and SDW traits. The Zn tolerance-related QTL located with the genes responsible for the translocation of metal in plants. The qZnSL3 (Chr.8: 6.9 Mb) gene was located 0.8 Mb away from the OsFRDL1 gene (citrate transporter), which is required for the translocation of iron (Yokosho et al., 2009). The qZnRL6 (Chr.6: 30.6 Mb) gene was located 264 kb away from the OsHMA2 gene (heavy metal ATPase 2), a P1B-ATPase that is involved in the root-to-shoot translocation of Zn and Cd in rice (Satoh-Nagasawa et al., 2012;Takahashi et al., 2012).. The qZnSL8 (Chr.9: 5.1 Mb) gene was located 1.1 Mb away from the OsZIP4 gene; a zinc-regulated zinc transporter (Ishimaru et al., 2005). The qZnSL4 (Chr.4: 30.1 Mb) gene was positioned at 102 kb away from the OsZIP3 gene; a zinc transporter (Sasaki et al., 2015).

Multi-Trait QTL for Tolerance to Fe, Zn, and Al Toxicity
Thirty-three of the identified QTL in the present study occurred in 12 regions on chromosomes 1, 2, 3, 4, 6, 8, and 9, each of which harbor QTL for more than two metal conditions (Figure 3). Five chromosomal regions were repeatedly covered by overlapping QTL for Fe and Zn stress. This is believed to be caused by the strong correlation observed between Fe and Zn stress tolerance in all the populations (Supplementary Table S2). Hence, it is plausible that the QTL for Fe and Zn were identified in the same QTL cluster. In previous QTL analyses it has been observed that the QTL for significantly correlated traits are usually located in the same chromosomal region. Two Zn toxicity tolerancerelated QTL, QSdw2a and QSdw5, were mapped together with the Fe toxicity tolerance QTL . One Fe toxicity tolerance QTL (qFRSDW11) and three Zn toxicity tolerance QTL (qZRRDW11, qZRTDW11, and qZRSDW11-1) were previously mapped together . In studies of tetraploid and hexaploid wheat, a grain Zn QTL on chromosome 2B colocalized with a grain Fe QTL, suggesting the possibility of simultaneous improvement of Fe and Zn tolerance (Velu et al., 2017). OsIRT1 (iron-regulated metal transporter) is not only a functional metal transporter for iron, and also involve as transporter for Zn and Cd (Lee and An, 2009). OsYSL2 is a rice metal-NA transporter that is responsible for the phloem transport of iron and manganese, including the translocation of iron and manganese into the grain (Koike et al., 2004). Expression of a rice NAS gene, OsNAS3, led to increased tolerance to Fe and Zn deficiencies and to excess metal (Zn, Cu, and Ni) toxicities (Lee and An, 2009). OZT1 confers plant tolerance to Zn and Cd ions (Lan et al., 2013). P1B-ATPases (also known as Heavy Metal ATPases: HMAs) are energized by ATP hydrolysis, and translocate heavy metals (Zn, Co, Cu, Cd, and Pb) out of cytoplasm (to plasma membrane and into vacuole) and thus play important roles in their transport, compartmentalization and detoxification. OsHMA2 transporter is involved in rootto-shoot translocation of Zn and Cd in rice (Satoh-Nagasawa et al., 2012;Takahashi et al., 2012). OsHMA3 can reduce the toxicity of Cd to rice seedling and maintains Zn balance in the rice stem (Sasaki et al., 2014). OsHMA4 is a causal gene for quantitative trait locus controlling Cu accumulation in rice grain (Huang et al., 2016). Heavy metal P-type ATPase, OsHMA5 and OsHMA4, involves in xylem loading of Cu in rice (Deng et al., 2013;Huang et al., 2016). In addition, more overlapping regions were detected for Fe and Zn tolerance QTL compared to previous studies Liu et al., 2016). A total of three and six multi-trait QTL were identified under three metal conditions with extremely significant values, p = 7.4E-7 and 5.1E-15, and were closely positioned (MT1.1 for SL and RL; MT1.2 for SL, RL, and SDW) at 35.5 Mb and 38.4 Mb of chromosome 1, respectively. The direction of the effect of all the clustered QTL detected for different tolerance traits were always equivalent, which may be caused by the pleiotropic effects of the same gene. Genetic overlaps for tolerance to different stresses have been reported in various stress conditions, including salt and drought (Wang et al., 2012) and blast resistance and drought (Xiong and Yang, 2003). These two aforementioned regions were previously reported to be associated with plant height using the same MAGIC populations (Meng et al., 2016a,b) as well as seedling shoot length (Abe et al., 2012). The 13 promising lines with MTS = 1 were selected under both Fe and Zn conditions, they are all from 8way populations. In our previous study, we also found that 8way was more powerful than the DC1 and DC2 populations for QTL identification useful for breeding (Meng et al., 2016a). The 13 ILs had an average higher RRL than average of eight parents by 29.5% (27.6%) under the Fe conditions and by 27.6% (22.3%) under Zn stress, respectively. The 13 IL S could be used as resistant donors in for breeding. For qFeMTS2.1, all lines had the same favorable allele from parental line A, C, D, E, F and G. The results indicated that the MAGIC populations are ideal material for genetic study and marker-assisted breeding (MAB), showing a tight integration of genetic research and breeding application in rice. Genetically overlapping loci can be easily exploited in the development of cultivars with tolerance to multiple metal toxicities, as it simplifies the process of MAB and is also associated with reduced costs.

CONCLUSION
The present study has identified the QTL associated with the toxicity tolerance of rice to three essential metal (Fe, Zn, and Al). As resistance cannot be measured directly, several parameters such as MTS and decreased biomass production were selected as indicators. A total of 21, 30, and 21 QTL were detected for traits related to Fe, Zn, and Al toxicity tolerance, respectively. Twenty-seven of the 72 QTL are in close proximity of previously reported QTL/genes for metal toxicity tolerance. Furthermore, three clusters MT1.1 (chr.1: 35.4-36.3 Mb),  were identified under all three metal conditions. The regions of the end of chromosomes 1 and 3 may play important roles in ion homeostasis and the maintenance of high biomass when rice is grown under different metal stress conditions. These results provide an opportunity to develop novel metal tolerant rice varieties via MAS.

AUTHOR CONTRIBUTIONS
QQ and GY designed the experiment; LM, XZ, and KP performed all the phenotypic evaluation; LM, BW, and XZ performed analysis and interpretation of the data; LM and GY drafted the manuscript; all authors revised the paper and approved the final version to be published.