Genome-Wide Association Studies Reveal the Genetic Basis of Ionomic Variation in Rice[OPEN]

A comprehensive study exploring the genetic architecture of the variation in the rice ionome sheds light on the mechanisms controlling mineral nutrient and trace element composition in plants. Rice (Oryza sativa) is an important dietary source of both essential micronutrients and toxic trace elements for humans. The genetic basis underlying the variations in the mineral composition, the ionome, in rice remains largely unknown. Here, we describe a comprehensive study of the genetic architecture of the variation in the rice ionome performed using genome-wide association studies (GWAS) of the concentrations of 17 mineral elements in rice grain from a diverse panel of 529 accessions, each genotyped at ∼6.4 million single nucleotide polymorphism loci. We identified 72 loci associated with natural ionomic variations, 32 that are common across locations and 40 that are common within a single location. We identified candidate genes for 42 loci and provide evidence for the causal nature of three genes, the sodium transporter gene Os-HKT1;5 for sodium, Os-MOLYBDATE TRANSPORTER1;1 for molybdenum, and Grain number, plant height, and heading date7 for nitrogen. Comparison of GWAS data from rice versus Arabidopsis (Arabidopsis thaliana) also identified well-known as well as new candidates with potential for further characterization. Our study provides crucial insights into the genetic basis of ionomic variations in rice and serves as an important foundation for further studies on the genetic and molecular mechanisms controlling the rice ionome.


INTRODUCTION
Plants require at least 14 essential mineral nutrients and several beneficial elements for growth, development, and resistance to biotic and abiotic stresses (Marschner and Marschner, 2012). Large quantities of macronutrient fertilizers (nitrogen [N], phosphorus [P], and potassium [K]) are used in modern agriculture to increase crop yields, often resulting in adverse impacts on the environment (Withers and Lord, 2002). Increasing the use efficiencies of macronutrient fertilizers is critical for environmental quality and agricultural sustainability. Meanwhile, plants also take up nonessential and toxic elements from the soil, which may cause phytotoxicity or enter the food chain, posing a risk to human health (Clemens et al., 2002;Williams and Salt, 2009). Cultivated rice (Oryza sativa) is one of the most important crops, as it feeds about half of the world's human population. For people who consume rice as a staple food, rice is a major dietary source of both essential micronutrients (e.g., iron [Fe] and zinc [Zn]) and toxic elements (e.g., cadmium [Cd] and arsenic [As]) (White and Broadley, 2009;Zhao et al., 2010;Clemens and Ma, 2016). It has been estimated that up to two billion people worldwide suffer from Fe and Zn deficiencies, particularly in populations using cereals as their staple food (White and Broadley, 2009;Swamy et al., 2016;Trijatmiko et al., 2016). Enhancing the accumulation of essential micronutrients and reducing the concentrations of potentially toxic elements in rice grain are of fundamental importance for food quality and human health.
The composition of mineral nutrients and trace elements in plants, defined as the plant ionome (Lahner et al., 2003;Salt et al., 2008), is determined by genetic, environmental, and developmental factors as well as their interactions. Many genes responsible for the uptake, translocation, and storage of mineral elements in plants have been identified. Based on studies using biparental populations, a large number of quantitative trait loci (QTLs) for mineral concentrations in rice have been reported (Lu et al., 2008;Garcia-Oliveira et al., 2009;Norton et al., 2010Norton et al., , 2012aDu et al., 2013;Zhang et al., 2014;Mahender et al., 2016;Ohmori et al., 2016). To date, only a small number of the mineral QTLs have been finely mapped, leading to the identification of the causal genes, including the transporter genes Os-HKT1;5 (or SKC1) for sodium (Na), Os-HMA3 for Cd, NRT1.1B for N, and Os-HMA4 for copper (Cu) (Ren et al., 2005;Ueno et al., 2010;Miyadate et al., 2011;Hu et al., 2015;. However, the gene networks controlling mineral accumulation and homeostasis are complex and largely remain to be elucidated (Lahner et al., 2003;Sasaki et al., 2016). A pivotal task in plant ionomic research is to unravel the genetic basis underlying the variations of the ionome across the different accessions of a plant species and crop cultivars. This knowledge is essential not only for understanding how plants adapt to their mineral environment but also for potential exploitation of alleles for breeding crop cultivars with improved nutrient use efficiencies and better quality.
Genome-wide association study (GWAS) is a powerful tool for unraveling the molecular basis for phenotypic diversity. GWAS has the power to genetically map multiple traits and provides a complementary strategy to classical mapping using biparental synthetic recombinant populations for dissecting complex traits . GWAS of the ionomic profiles of Arabidopsis (Arabidopsis thaliana) accessions has led to the identifications of a number of loci as well as genes related to mineral accumulation, such as HKT1;1 for Na (Baxter et al., 2010;Segura et al., 2012), MOLYBDATE TRANSPORTER1 (MOT1) for molybdenum (Mo) (Shen et al., 2012;Forsberg et al., 2015), HMA3 for Cd , and HAC1 for As (Chao et al., 2014). Rice landraces have evolved from their wild progenitors and show both high genotypic and phenotypic diversity . It is possible, therefore, to perform GWAS to identify the genetic basis of much of this phenotypic variation in rice. A large number of traits in rice, including numerous agronomic characteristics and metabolites, have been studied using GWAS (Huang et al., 2011;Chen et al., 2014;Matsuda et al., 2015;Yano et al., 2016). Extensive heritable variation is known to exist in the rice grain ionome for 1,763 rice accessions of diverse geographic and genetic origins (Pinson et al., 2015), and GWAS has been attempted to analyze the concentrations of four minerals (As, Cu, Mo, and Zn) in rice grain (Norton et al., 2014), eight elements (Zn, Fe, manganese [Mn], Cu, P, calcium [Ca], K, and magnesium [Mg]) in brown rice (Nawaz et al., 2015), and on aluminum (Al) tolerance (Famoso et al., 2011). These studies show the promise of GWAS for examining ionomic traits, but they are limited by the low density of single nucleotide polymorphisms (SNPs) and bias caused by the severe population structure inherent in the diversity panel containing both indica and japonica rice cultivars.
A common problem encountered in studies of ionomic QTLs is the large effect of the environment, especially soil properties and conditions that affect mineral supply and availability (Pinson et al., 2015;, which may mask the genetic variation or render QTLs unreproducible. Therefore, it is important to conduct ionomic studies on rice diversity panels across different field conditions to help reveal both the constitutive and environment-specific QTLs. Here, we describe a comprehensive study of the rice ionome based on GWAS using 529 diverse rice accessions, each genotyped at ;6.4 million SNP loci. The rice diversity accessions were grown under five different field or nutrient conditions. Vegetative tissues at the heading and maturity stages and rice grains were analyzed for 17 mineral elements (N, P, K, Ca, Mg, Fe, Mn, Mo, boron [B], Cu, Zn, cobalt [Co], Na, Cd, As, lead [Pb], and chromium [Cr]) using high-throughput inductively coupled plasma mass spectrometry. We identified large numbers of significant marker-trait associations and multiple potentially causal candidate genes, three of which we validated further. We also compared rice and Arabidopsis GWAS data and identified some conserved features of the genetic architecture underlying ionomic variations across plant taxa. This study provides crucial insights into the genetic basis of natural ionomic variations in cultivated rice and important information for comparative ionomic mapping in other plant species.

Variations in the Ionome among 529 Cultivated Rice Accessions
To investigate the variation in the ionome in cultivated rice, a population of 529 cultivated rice accessions representing large variations in the geographical origins, the genetic diversity, and the usefulness for rice improvement was used in this study (Supplemental Data Set 1). As described in our previous works Xie et al., 2015;Zhao et al., 2015), this population has been sequenced using the Illumina HiSeq 2000 system, generating a total of 6.7 billion reads. We included the sequences from 950 rice varieties generated by Huang et al. (2011) and selected SNPs with missing rates of less than 20%. We obtained a total of 6,428,770 SNPs with a missing data rate (SNP call rate) of ;0.38% for use in this study. Comparisons between the imputed genotypes and relevant high-quality genome sequences in the database as well as our array-based genotypes showed a high accuracy (>99%) of the imputed genotypes . The SNPs and imputed genotypes can be queried at RiceVarMap (http://ricevarmap.ncpgr.cn/v1/)  The ionome is influenced by genetic and environmental variations as well as developmental stage . N and P are two major nutrients that have large impacts on plant growth and development. In this study, we first investigated the variation in the rice ionome among 529 accessions grown in three adjacent paddy fields in Wuhan, China, in the same year varying in the status of N and P supply, including a normal field (NF) receiving standard N and P fertilization as well as a low-N (LN) field and a low-P (LP) field where either N or P fertilizers have been withdrawn for the last decade. The ionomic compositions (N, P, K, Ca, Mg, Fe, Mn, Mo, B, Cu, Zn, Co, Na, Cd, As, Pb, and Cr) of the shoots at the heading stage, straw, and brown rice at maturity were determined.
We found that the scale of ionomic variation depends on both the element and the plant tissue, with the essential macronutrients (N, P, K, Ca, and Mg; relative standard deviations varying between 10 and 38%) showing smaller variations than micronutrients (Fe, Zn, Cu, Mn, and B; relative standard deviations varying between 20 and 191%) or nonessential trace elements (Cd, As, Co, Cr, and Pb; relative standard deviations varying between 28 and 196%) and the reproductive tissue (brown rice) showing smaller variations than the vegetative tissues (shoots or straws) (Supplemental Figure 1, Supplemental Data Set 2). The broad-sense heritability (H 2 ) of the 17 elements, averaged across different field conditions and tissue types, ranged from 0.31 to 0.88 (Supplemental Table 1), consistent with previous reports that a significant portion of ionomic variation has a genetic basis (Pinson et al., 2015). In particular, elements such as N, P, K, Mg, Mo, Zn, As, and Cd showed high average H 2 ($0.68). The ionomic compositions also were influenced by N and P treatments. Low N resulted in considerable decreases in the concentrations of most elements, while low P also led to lower concentrations of some elements (e.g., N, Na, Fe, As, Cd, and Pb). In general, LN had a larger impact on the rice ionome compared with LP, as revealed by multivariate analyses. Principal component analysis based on the shoot or straw ionome revealed a reasonable separation of accessions between NF and LN field conditions but not between NF and LP field conditions (Supplemental Figure 2). Similar to the principal component analysis results, hierarchical clustering analysis based on the shoot or straw ionome showed that accessions grown in NF or LN fields were generally clustered together, while such clustering was not observed between accessions from NF and LP fields (Supplemental Figure 3). Indica and japonica are two main subspecies of cultivated rice showing significant diversity in phenotypic variation. In general, japonica accessions accumulated significantly higher levels of most elements than indica accessions (Supplemental Table 2). Median concentrations of Na, Mn, Co, Zn, and Mo in japonica accessions were consistently higher than those of indica accessions in all three tissues under three field conditions, with Mo concentrations showing the largest difference (1.4-3.9 times) (Supplemental Table 2). By contrast, the median Cd concentration in brown rice of japonica accessions was considerably lower than that of indica accessions, which is in agreement with previous reports (Arao and Ae, 2003;Uraguchi and Fujiwara, 2013;Pinson et al., 2015). The finding that the median shoot Cd concentration at the heading stage was larger in japonica than in indica suggests a lower translocation of Cd from the shoots to the grains in the former. Figure 1 summarizes the influence of field condition, plant tissue, and subspecies (indica versus japonica) on the ionome. The type of plant tissue has the most distinctive effect on the ionome, with N, P, Mo, Cu, and Cr being distributed preferentially to brown rice and forming a cluster and the 12 other elements in a second cluster showing an opposite pattern. The brown rice ionome appears to be separated according to the subspecies, whereas the ionomes of shoot and straw are separated by the N and P status of the fields.
To investigate further the effects of environment on the rice ionome, we grew the rice population in Youxian, Hunan province, for two consecutive years. The field site at Youxian had a lower soil pH and a higher Cd concentration than the soil at Wuhan (Supplemental Table 3). Ionomic analysis was conducted on both the vegetative tissue and dehusked brown rice harvested at crop maturity. Similar to the results from the trials at Wuhan, most mineral elements showed relatively high H 2 across the two years at Youxian (Supplemental Table 4). The differences in the mineral concentrations between the indica and japonica subpopulations also were generally consistent between the Youxian and Wuhan sites (Supplemental Tables 2 and 5).

Genetic Basis of Variations in the Rice Ionome
Based on the ionomic and sequence data for the 529 rice accessions, we performed GWAS of the indica and japonica subpopulations separately, using both simple linear regression and a linear mixed model (LMM). Because LMM results in fewer false positives , the GWAS results based on LMM are presented here. P-values of 1.8 3 10 26 and 4.1 3 10 26 were set as the significance thresholds for indica and japonica, respectively, after Bonferroni correction. For each element, a series of SNPs in a defined region with P-values lower than the threshold were detected. On the basis of the median distance of the linkage disequilibrium (LD) decay in indica and japonica (93 and 171 kb, respectively) of this rice population , the significant SNPs within a 300-kb region are considered to represent a locus, and the SNP with the lowest P-value in a locus is defined as the SNP in the closest linkage to the causal gene. Based on this criterion, a total of 373 loci (or associations) were detected across subpopulations, tissues, and trials, among which 41 were detected in at least two field conditions for the same tissue of the same subpopulation at Wuhan, 14 loci were detected in the two consecutive years at Youxian (Supplemental Data Sets 3 and 4), while 32 loci were scanned across both Wuhan and Youxian (Supplemental Data Set 5). A total of 72 loci were detected in at least two of the five trials, 13 in at least three trials, and 2 in all five field trials (Figures 2 and 3, Supplemental Data Set 6, Supplemental Figure 4). Of the 17 mineral elements analyzed, all elements had at least one GWAS locus and 9 elements had at least 4 loci (Supplemental Data Set 6).
The 72 loci presented in Figure 3 are based on repeatability in at least two of the five trials. Another way to reduce the environmental effect is to perform GWAS using the best linear unbiased prediction (BLUP). BLUP-based GWAS produced 60 and 69 loci for the Wuhan and Youxian sites, respectively, with 7 loci being common at both sites (Supplemental Data Set 7). Among the 122 loci identified by the BLUP method, 33 were the same or very close to the loci identified based on the repeatability method. Moreover, 11 of the 13 loci detected in at least three of the five trials by the repeatability method also were detected by the BLUP method, suggesting a good agreement between the methods.

Functional Interpretations of the Ionome GWAS Results
The ionomic loci identified by GWAS provide important clues for understanding the genetic architecture of the observed variations in the rice ionome. To identify candidate genes responsible for each ionomic locus, we extracted all genes within 300 kb of the most significant SNPs and considered their annotations, functions of homologous genes, and distances from the peak SNPs (Supplemental Data Set 8, Supplemental Figure 5). By applying this approach, we obtained a list of genes that represent plausible candidates for the causal gene for each of the loci controlling elemental concentrations in rice (Supplemental Data Set 9).
We chose two loci for further investigation with the aim of identifying the causal genes. We found that Os-HKT1;5 is close (;14 kb upstream) to the most significant SNP (sf0111478828) of an Na-related locus, which was detected at both the Wuhan and Youxian locations ( Figures 4A and 4B, Supplemental Data Set 9). Os-HKT1;5, encoding an Na transporter, was isolated previously in rice by map-based cloning using parents of Nona Bokra (a salttolerant indica variety) and Koshihikari (a susceptible elite japonica variety) (Ren et al., 2005), suggesting that allelic variation in Os-HKT1;5 explains the phenotypic variation in salt tolerance in the two rice varieties. Furthermore, the study showed that the phenotypic variation is caused by the difference in the transport activity of Os-HKT1;5 (Ren et al., 2005), not by the difference in the Os-HKT1;5 expression level. Based on our sequence data, we found nine SNPs resulting in nine amino acid changes in the coding sequence of Os-HKT1;5 ( Figure 4C), including the four amino acids that were suggested to be responsible for the functional difference between the Nona Bokra and Koshihikari alleles at Os-HKT1;5 (Ren et al., 2005). We investigated the linkage relationship between these nine SNPs and the most significant SNP at this locus. Interestingly, only one SNP (sf0111461701) is in strong linkage with the peak SNP ( Figure 4D). This SNP displays a substitution at the 184th amino acid of Os-HKT1;5 (H184R), with the variation occurring mainly in the indica subspecies ( Figure 4E). The two allelic groups separated by the base type of this SNP showed significant differences in straw Na concentration under all three field conditions, with the T allele (corresponding to His at the 184th position of the amino acid) containing lower Na levels than the C allele (corresponding to Arg at the 184th position of the amino acid) ( Figures 4F to 4J). Because the nine SNPs mentioned above are not closely linked with each other, there are likely to be more than two types of Os-HKT1;5 alleles. According to these nine SNPs, there are 21 haplotypes of Os-HKT1;5 in all 1479 accessions (529 accessions used in this study and 950 accessions sequenced by Huang et al. [2011]) (Supplemental Figure 6). Both the haplotypes of HKT1;5 NB (Os-HKT1;5 from Nona Bokra, designated as Hap2) and HKT1;5 K (Os-HKT1;5 from Koshihikari, designated as Hap3) were found in the accessions (Supplemental Figure 6). In fact, only six haplotypes (Hap1 to Hap6) exist in more than 10 accessions among the 1479 accessions. We analyzed the Na concentrations of the six haplotypes in our collection. We found that Hap2, the only haplotype representing the His amino acid at the 184th position, showed the lowest Na concentration in straws compared with other haplotypes in each subpopulation ( Figure 5). These results suggest that Hap2 is the strongest allele for restricting Na accumulation in the shoots and for Na tolerance. It has been suggested that Os-HKT1;5 functions in the roots to recycle Na out of the xylem to limit the translocation of Na to the shoot. Loss or reduced function of Os-HKT1;5, therefore, leads to increased shoot Na concentration (Horie et al., 2009). Taken together, our results indicate that H184R is likely the key substitution that causes functional variation in Os-HKT1;5. Unlike its homologous gene HKT1;1 in Arabidopsis (Baxter et al., 2010), Hap2 of Os-HKT1;5 shows no obvious difference in the geographical distribution from other haplotypes (Supplemental Figure 7), perhaps because rice is domesticated and, therefore, the link between genotype and the environment is now lost. Interestingly, Os-HKT1;5 is relatively conserved in japonica, with most accessions (93.5%) possessing the weak allele Hap3 (Supplemental Figure 6). Only one japonica accession has the Hap2 allele. Therefore, the Hap2 allele of Os-HKT1;5 has the potential to enhance salt tolerance of japonica as well as indica cultivars without this allele. Furthermore, indica II accessions contained significantly lower Na concentrations in the shoots compared with other subpopulations with the same Os-HKT1;5 haplotype (Supplemental Figure 8), regardless of whether the Os-HKT1;5 haplotype is strong or weak, suggesting that the indica II accessions harbor additional genes involved in regulating Na accumulation.
Because SNP sf0800123053 is strongly associated with Mo concentration ( Figures 6A and 6B), we investigated the causal gene for the locus. SNP sf0800123053 is located 37 kb from LOC_Os08g01120 (encoding a putative sulfate transporter named Os-MOT1;1), which is the closest homolog of At-MOT1 (a major Mo transporter in Arabidopsis) (Tomatsu et al., 2007;Baxter et al., 2008). Analysis of GUS expression driven by the Os-MOT1;1 promoter revealed that Os-MOT1;1 was expressed mainly in the lateral roots, with weak expression in the leaf blade, leaf sheath, The first part of the title for each panel indicates the element type; the second part represents the location from which it was scanned; the third part shows the tissue type (for Wuhan, the first letter indicates the field type [C, N, and P indicate NF, LN field, and LP field, respectively]; the second and third letter represent the tissue type [HS, shoots at the heading stage; MS, shoots at the maturation stage; MG, brown rice]); the last part represents the subpopulation type (Ind, indica; Jap, japonica). Red arrows point to the loci repeatedly scanned in different field trials. culm, node, and grain (Supplemental Figure 9). We obtained a knockdown mutant of Os-MOT1;1 and found that the mutant accumulated significantly lower concentrations of Mo in both roots and shoots than wild-type plants (Supplemental Figure 10), supporting the suggestion that Os-MOT1;1 is the causal gene underlying this Mo locus. Within this locus, significant SNPs were found almost exclusively in indica subpopulations. One of these SNPs, sf0800085089, was found to land directly on the promoter region of Os-MOT1;1 ( Figure 6C). Moreover, accessions separated by SNP sf0800085089 in indica (hereafter called Ind-C and Ind-T) showed significant differences in straw Mo concentrations in all five field trials (Figures 6E to 6I, Supplemental File 1), suggesting that Os-MOT1;1 contributes to the variation of Mo concentration in indica. The number on the left side of each column represents the physical location (Mb) of each lead SNP. The letters to the right of the column represent the corresponding elements. Colors represent the locations from which the loci were scanned: blue, Wuhan; green, Youxian; magenta, both locations. Underlining indicates that the locus was scanned at Wuhan in all three fields or at Youxian in both fields. One and two asterisks indicate that the locus was scanned across locations and was repeatedly scanned within one and two locations, respectively. The Venn diagram shows the numbers of repeated loci within or across locations.
We then analyzed the nucleotide diversity of the Os-MOT1;1 genomic sequence in indica. Several SNPs or insertions/deletions (INDELs) were found to show a close linkage to SNP sf0800085089 (r 2 > 0.65) ( Figure 6D). One INDEL locating at the 59 untranslated region of Os-MOT1;1 with a 10-bp insertion in the Ind-C group (higher Mo group) and one SNP locating at the open reading frame of Os-MOT1;1 resulting in a one-amino acid change were found and confirmed by fully sequencing 20 randomly selected indica accessions ( Figure 6C). Additionally, the fully sequenced data revealed a number of new variation sites in the promoter region, including two insertions of 999-and 226-bp fragments and some deletions of short fragments ( Figures 6J and 6K), but no additional variation was found in the Os-MOT1;1 coding sequence. To test whether the variations in the promoter region affect the expression level of Os-MOT1;1 among the indica accessions, we measured Os-MOT1;1 transcript levels in the roots of 10 accessions each of the . Box plots represent the interquartile range, the thick line in the middle of each box represents the median, the whiskers represent 1.5 times the interquartile range, and the dots represent outlier points. The data are based on two (for Wuhan) or three (for Youxian) biological replicates. DW, dry weight.
Ind-C and Ind-T groups under both Mo-sufficient and Mo-deficient conditions by qRT-PCR. The Ind-C group with higher Mo concentration had significantly higher levels of Os-MOT1;1 transcript than the Ind-T group under both conditions ( Figures 6L and 6M). At least four variants were found in the Os-MOT1;1 protein among the 529 rice accessions ( Figure 7A). To investigate whether there are functional differences among these variants, we overexpressed the four Os-MOT1;1 types in rice plants (japonica cv Zhonghua 11) and selected six transgenic lines for each Os-MOT1;1 type, three with relatively high expression and three with relatively low expression, for further study ( Figure 7B). We found that overexpression of Os-MOT1;1 increased Mo concentrations in both roots and shoots in an expression level-dependent manner ( Figures 7C and 7D). The relationships between expression level and root Mo concentration were similar among the four types of Os-MOT1;1 overexpressors ( Figure 7E). For shoot Mo concentration, the type 2 overexpressing plants appeared to be more effective than the other three types, but there was no significant difference among types 1, 3, and 4, the latter two types corresponding to the accessions of Ind-T and Ind-C, respectively ( Figure 7F). Taken together, our evidence is consistent with the notions that Os-MOT1;1 is the causal gene for the Mo locus with the lead SNP sf0800123053 and that the variation in Mo concentration is caused by the allelic variation in the promoter region leading to variable expression of Os-MOT1;1.

Comparative GWAS for Ionomic Traits between Rice and Arabidopsis
It has been suggested that natural ionomic variation in different species tends to be controlled by genes from the same gene families . In recent years, GWAS has been performed extensively in Arabidopsis, and several significant ionomic loci have been identified using a set of ;349 wild accessions (Baxter et al., 2010;Chao et al., 2012Chao et al., , 2014Forsberg et al., 2015). Comparing the rice ionomic GWAS data described here with those obtained previously using Arabidopsis should allow us to explore more systematically to what degree the genetic architecture controlling ionomic variation is conserved across taxa between rice and Arabidopsis. This comparative GWAS approach also could provide further evidence in support of candidate genes. We reran the GWAS using the publicly available leaf ionomic data from the set of 349 Arabidopsis accessions with an updated fully imputed SNP data set with 10,707,430 biallelic SNPs using an accelerated mixed model . The ionomic data were composed of the leaf concentrations of Na, Ca, Mg, B, P, K, Mn, Fe, Co, Cu, As, Zn, Cd, and Mo. All significant SNPs [2log(P) > 5] and corresponding candidate genes were organized based on their P-values per phenotype, as listed in Supplemental Data Set 10. As shown previously (Baxter et al., 2010;Chao et al., 2012Chao et al., , 2014Forsberg et al., 2015), the leaf concentrations of Na, Cd, As, and Mo showed the strongest associations with causal genes that have been established as At-HKT1;1 (AT4G10310), At-HMA3 (AT4G30120), At-HAC1 (AT2G21045), and At-MOT1;1 (AT2G25680), respectively. We also found several additional strong SNPs that are highly associated with leaf Co or Zn concentrations (P < 10 29 ) (Supplemental Data Set 10). The most highly associated SNP for the variation in leaf Co concentration was at Chr5: 902186 [2log (P) = 29.08] within ;1 kb of the gene IRON-REGULATED2 (At-IRGT2/FPN2; AT5G03570). This gene was shown previously to control leaf Co Letters highlighted with a colored background are single-letter codes for the amino acid residues at the corresponding positions. Red and blue indicate amino acid (AA) residues the same as and different from that in Hap2, respectively. Box plots represent the interquartile range, the thick line in the middle of each box represents the median, the whiskers represent 1.5 times the interquartile range, and the dots represent outlier points. The data are based on two (for Wuhan) or three (for Youxian) biological replicates. The box in each row in the box plot corresponds to the haplotype (Hap) of the same row on the left side. Significant differences at P < 0.05 within each group are indicated by different letters (one-way ANOVA test). DW, dry weight; Ind_int, indica intermediate; Jap_int, japonica intermediate; TeJ, temperate japonica; TrJ, tropical japonica. concentration using linkage mapping in a biparental F2 population derived from the accessions of Col-0 and Ts-1 (collected in Tossa de Mar, Spain) (Morrissey et al., 2009). In addition, a number of SNPs with 2log(P) > 5 for the leaf concentrations of Ca, Mg, B, P, K, Mn, Fe, or Cu and candidate genes also were identified (Supplemental Data Set 10).
We then extracted all genes within 300 kb of the peak SNPs in the rice GWAS data set for each element to look for orthologous genes with candidates obtained from the Arabidopsis GWAS data using the PLAZA Comparative Genomics Platform (PLAZA v2.5) (Van Bel et al., 2012;Tomcal et al., 2013). Only five pairs of orthologous genes were obtained for associations with As, Mn, Mo, and Na (Supplemental Data Set 11). The rice genes Os-HKT1;5 and Os-MOT1;1 associated with variation in rice shoot Na and Mo, respectively, were found to be paired with At-HKT1;1 and At-MOT1, which were determined previously to control the variation in these same elements in Arabidopsis leaves. This provides strong evidence that the variations in Na and Mo levels are governed by the same genes across monocots and dicots. Therefore, the comparative GWAS approach also provides a useful and independent method to further validate the roles of candidate genes identified using GWAS.

The Rice Ionome Is Affected by Heading Date
Heading date is an important trait for the adaptation of rice plants to different growth environments (Jung and Müller, 2009). In general, rice accessions with later heading dates have a longer vegetative phase and, consequently, a larger yield potential. To test whether heading date influences the rice ionome, we analyzed the correlations between elemental concentrations and heading dates (Supplemental Figure 11). We found significant negative correlations between heading date and the concentrations of N, P, K, Ca, B, Cu, Fe, and As in some of the three plant tissues or field conditions at Wuhan. By contrast, there were significant positive correlations between heading date and Cd concentrations in all tissues from all three field conditions as well as with Pb and Mn concentrations in some tissues and field conditions. Negative and positive correlations between heading date and As and Cd, respectively, also were observed in a panel of 467 locally adapted rice cultivars in south China (Duan et al., 2017).
N is one of the most important nutrients limiting crop productivity, whereas the overuse of N fertilizers can cause serious environmental damage. Increasing N-use efficiency in crops, therefore, is an important goal to enhance agricultural sustainability. The locus defined by SNP sf0709172004 was found to be strongly associated with the concentrations of N in shoots at the heading stage, while the nearby SNP sf0709177919 defined the locus for heading date (Figures 8A and 8B). Furthermore, the locus associated with N concentrations was still detected when heading date was used as a covariate ( Figure 8C). Therefore, we selected this locus for further analysis. We found that the gene LOC_ Os07g15770, known as Ghd7 (Grain number, plant height, and heading date7) and having a large effect on rice heading date, is located at this locus (;17 kb away from the lead SNP) (Xue et al., 2008). Because N concentration is strongly correlated with heading date (Supplemental Figure 11), we suggest that Ghd7 is a likely candidate. According to the sequence data, there are 45 haplotypes of Ghd7 in all 529 accessions (Supplemental Data Set 12). Within the indica accessions, there are 5 main haplotypes (Hap1 to Hap5) ( Figures 8D and 8E, Supplemental Data Set 12), among which Hap2 showed significantly shorter heading date but higher N concentrations than the other haplotypes ( Figures 8F to  8I). To investigate the effect of Ghd7, we compared N concentrations and heading dates among four near-isogenic lines (NILs) developed previously in the genetic background of Zhenshan 97 with the introgressed segment around the Ghd7 region from different rice varieties (Xue et al., 2008). Three lines containing functional alleles of Ghd7 [NIL(mh7), NIL(nip7), and NIL(tq7)] contained significantly lower N concentrations in various tissues than the line NIL(zs7), carrying a compete deletion in the corresponding region (Figures 9A to 9C). Furthermore, overexpressing a functional allele (mh7) of Ghd7 under the control of the ubiquitin promoter in the Zhenshan 97 background generally decreased N concentrations in shoots, straw, and brown rice compared with the wild type ( Figures 9D to 9F, Supplemental Figure 12). The wild type was used for comparison because there was no significant difference between the wild type and the null segregants in the T1 generation for plant height, heading date, or straw N concentration.

Loci Involved in Rice Growth under Low-N or Low-P Conditions
It has been suggested that N uptake explained most of the variation for N-use efficiency and grain yield under low N supply (Hirel et al., 2007). Variation in plant biomass production under low N or P supply reflects the variation in the ability to tolerate low N or P levels, which are important traits for breeding rice cultivars that can produce acceptable levels of yield with reduced fertilizer inputs. In this study, we found large variations among rice accessions in both the aboveground biomass and grain yield under low-N or low-P conditions (Supplemental Figure 13). We performed GWAS on aboveground biomass at both the heading and mature stages and grain yield. In our collection, biomass at both the heading and mature stages correlated negatively with N concentrations but positively with the total amount of N in the plant (Supplemental Figure 14). Grain yield also correlated positively with total N levels in the plant, especially in the LN field. Loci that were associated with variation in biomass and grain yield in LN fields are likely to be involved in controlling rice growth under low-N fertilization. We explored the loci identified from at least two traits (for instance, biomass at the heading stage and grain yield) in the LN field, LP field, as well as NF in Wuhan (Supplemental Table 6). Meanwhile, loci that were identified in the same tissue and subpopulation in at least two fields for each trait also were obtained. Based on the above criteria, we obtained 11 loci responsible for biomass and grain yield (Supplemental Table 6). Of these loci, two are colocated with either Ghd7 or Days to heading8/Ghd8, which are known to control the variation in biomass and grain yield (Xue et al., 2008;Wei et al., 2010;Yan et al., 2011). Interestingly, we found that five loci were detected exclusively in the LN or NF, whereas only two loci were common in both fields, suggesting different genetic bases for controlling biomass production under LN and NF conditions. By contrast, the LP field shared all loci with NF or LN.
Furthermore, a locus defined by SNP sf1125356264 was found under all three field conditions, which may contribute to variations in biomass or grain yield under various N or P fertilization conditions.
Genetic Basis of the Effects of Low N and Low P on the Ionome As described above, the concentrations of many elements were greatly affected by the N or P status (Supplemental Figure 1, Supplemental Data Set 2). This may reflect the crosstalk between N or P and other elements and the adjustment of the ionome in response to low-N or low-P conditions. Moreover, as N and P are essential macronutrients, N or P deficiency may affect the (G) to (I) Box plots for shoot N concentrations of different Ghd7 haplotypes at the heading stage in the NF (G), LN field (H), and LP field (I) in Wuhan. Significant differences at P < 0.05 within each group are indicated by different letters (one way ANOVA test). DW, dry weight. Box plots represent the interquartile range, the thick line in the middle of the box represents the median, the whiskers represent 1.5 times the interquartile range, and the dots represent outlier points. The data are based on two biological replicates. transcription levels of many genes and alter root morphology (Hermans et al., 2006). It is likely that these effects are genotype dependent. Therefore, investigating the genetic basis of changing ionomic profiles under low-N or low-P conditions may provide some clues to understanding the responses of the ionome to the environment.
Based on the GWAS results of the rice ionome performed above, we extracted 32 loci associated with the concentrations of 11 elements (Supplemental Data Set 13), which were scanned across at least two tissues in one field condition. It is interesting that the majority of these loci (19 of 32) also were common in two or more field conditions. By contrast, 13 loci were specific to field conditions, including 7 for NF, 5 for LN field, and 1 for LP field. Identification of the causal genes underlying these loci would provide insights into how N and P affect the rice ionome and also would benefit rice breeding aimed at tailoring varieties to local conditions.

DISCUSSION
As the ionome is an essential component of all living systems, natural ionomic variation has important effects on evolution and adaptation (Lahner et al., 2003;. Over the last decade, numerous genes and gene networks have been shown to be involved in controlling mineral nutrient and trace element homeostasis in plants. However, the loci that control natural ionomic variation are still largely unknown, especially in rice. By performing GWAS on the concentrations of 17 elements in different tissues of 295 indica and 156 japonica rice accessions grown under different locations, years, and field conditions, we identified 72 locus-element associations with high reproducibility across plant tissues, growth stages, and field conditions (Supplemental Data Sets 3 to 6).
GWAS of the rice ionome requires a number of limiting factors to be overcome, including the relatively small variation in plant elemental concentrations due to the homeostasis of essential mineral nutrients and the usually large environmental variation. Consistent with previous reports on Arabidopsis  and rice (Pinson et al., 2015), we found that the levels of trace elements, especially the nonessential toxic elements, showed greater variations than macronutrient levels in rice tissues (Figure 1, Supplemental Data Set 2), suggesting that the latter are more tightly regulated. It is possible that plants have evolved stronger mechanisms to control the internal fluctuation of essential nutrients, especially macronutrients (e.g., N, P, and K), than for nonessential elements, because the concentrations of macronutrients need to be maintained within relatively narrow ranges that are optimal for growth and development. For nonessential elements, there would be little selection pressure to control their concentrations until they reach the levels of toxicity. Despite the different ranges of variation, 13 of the 17 elements measured in our study were found to have a relatively high heritability within the same location (average H 2 $ 0.5) ( Supplemental Tables 1 and 4), indicating a good genetic basis for their variations, which is in agreement with previous studies on both Arabidopsis and rice Pinson et al., 2015). By contrast, elements including B, Fe, Cr, and Pb showed low heritability; the latter three elements are known to be prone to soil or dust contamination,  Data are shown as means of three biological replicates 6 SD; the corresponding tissues of five plants were mixed to make one replication. Two asterisks indicate significant differences (P < 0.01) compared with the wild type (t test). DW, dry weight.
inflating their concentrations in plant tissues. However, the heritability for most elements was lower when plants were grown at different locations (Supplemental Table 4), which is consistent with previous studies showing that the genetic variation of elements was relatively strong across years at the same field location but relatively weak between field locations (Norton et al., 2012b). To overcome the limitations associated with genetic versus environmental variations, we conducted five field trials across different locations and years. This approach allowed us to identify loci that are reproducible across environments (Norton et al., 2012b(Norton et al., , 2014. We obtained 32 loci (out of 72) that were common across field locations (Supplemental Data Set 5). Of the 72 loci obtained based on repeatability across different field conditions, 33 loci also were confirmed by BLUP-based GWAS. Moreover, 38 loci are close to the QTLs associated with specific elemental concentrations in rice reported in previous studies using biparental crosses or association mapping (Supplemental Data Set 6). These loci represent genetic variations that are stable across different environments and would be more useful for breeding purposes.
The high density of SNPs (;17 SNPs per kb on average) in our GWAS panel also facilitates high-resolution mapping. The resolution (within ;100 kb) (Si et al., 2016) is much higher than the loci generally defined using the interval QTL mapping approach (Supplemental Data Set 6). Due to the slow LD decay in rice and the relatively complex genetic architecture of elemental traits , a locus in rice identified by GWAS typically covers more than 20 genes (International Rice Genome Sequencing Project, 2005). Using information about functional gene annotation or their orthologous genes in other plant species, plausible candidate genes for a number of loci identified by GWAS could be established (Supplemental Data Set 9). We further provide strong evidence that Os-HKT1;5 and Os-MOT1;1 are the causal genes for the loci with the lead SNP sf0111478828 and sf0800123053, respectively, which cause variations in Na and Mo accumulation (Figures 4 to 7). Os-HKT1;5 is a known Na transporter involved in salt tolerance in rice (Ren et al., 2005;Cotsaftis et al., 2012;Negrão et al., 2013;Platten et al., 2013). However, the relationship between Os-HKT1;5 haplotypes and Na accumulation remains unclear. For example, Cotsaftis et al. (2012) suggested that L395V is an essential substitution affecting pore rigidity based on the 3D model of OsHKT1;5 protein, whereas Negrão et al. (2013) proposed that P140A and R184H, but not L395V, are the causal alterations according to their associations with salt stress tolerance. Negrão et al. (2013) genotyped 392 rice accessions by EcoTILLING and presented 15 haplotypes of Os-HKT1;5, and then phenotypically characterized the 59 most highly representative accessions for association analysis. However, due to significant differences in Na accumulation between subpopulations (even between subsubpopulations, e.g., indica I and indica II; Supplemental Figure 8), comparing the ability of different haplotypes (or SNPs) across subpopulations may not be appropriate. The strength of our study lies in the large number of rice accessions within each subspecies, so that the effects of population structure can be excluded, as well as the reproducibility of the haplotype differences across different sites, years, and conditions. Based on an analysis of data from five field trials in each subpopulation, our results show that neither L395V nor P140A had a consistent effect on Na concentration ( Figure 5), suggesting that they probably have no or a weak function in Na variation. By contrast, we found that H184R is a key substitution separating strong and weak alleles of Os-HKT1;5 and that the weak alleles are present in almost all japonica and most indica accessions (Figures  4 and 5).
The second GWAS locus for which the causal gene was identified in our study is Os-MOT1;1, which controls Mo accumulation. Unlike Os-HKT1;5, which affects Na accumulation via coding region variation, we found that the variation in Mo accumulation in rice is caused by variation in the expression levels of Os-MOT1;1 (Figures 6 and 7). This finding is similar to that for MOT1;1 in Arabidopsis . Further efforts will be needed to identify the causal SNPs/INDELs resulting in expression level variation in the promoter region of Os-MOT1;1.
In addition to Os-HKT1;5 and Os-MOT1;1, we identified candidate genes for a number of loci controlling the concentrations of each of the 17 elements measured in rice (Supplemental Data Set 9). For example, OsNRAMP5, encoding a major Mn transporter in rice (Sasaki et al., 2012;Yang et al., 2014), was scanned for Mn concentrations at both the Wuhan and Youxian locations (Supplemental Figure 15). Similarly, OsNRAMP5, which was identified recently by QTL mapping, contributes to Mn variation among rice accessions via variation at the gene expression level . Further investigations are needed to verify these candidate genes. Various approaches, including expressing genes in yeast mutants, knocking out and overexpressing candidates, or developing platforms to examine gene expression variations by RNA sequencing, could be used to select and validate candidate genes in future work. It is also noteworthy that, for most elements, repeated loci responsible for the same element are closely located on the same chromosome (Figure 3, Supplemental Data Sets 3 to 5). These close loci should be considered together to aid the validation of candidate genes.
Indica and japonica are two major subspecies of cultivated rice that exhibit large genetic differentiation, resulting in a variety of phenotypic differences . In this study, we found that this large diversity also was reflected in the ionomic phenotypes, with japonica accessions generally accumulating higher levels of most elements than indica accessions (opposite pattern for Cd; Supplemental Tables 2 and 5). Interestingly, we found that the GWAS loci for the same mineral element were mostly not colocalized in the indica and japonica subpopulations (Supplemental Data Sets 3 to 5), suggesting that different genetic components underlie the variations in ionomic traits between these two subspecies. This hypothesis is supported by the characteristics of the genes that control the variations of the rice ionome in our study. For example, the strong haplotype (or Hap2) of Os-HKT1;5 was quite common in indica but was extremely rare in the japonica subpopulation. Os-MOT1;1 was scanned specifically in indica GWAS panels as a possible reason for the much lower Mo concentrations in indica accessions (Supplemental Figure 4, Supplemental Tables 2 and 5). In both cases, indica subspecies harbor greater allelic and/or functional diversity than japonica subspecies. These subspecies differences suggest a great potential to breed across subpopulations for better nutrient uptake or salt tolerance.
For a number of agronomic traits, including seed size, shattering habit, and flowering time, as well as metabolic traits, the locations of QTLs have been shown to be conserved across rice, sorghum (Sorghum bicolor), and maize (Zea mays) (Paterson et al., 1995;Chen et al., 2016), suggesting a convergent domestication underlying these phenotypes across cereals. Common causal or candidate genes also were found to control the variations in Na (Os-HKT1;5/At-HKT1), Mo (Os-MOT1;1/At-MOT1), and Cd (Os-HMA3/At-HMA3) levels in rice and Arabidopsis (Supplemental Data Set 11) (Ueno et al., 2010;Chao et al., 2012), suggesting the conservation of genetic variation at specific loci across monocots and dicots. Moreover, Tm-HKT1;5-A in wheat (Triticum turgidum) (Munns et al., 2012), Zm-MOT1 (an ortholog of MOT1s) in maize (Asaro et al., 2016), and Tc-HMA3 in Thlaspi caerulescens (now Noccaea caerulescens) also are responsible for controlling variations of Na, Mo, and Cd accumulation, respectively (Ueno et al., 2011). However, the functional SNPs/INDELs within each ortholog group are always different. For example, a 53-bp deletion 27 nucleotides upstream from the translation start site of At-MOT1 was suggested to be responsible for controlling the variation of gene expression and the shoot Mo concentration in Arabidopsis , but this deletion was not found in the promoter region of Os-MOT1;1 in rice, suggesting that other polymorphisms control the expressional variation in rice ( Figure 6K). Such a difference suggests that the variations in these loci were obtained independently by different plant species during their evolutionary selection, leading to similar ionomic consequences. It is unclear why selection was acting on the same genes across species, but such conservation provides a useful pathway to combine the discoveries of ionomic loci from different plant species and help mine candidate genes. It also means that the ionomic GWAS loci in rice identified in this study can serve as a useful resource for other crop species.
Ghd7 is an important regulator of heading date, plant height, and grain yield and contributes to the variations of these traits in rice germplasm collections (Xue et al., 2008;Weng et al., 2014;Zhang et al., 2015). In this study, we identified Ghd7 as the causal gene for variations in the N concentrations in both shoots and straws (Figures 8 and 9). Interestingly, Ghd7 has been observed in a QTL related to N-deficiency tolerance in rice . However, due to the small log of the odds value (3.01) and the large confidence interval (25.6 cM) of this QTL, the relationship between Ghd7 and N accumulation was not established in the study of Wei et al. (2012). Here, we found that N concentration in the rice collections was highly correlated with heading date (Supplemental Figure 11), and overexpressing a functional allele of Ghd7 in rice plants significantly decreased N accumulation (Figure 9). Consistent with the observation that N metabolism is coupled with the rate of photosynthetic C fixation, the enhanced expression of Ghd7 also lowered the chlorophyll content by down-regulating the expression of genes involved in the biosynthesis of chlorophyll and other chloroplast components . Interestingly, another important gene for rice yield, GRF4, was shown recently to affect both N and C metabolism by regulating the expression levels of genes involved in the uptake and assimilation of N, photosynthesis, and C assimilation (Li et al., 2018). Ghd7 may utilize a similar molecular mechanism affecting both N accumulation and photosynthesis in rice. DEP1 (DENSE AND ERECT PANICLE1), a key regulator of grain number, also regulates both yield and N-use efficiency in rice (Sun et al., 2014). Mutation of DEP1 resulted in decreased plant height and increased N concentration but did not affect heading date. We suggest that these key genes regulating agronomic traits and N levels form a complex relationship in regulating plant growth and N accumulation. In our collection, we found significant negative correlation between biomass and N concentrations at both the heading and mature stages but a significant positive correlation between biomass and the total amount of N accumulated (Supplemental Figure 14). These opposite relationships suggest that cultivars with higher biomass are able to acquire more N from the soil and that there is a dilution effect of increasing biomass on the N concentration. Therefore, it is important to identify loci that are associated with biomass, grain yield, and N concentrations together (Supplemental Table 6). The availability of these loci would provide more chances to investigate the relationship between N concentrations with biomass as well as grain yield and also would facilitate breeding rice with high N-use efficiency as well as high yield.
Reducing the inputs of N and P fertilizers while maintaining high yield and quality is crucial for agricultural sustainability. Crop breeding is an important tool for realizing this goal. This study has shown that N and P deficiencies can affect the ionomic profile of rice plants, thus affecting grain quality by altering the concentrations of essential and toxic elements (Supplemental Figure 1, Supplemental Data Set 2). Moreover, the extent of the influence of nutritional status on different elements is obviously different. Analyses of the genetic basis of variation in the ionomic profile under different N and P supply conditions have revealed both the conservative and nonconservative genetic components involved in controlling the ionomic profile in response to the altered nutrient status. Moreover, we have identified loci for element concentrations that are highly dependent on the N or P supply (Supplemental Data Set 13). Isolation of the causal genes for these loci would help to unravel the complex crosstalk between N or P and other elements in rice plants.
In conclusion, we have described a comprehensive study of the rice ionome combining genetic methodologies with highthroughput elemental profiling to dissect the genetic basis of variations of 17 mineral elements in rice. Our study revealed at least 72 common locus-element associations showing high reproducibility under different field conditions and a large number of specific loci corresponding to different field conditions. We have identified strong candidates for three of the loci, including strong evidence that Os-MOT1;1 is involved in the variation of Mo accumulation and the effect of the heading date gene Ghd7 on N concentration in rice. The information about locus-element associations obtained in this study provides an important foundation for future studies on the genetic and molecular mechanisms controlling the rice ionome.

Plant Materials and Sequencing Data
A set of 529 rice (Oryza sativa) accessions was collected and sequenced as described previously . The set included 192 accessions from a core/minicore collection of rice in China , 132 parental lines used in the International Rice Molecular Breeding Program (Yu et al., 2003), 148 accessions from a minicore subset of the U.S. Department of Agriculture rice gene bank (Agrama et al., 2009), 15 accessions used for SNP discovery in the OryzaSNP project (McNally et al., 2009), and 46 additional accessions from the Rice Germplasm Center at the International Rice Research Institute. Information about the accessions, including names, countries of origin, geographical locations, and subpopulation classification, is listed in Supplemental Data Set 1. The 529 accessions were sequenced using the Illumina HiSeq 2000 platform in the form of 90-bp paired-end reads to generate high-quality sequences of more than 1 Gb per accession (>2.53 per genome, total 6.7 billion reads). The assembly release version 6.1 of genomic pseudomolecules of japonica cv Nipponbare, downloaded from the rice annotation database of Michigan State University, was used as the reference genome. Imputation was performed using an in-house-modified k nearest neighbor algorithm. The detailed procedures of data analysis were described previously Xie et al., 2015). The SNP information is available on the RiceVarMap website .
The wild-type, mutant, and transgenic plants related to Os-MOT1;1 used in this study were all in the cv Zhonghua 11 (japonica) background. The osmot1;1 mutant is a T-DNA insertion mutant (ID: 03Z11DK23) identified from the Rice Mutant Database (Wu et al., 2003;Zhang et al., 2006).

Field Trials and Sample Preparation
Field trials were conducted at two locations in China: the experimental farm of Huazhong Agricultural University in Wuhan, Hubei province, and an experimental farm in Youxian, Hunan province. At Wuhan, seeds of the 529 rice accessions were germinated in a seed bed in mid-May 2012 and transplanted to three adjacent paddy fields with different nutritional status in mid-June 2012. The three different fields, NF receiving standard fertilizer inputs, LN, and LP, were constructed by controlling fertilizer applications over the past decade. Fertilizers were applied (per hectare) as follows: 90 kg of N, 45 kg of P, and 72 kg of K in the NF; 45 kg of P and 72 kg of K in the LN field; and 90 kg of N and 72 kg of K in the LP field. At Youxian, the field trials were performed for two consecutive years (2014 and 2015). The seeds of the 529 rice accessions were germinated in mid-May, and ;25-d-old seedlings were transplanted to the field in each year. Before rice transplanting, 11 to 13 soil samples were collected randomly from each field for the analysis of soil nutrient (and heavy metals/metalloids/pH values) status (Supplemental Table 3). The irrigation water also was analyzed (Supplemental Table 7). Each accession was planted in two replicates in each field at Wuhan with 30 plants (three rows of 10 plants) per replicate and in three replicates for each year at Youxian with 20 plants (two rows of 10 plants) per replicate in a randomized complete block design. The planting density was 16.5 cm between plants in a row and 26 cm between rows. The trials were managed according to normal agricultural practices with regard to crop protection and paddy water management (considering the varied heading dates in different accessions, we drained the paddy water in early September for ;15 d, then irrigated again for about 1 week, and drained the water in late September). Five plants of each accession were harvested at the heading stage, and another five plants were harvested at the maturity stage (40 d after flowering for each line). Plants inside each row were used for sampling to avoid the edge effect. Five plants were combined to represent one biological replicate. Plant samples at maturity were separated into straw and brown rice (dehusked).

Elemental Analysis
Plant samples were dried at 80°C for 3 d and ground to fine powders. For the determinations of N and P concentrations, 0.2 g of dried powder was digested with 5 mL of 98% H 2 SO 4 (w/w) and 5 mL of 30% H 2 O 2 (w/w). After cooling, the digested sample was diluted to 100 mL with distilled water. The N concentration in the solution was determined colorimetrically at 660 nm using a modified Berthelot reaction with salicylate, dichloroisocyanurate, and complex cyanides on an automated discrete analyzer (SmartChem 200). The P concentration in the solution was measured using the molybdate blue method with absorbance read at 880 nm on an automated discrete analyzer. For the determinations of other mineral elements, 0.2 g of dried powder of each sample was digested in 65% nitric acid in a MARS6 microwave (CEM) with a gradient of temperatures from 120°C to 180°C for 45 min. After dilution in deionized water, the concentrations of the 15 elements were determined by inductively coupled plasma mass spectrometry (Agilent 7700 series).

Genome-Wide Association Analyses
SNPs with minor allele frequency $ 0.05 and the number of rice accessions with the minor allele $ 15 in the population were used to carry out GWAS. In total, 2,767,191 and 1,857,866 SNPs were used in GWAS for subpopulations of indica and japonica, respectively. To control spurious associations, genetic relatedness was modeled as a random effect in LMM using the kinship matrix. The evenly distributed random SNP set used to analyze population structure was used to calculate kinship. GWAS was performed using LMM and simple linear regression provided by the FaST-LMM program. According to a modified Bonferroni correction described by Li et al. (2012), the effective number of independent SNPs in each population was calculated and then used to replace the total number of SNPs to determine the genome-wide significance thresholds of the GWAS. The effective number of independent SNPs was calculated as 571,843 and 245,348 for the subpopulations of indica and japonica, respectively, and the suggestive P-values were specified as 1.8 3 10 26 in indica and 4.1 3 10 26 in japonica. The corresponding thresholds then were set to identify significant association signals by LMM. To obtain independent association signals, multiple SNPs exceeding the threshold in a 5-Mb region were clustered by r 2 of LD $ 0.25, and SNPs with the lowest P-value in a cluster were considered to represent lead SNPs. The detailed method was described previously . Before performing GWAS, the Box-Cox procedure was used to normalize the elemental traits.
For comparative GWAS of rice with Arabidopsis (Arabidopsis thaliana), the publicly available leaf ionome data for ;349 wild Arabidopsis accessions for the elements Na, Ca, Mg, B, P, K, Mn, Fe, Co, Cu, As, Zn, Cd, and Mo were used (Baxter et al., 2010;Chao et al., 2012Chao et al., , 2014Forsberg et al., 2015). The leaf ionome data were run against the fully imputed SNP data set containing 10,707,430 biallelic SNPs using the online webtool GWAPP accessible at http://gwas.gmi.oeaw.ac.at ). An Accelerated Mixed Model approach was applied on the untransformed data, and the significant SNPs with 2log(P) > 5 were filtered and analyzed for associations.

Vector Construction and Rice Transformation
To generate the overexpression constructs of Os-MOT1;1, four types (types 1-4) of full-length cDNA of Os-MOT1;1 were amplified using primers Os-MOT1;1-OX-F/Os-MOT1;1-OX-R from four different rice accessions: W024, W102, C020, and C055. The amplified cDNA was first introduced into the Gateway vector pDONR207 and then transferred into the destination vector pJC034 using the Gateway recombination reaction (Invitrogen). The constructs were transformed into cv Zhonghua 11 by Agrobacterium tumefaciens-mediated transformation (Hiei et al., 1994;Lin and Zhang, 2005). To construct the MOT1;1-promoter:GUS plasmid, 2.5 kb of genomic sequence located upstream of the Os-MOT1;1 initiation codon was amplified using primers Os-MOT1;1-P-F/Os-MOT1;1-P-R from cv Nipponbare genomic DNA. The amplified promoter fragment then was cloned into pDONR207 and transformed into the Gateway-compatible GUS fusion vector pGWB3 by Gateway recombination reaction. Callus of cv Zhonghua 11 was transformed with this construct. The transgenic plant tissues were incubated in 5-bromo-4-chloro-3-indolyl-b-glucuronic acid staining buffer at 37°C for 4 h.
To generate the overexpression constructs of Ghd7, the full-length cDNA of a functional allele of Ghd7 was amplified using primers Ghd7-OX-F/Ghd7-OX-R from cv Minghui 63. The amplified cDNA was first introduced into the pGEM-T vector (Promega) and then cloned into a modified overexpression vector PU1301 GPF with XhoI and BamHI sites. All primers used for vector construction are listed in Supplemental Table 8.

Hydroponic Experiments
For expression analysis of Os-MOT1;1 and the determination of Mo concentrations in roots and shoots of both transgenic and wild-type plants, hydroponic experiments were performed using a standard rice culture solution [1.44 mM NH 4 NO 3 , 0.3 mM NaH 2 PO 4 , 0.5 mM K 2 SO 4 , 1.0 mM CaCl 2 , 1.6 mM MgSO 4 , 0.17 mM Na 2 SiO 3 , 50 mM Fe-EDTA, 0.06 mM (NH 4 ) 6 Mo 7 O 24 , 15 mM H 3 BO 3 , 8 mM MnCl 2 , 0.12 mM CuSO 4 , 0.12 mM ZnSO 4 , 29 mM FeCl 3 , and 40.5 mM citric acid with pH adjusted to 5.5 by sulfuric acid] (Yoshida et al., 1976). Rice plants were grown in full nutrient solution or Mo deficiency solution (full nutrient solution without Mo) for approximately 4 weeks after germination. The nutrient solution was renewed every 5 d. The plants were grown in a greenhouse under natural light and temperature conditions (June/July, Wuhan). For both expressional and elemental analysis, three plants (or tissues from three different plants) were mixed in one replication, and three biological replicates were conducted for each line.

RNA Extraction and qRT-PCR
Total RNA was extracted using Trizol reagent (Invitrogen). The first-strand cDNAs were synthesized with 3 mg of total RNA in 20 mL of reaction mixture using SuperScript III reverse transcriptase (Invitrogen) according to the manufacturer's instructions. qRT-PCR was performed using SYBR Premix Ex Taq (TaKaRa) on an ABI 7500 PCR instrument (Applied Biosystems). The PCR device was programmed as follows: 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 40 s. The rice Ubiquitin5 gene was used as the internal control with primers qUbq-F/qUbq-R. The expression measurements were made using the relative quantification method (Livak and Schmittgen, 2001). All primers used for qRT-PCR are listed in Supplemental Table 8.

Statistical Analyses
H 2 was calculated using the following equation by treating accessions as a random effect and the biological replication as a replication effect using one-way ANOVA as described previously : where var (G) and var (E) are the variances derived from genetic and environmental effects, respectively. The BLUPs were implemented in the lmer function of the R package lme4 (Bates et al., 2014), using variety and field conditions (for Wuhan) or variety and year (for Youxian) as random effects. BLUPs were performed both on the phenotypic means of two or three biological replicates within the randomized complete block design after Box-Cox transformation and on the values of individual biological replicates using the replicate as a random effect to account for field-scale heterogeneity. Because both approaches produced similar results, only the results based on the phenotypic means are presented. LD was estimated using standardized disequilibrium coefficients, and squared allelefrequency correlations (r 2 ) for pairs of SNP loci were determined according to the TASSEL program (Bradbury et al., 2007). LD plots were generated in Haploview (Barrett, 2009), indicating r 2 values between pairs of SNPs multiplied by 100 and color coded with white (r 2 = 0), shades of gray (0 < r 2 < 1), and black (r 2 = 1).

Accession Numbers
Sequence data used in this article can be found in the GenBank/EMBL database under the following accession numbers: Os-HKT1;5 (Os01g0307500), Os-MOT1;1 (Os08g0101500), Ghd7 (Os07g0261200), and Os-NRAMP5 (Os07g0257200). The raw sequence data can be obtained in the National Center for Biotechnology Information under BioProject accession number PRJNA171289. The SNPs and imputed genotypes can be queried at RiceVarMap (http://ricevarmap.ncpgr.cn).