The rice population was from a global-wide rice diversity panel consisting of 413 inbred accessions (Zhao et al., 2011). In all, 276 accessions with suitable heading date time and good germination performance were selected for field experiments (Supplementary Table S1). These accessions were previously genotyped using three types of SNP arrays (1536, 44 K, and 700 K SNP chips) (Zhao et al., 2010, 2011; McCouch et al., 2016). We successfully extracted 415 754 high-quality SNPs (missing rate, <0.2; minor allele frequency, >0.05) from these array data for the 276 accessions. About 45% of the SNPs were located within genes, half of which were distributed within exons; the remaining 55% of SNPs were located in intergenic regions, half of which were distributed in putative regulatory regions (2 kb before a transcriptional start site) (McCouch et al., 2016). We also extracted high-quality SNPs for each investigated subpopulation (190,616, 198,790, 117,880, and 147,039 SNPs, respectively). The SNP data were processed using PLINK version 1.07 (Purcell et al., 2007). Based on these 416 K SNPs, we conducted PC analysis by using GCTA version 1.91 (Yang et al., 2011). Population structure was analyzed using fastSTRUCTURE, and the optimum subpopulation number was determined to be K = 5 (Raj et al., 2014; McCouch et al., 2016).

The rice plants were cultivated in farms, which were reported to be heavily contaminated by mining tailings from nearby abandoned mines (Wang et al., 2008; Wu et al., 2013) in Shangyu, Shaoxing City, Zhejiang Province, China (N30°00′14′′, E120°46′39′′) in 2014 and 2015. The soil levels of As, Cd, and Pb were about 3,000, 4.0, and 15,000 mg/kg, respectively, which are considerably higher than the environmental quality standard for soils in China, as well as soil background levels of THMMs in Zhejiang Province (Wang et al., 2008; Wu et al., 2013). For each year, each accession was arranged in a randomized complete block design into two-row plots with 12 plants per row at spacing of 20 cm × 17 cm with three replications. The farms were managed under flooded condition. Seeds for each accession were bulk-harvested from the eight central plants from each replication at the maturing stage. Before the grain THMM concentrations were determined, the seeds were dried for 7 days and then stored for 3 months.

The de-hulled seeds (brown rice) of each accession for each year were milled into fine flour, followed by filtering (0.180 mm). For each sample, 0.25 g filtered powder was placed into a 50 ml polypropylene digestion vessel containing 8 ml nitric acid, and then mixed in a graphite furnace (Agilent Technologies, CA) at 65°C for 30 min. The mixture was digested at 120°C until 3 ml was left. After cooling, 3 ml hydrogen peroxide was added to each digestion vessel. The mixture was again digested until 1 ml was left. After digestion, the solution was diluted to 50 ml with ultrapure water and mixed to uniformity. The diluted solutions were analyzed for As, Cd, and Pb by using inductively coupled plasma mass spectroscopy (Agilent 7500A; Agilent Technologies, CA). The rice reference material (GBW10045) was used as the standard. Three technological repetitions were performed for each sample, and three biological repetitions were used for each accession. The averages of 2-year repeats were used for data analyses. Correlation tests showed a significant correlation between each pair of 2-year repeats (P = 4.52 × 10^-4 to 2.63 × 10^-53; Supplementary Table S2). Two-way ANOVA was employed to examine the genotype and year effects on grain THMM concentrations (Supplementary Table S3).

We examined the LD decay for the 276 lines by using 416 K SNPs and PopLDdecay version 3.4 (Zhang et al., 2018). The following parameters were used: maximum distance between two SNPs, 500 kb; missing rate, <0.2; minor allele frequency, >0.05. We also calculated LD decays for each subpopulation. The widely accepted linear mixed model was used for genome-wide association study and was implemented in GEMMA version 0.97 (Zhou and Stephens, 2012). The population structure was controlled using a relatedness matrix that was calculated using all high-quality SNPs mentioned above. In addition to association mapping for the entire population, we performed association mapping for the respective subpopulations (not including the admixed subpopulation and the very limited aromatic subpopulation). For association mapping for subpopulations, high-quality SNP datasets were obtained for each subpopulation as mentioned above, and a corresponding relatedness matrix was constructed using the respective SNP dataset. Null hypothesis of no association was tested, and the statistical probability values were calculated using Wald test in GEMMA. The threshold for significance was determined according to a previous study (McCouch et al., 2016). In the case of a strong background in some mapping results (for grain Cd and Pb concentrations in the AUS subpopulation), we adjusted the threshold for significance by using Bonferroni correction for multiple comparisons at a false discovery rate of 0.05 (P-value = 0.05/SNP number).

Considering the balance between sensitivity and reliability, the potential candidate genes were predicted within a 400 kb genomic region for each locus (±200 kb of the peak SNPs). For the loci with multiple significant SNPs, the potential candidate genes were predicted within ±200 kb of the borders of the significant SNPs. Pairwise LD was determined by calculation of r² (the square of the correlation coefficient between SNP states) using PLINK. LD blocks were analyzed using LDheatmap version 0.99 (Shin et al., 2006), and candidate LD blocks were estimated using r² > 0.8. The transcriptomic data of rice under THMM stress were extracted from the supplemental data of the related references and NCBI GEO database (Huang et al., 2012; Yu et al., 2012; Ogo et al., 2014). The expression profile of the candidate genes was analyzed using the Bio-Analytic Resource Plant Biology database¹. The differences in THMM concentration between the haplotypes were statistically analyzed using Welch’s t-test (two haplotypes) or Tukey’s test (more than two haplotypes). All SNP variations across genes were extracted from the whole-genome resequencing data of 3 000 Rice Genomes Project (3 K filtered dataset) (Mansueto et al., 2017; Wang et al., 2018).

Significant differences in the concentration of the three THMMs (As, Cd, and Pb) in grain (brown rice) were noted among the 276 global representative rice accessions (Supplementary Table S1) grown in a multi-contaminated farm with high soil As, Cd, and Pb levels, suggesting vast variations of grain THMM accumulation abilities of rice (Figure 1 and Table 1). The means for grain As, Cd, and Pb concentrations were 0.34 ± 0.11, 0.12 ± 0.06, and 0.78 ± 0.45 mg/kg, respectively (Table 1). The grain As concentration showed minimum variation with a 6.25-fold difference, whereas grain Cd concentration showed maximum variation with a 23.57-fold difference (Table 1). The grain As and Pb concentrations of most accessions were over the maximum permitted level for grain commercialization, but the grain Cd concentration of most accessions were below the permitted level (Figure 1). Flooded condition during cultivation could be responsible for the relatively low accumulation of Cd in grain in our study. The accessions with low grain THMM accumulation are listed in Supplementary Table S4 for rice breeding.

Compared to the normal distribution, the frequency distributions of all three grain THMM concentrations, especially that for Pb, showed unique stronger peaks and longer right tails (Figure 1 and Table 1). The broad-sense heritability (H²) of grain As, Cd, and Pb concentrations were 0.73, 0.35, and 0.21, respectively, suggesting that environmental factors have an important impact on grain THMM accumulation, especially for Cd and Pb (Table 1). Correlation tests showed a significant correlation between each pair of 2-year repeats for grain THMM concentrations (all P ≤ 4.52 × 10^-4, Supplementary Table S2), enabling the examination of the genetic basis for grain THMM concentrations. Two-way ANOVA was then performed to examine both the genotype effect and year effect on grain THMM concentrations. Significant genotype effects were observed (all P ≤ 2.37 × 10^-4, Supplementary Table S3), while no significant year effect was found for all three grain THMM concentrations (all P ≥ 0.41, Supplementary Table S3). These results together suggest that the environmental factors affecting grain THMM concentrations are complex and remain unclear. Correlation analysis performed among the three grain THMM concentrations revealed a significant positive correlation between grain Cd and Pb concentrations (r = 0.63, P = 4.87 × 10^-34) and a negative correlation between grain As and Cd concentrations (r = -0.15, P = 0.02). This suggests that grain Cd accumulation could be impacted by similar factors as those for grain Pb accumulation, but different from those for grain As accumulation. No significant correlation was noted between grain As and Pb concentrations.

Previous PC analysis revealed a clear population structure of the rice diversity panel (Zhao et al., 2011). We extracted 415,754 high-quality SNPs (missing rate, <0.2; minor allele frequency, >0.05) from three SNP arrays for the selected 276 accessions, with a high density of 0.90 kb/marker (Table 2). The PC analysis based on these 416 K SNPs also indicated the same population structure for the selected accessions (Supplementary Figure S1). Next, we performed correlation analysis between grain THMM concentrations and the top four PCs. We found a positive correlation between PC3 and all three grain THMM concentrations and a negative correlation between PC4 and grain As concentration, suggesting that genetic diversity plays a role in grain THMM accumulation abilities (Table 3).

Population structure analysis of these 276 accessions determined the optimum subpopulation number to be five (McCouch et al., 2016). The population consisted of five distinct subpopulations [aromatic (ARO), aus (AUS), indica (IND), temperate japonica (TEJ), and tropical japonica (TRJ)] and an admixed (ADM) subpopulation (Supplementary Table S1; Zhao et al., 2011). The top four PCs clearly separated these subpopulations except the ADM subpopulation (Supplementary Figure S1). Different subpopulations showed variations in grain THMM concentrations (Figure 2). The ARO subpopulation showed the lowest grain As and Pb concentration, and the TRJ subpopulation had the lowest grain Cd concentration; TEJ, ARO, and IND subpopulations had the highest grain As, Cd, and Pb concentrations, respectively (Figure 2). These data indicate a strong relationship of population structure with grain THMM concentrations. Further, we noted remarkable variations in grain THMM concentrations within some subpopulations, for example, AUS, IND, and TEJ subpopulations showed large variations of grain As, Cd, and Pb concentrations, respectively (Figure 2), suggesting the existence of major QTLs within the rice subpopulations.

The LD decay analysis yielded similar results as those reported previously (Supplementary Figure S2) (McCouch et al., 2016). We used linear mixed model for genome-wide association mapping, which has been widely used and validated to be efficient (Yang et al., 2014). The impact of population structure in association mapping was determined by estimating the genetic relatedness. The association mapping performed using 416 K high-quality SNPs identified a total of 22 QTLs (21 non-redundant loci) reaching the thresholds (P < 1.00 × 10^-5; Figure 3–6 and Tables 4–6). The detection power was the highest for grain As concentration, and 10 QTLs were detected in the rice population (Figure 3 and Table 4). Seven QTLs were found to be responsible for grain Cd concentration, and the remaining five were responsible for grain Pb concentration (Figure 4, 5 and Tables 5, 6). A QTL in chromosome 5 was detected for both grain As and Pb concentrations (qGAS7 and qGPB3; Figure 3, 5, 6 and Tables 4, 6).

Although the impact of population structure has been considered during association mapping, we still performed association mapping for each of the subpopulations considering a strong relationship between population structure and grain THMM concentrations. The association mapping was not possible for the ARO subpopulation because of very limited accessions and was not performed for the ADM subpopulation. For the AUS, IND, TEJ, and TRJ subpopulations, 190,616, 198,790, 117,880, and 147,039 high-quality SNPs (missing rate, <0.2; minor allele frequency, >0.05) were successfully extracted, respectively. We detected four QTLs in each of the AUS, IND, and TRJ subpopulations for grain As concentration (totally 12 QTLs; Figure 3 and Table 4). For grain Cd concentration, four, four, and two QTLs were detected in the IND, TEJ, and TRJ subpopulations, respectively (totally 10 QTLs; Figure 4 and Table 5). For grain Pb concentration, four, two, two, and eight QTLs were identified in the AUS, IND, TEJ, and TRJ subpopulations, respectively (totally 16 QTLs; Figure 5 and Table 6). Consequently, 8, 10, 6, and 14 QTLs were discovered in the AUS, IND, TEJ, and TRJ subpopulations, respectively. Three QTLs (qGAS16, qGAS17, and qGPB6) were co-localized with the QTLs detected in the entire rice population (Figure 6). Additionally, two QTLs for grain As concentration, one in the AUS subpopulation (qGAS18) and the other in the TRJ subpopulation (qGAS22), were co-localized on the chromosome 11 (Figure 6).

We determined whether the identified QTLs harbored known THMM-related genes. qGAS8 and qGAS17 contain OsPIP2;7, encoding an aquaporin involved in arsenite permeability and transport (Mosa et al., 2012). qGAS12 contains OsARM1 and OsPT23. The former encodes a MYB transcription factor responsible for As root-to-shoot translocation (Wang et al., 2017), and the latter encodes a phosphate transporter involved in the uptake of arsenate (Yu et al., 2012; Clemens and Ma, 2016). Next, we predicted the possible candidate genes in the regions of QTLs (±200 kb of the significant SNPs; Supplementary Table S5). Genes encoding proteins with transport functions are known to play important roles in THMM accumulation in plants (Clemens and Ma, 2016). We found that 61.8% of the identified QTLs contained genes encoding transporter-related proteins (Supplementary Table S5). Notably, high concentrations of THMMs in the soil used in this study could induce toxic effect and stress for rice plants. Interestingly, OsPIP2;7 and OsARM1 harbored in the identified QTLs (qGAS8, qGAS12, and qGAS17) for grain As accumulation were reported to be involved in As tolerance and detoxification (Mosa et al., 2012; Wang et al., 2017). By examining the THMM stress-responsive genes previously identified using transcriptome analysis in rice (Huang et al., 2012; Yu et al., 2012; Ogo et al., 2014), we found that 12 and six QTLs for grain As and Cd concentrations contain genes responding to As and Cd stress, respectively, including OsARM1 (Supplementary Table S6).

Subsequently, we focused on qGAS1, which had the best P-value among the QTLs detected in the entire population. Analysis of LD blocks revealed two blocks in the region of qGAS1 (Figure 7). Almost all significant SNPs were located in the second block (Figure 7). We found four transporter-related genes in this block. According to protein homology, LOC_Os01g55500 encodes a nucleobase–ascorbate transporter, LOC_Os01g55600 and LOC_Os01g55610 encode nitrate transporters, and LOC_Os01g56050 encodes a MATE (also known as multi-antimicrobial extrusion) transport protein. The expression profiles of the four genes showed that they were mainly expressed in seeds/inflorescences, except LOC_Os01g55500 (which was mainly expressed in the shoot apical meristem; Figure 8). Next, we performed haplotype analysis for these genes (Figure 9A). LOC_Os01g55500 and LOC_Os01g56050 showed different grain As concentrations between some of the haplotypes (Figure 9B). Taken together, our findings suggest that LOC_Os01g56050 might be a possible candidate gene for qGAS1.

In our SNP dataset, only three SNP markers were found on the promoter (within 1 kb) of LOC_Os01g56050, which divided the gene alleles into three haplotypes (Figure 9A). Haplotype A with the highest grain As concentration was primarily distributed in the TEJ subpopulation, whereas haplotype C with the lowest grain As concentration was primarily distributed in the TRJ subpopulation (Figure 9B,C). We then extracted all the SNPs on LOC_Os01g56050 from the whole-genome resequencing data of rice population (Wang et al., 2018). The eight SNPs found on this gene also divided the gene alleles into three major haplotypes (Supplementary Table S7). Two of the SNPs caused non-synonymous changes between the haplotypes (valine²¹³ to isoleucine²¹³ and lysine⁴⁵⁹ to isoleucine⁴⁵⁹; Supplementary Table S7). The other SNPs were located in the introns or caused synonymous changes (Supplementary Table S7). Both valine and isoleucine are hydrophobic amino acids and have similar biochemical characters, whereas lysine and isoleucine have distinct biochemical characters (lysine is a basic amino acid). Furthermore, lysine⁴⁵⁹ is a relatively conserved residue in the MATE domain (Supplementary Figure S3). Therefore, the difference of lysine⁴⁵⁹ (haplotypes A and B) and isoleucine⁴⁵⁹ (haplotypes C) could have an impact on the protein activity of LOC_Os01g56050 alleles.