Multi‐Experiment and Multi‐Locus Genome‐Wide Association Mapping for Grain Arsenic in Rice Population

Abstract

ABSTRACT

Rice is a globally important crop and is particularly efficient at assimilating arsenic (As). Identifying QTLs and genes associated with grain As is essential for breeding low‐As rice cultivars. In this study, data on As accumulation in grains of Rice Diversity Panel 1 in five field environments at four diverse geographic sites were reanalyzed to compare genome‐wide association (GWA) methods. Two single‐locus (EMMAX for single trait and GEMMA for multi‐experiments) and six multi‐locus (FASTmrEMMA, ISIS EM‐BLASSO, mrMLM, pKWmEB, pLARmEB, and FASTmrMLM) GWA methods were used. A total of 90 and 111 QTLs were detected using EMMAX and GEMMA, respectively. A total of 2, 11, 12, 19, 23, and 25 QTNs were identified by FASTmrEMMA, ISIS EM‐BLASSO, mrMLM, pKWmEB, pLARmEB, and FASTmrMLM, respectively. Among these, 22 QTLs/QTNs were co‐detected by single‐locus and multi‐locus GWAS methods. From these QTLs/QTNs, a total of 10 candidate genes were identified. Analysis of the haplotype variants of one candidate genes, OsABCC1, and one cluster of the plasma membrane intrinsic proteins genes revealed that a greater than 10% reduction in grain As could be achieved. The QTLs/QTNs and candidate genes identified give insight into the molecular mechanisms regulating As accumulation in rice and serve as breeding targets for developing low grain As rice cultivars.

Introduction

Rice is one of the most important crops globally, providing the calorific needs of millions of people daily (McNally et al. 2009). However, rice also accumulates a higher concentration of arsenic (As) than other staple cereals such as wheat and barley (Williams et al. 2007). Of note, due to expense, most studies measure total As concentration that includes both inorganic As, the more toxic form, and organic As, which is less toxic to humans. The World Health Organization (WHO) has established FAO/CODEX limits only for inorganic As (Codex Alimentarius 2017). Inorganic As includes both the trivalent form arsenite, As(III), which is highly soluble in water, and the more toxic pentavalent form arsenate, As(V), which is less water soluble. Rice is efficient at As accumulation as it is predominantly grown in flooded paddy soils, which leads to As(III) mobilization, and the inadvertent yet sometimes high uptake of As(III) through the highly effective silicon transport pathway found in rice (Etesami et al. 2023; F. Zhao et al. 2010). Arsenic contamination of soil due to irrigation with As-rich groundwater or mining activities in Asia has resulted in further elevation of As levels in rice (Meharg and Rahman 2003; Williams et al. 2006; Zhu et al. 2008). While not all rice contains concerning levels of As, long-term consumption of rice contaminated with high levels of As has been linked to multiple health problems such as cardiovascular disease risk, skin lesions, and squamous cell skin cancers (Karagas et al. 2019; Rahman et al. 2023). Both preharvest and postharvest mitigation approaches have been proposed to decrease As content in cooked rice, but leveraging the natural variation that exists within rice to develop cultivars that limit accumulation of As in rice grains may also prove effective.

Several quantitative trait loci (QTLs) for As accumulation in rice have been identified (Dong et al. 2023; Jian et al. 2008; X. Liu et al. 2019; H. Liu et al. 2021; Murugaiyan et al. 2019; Norton et al. 2014, 2019; Pinson et al. 2022; Yang et al. 2018; J. Zhang et al. 2008). For example, a QTL on Chromosome 2 for As concentration in the shoots and two QTLs on Chromosomes 6 and 8 for As concentration in the grain, identified in F1 plants from a cross between a japonica cultivar, CJ06, and an indica cultivar, TN1 (J. Zhang et al. 2008). In addition, two QTLs for As concentration in roots have been mapped on Chromosome 8, and six QTLs mapped to Chromosomes 2, 5, 6, and 9, from a backcross breeding population (Murugaiyan et al. 2019). A total of 22 QTLs have been identified for grain As based on genome-wide association (GWA) mapping of 276 rice accessions from the Rice Diversity Panel 1 (RDP1), a collection of rice cultivars from 10 geographic regions where rice is grown (X. Liu et al. 2019) and 74 QTLs for grain arsenic identified by GWA mapping using the Bengal and Assam Aus Panel (BAAP) (Norton et al. 2019). Moreover, genes involved in the uptake and transport of As in rice have been reported, such as Lsi1 and Lsi2, which are transporters belonging to the nodulin 26-like intrinsic proteins (NIPs) subfamily of aquaporins. OsLsi1 (OsNIP2;1) is a passive aquaporin channel with permeability to several substrates including silicon and As, and mutation of Lsi1 significantly decreased As uptake (Jian et al. 2008; Mitani et al. 2008; Mitani-Ueno et al. 2011). OsLsi2 is an active efflux transporter for silicon and As, where mutation of Lsi2 significantly decreased As transport to the xylem and accumulation in shoots and grains (Jian et al. 2008). Additionally, several NIPs (OsNIP1;1, OsNIP2;2 (OsLsi6), OsNIP3;1, OsNIP3;2, and OsNIP3;3) have been reported to have an impact on As(III) accumulation (Jian et al. 2008; Sun et al. 2018). Members of the plasma membrane intrinsic proteins (PIP) subfamily, OsPIP2;4, OsPIP2;6, and OsPIP2;7, have also been implicated in As(III) tolerance and transport (Mosa et al. 2012). Phosphate transporters OsPT1, OsPT4, and OsPT8 are also capable of As(V) uptake (Kamiya et al. 2013; Ye et al. 2017; P. Wang, Zhang, et al. 2016). After being taken up by rice roots, As(V) can be readily reduced to As(III) by OsHAC1;1, OsHAC1;2, and OsHAC4 (Shi et al. 2016; J. Xu et al. 2017). The rice ABC transporter OsABCC1 has been reported to sequester As(III) into the vacuoles of phloem companion cells in stem nodes, and thus reduce As accumulation in rice grain (Song et al. 2014). It is important to appreciate that post uptake compartmentalization of arsenic is critically important in determining grain arsenic content (Pinson et al. 2022).

GWA mapping identifies associations between nucleotide polymorphisms and phenotypic variance using a diverse population. Based on these genomic associations, candidate genes associated with phenotypic traits have been identified in many crops, such as maize (M. Wang et al. 2012), rice (Huang et al. 2010, 2012), and soybean (Z. Liu et al. 2020). GWA mapping for phenotypic traits in plants typically identifies associations between a single trait and single locus (univariate) based on an individual experiment (Chen, Norton, and Price, 2020, Chen et al., 2021; Norton et al. 2019), and numerous statistical models have been developed, such as the Efficient Mixed-Model Association eXpedited (EMMAX) (Kang et al. 2010) and Genome-Wide Efficient Mixed-Model Association (GEMMA) (Zhou and Stephens 2014). However, a limitation of this approach is inconsistent QTL identification due to high environmental variation between experiments (e.g. temperature, humidity, and soil conditions) and the complex genetic architecture of the population, for example, rare variants in genes of large effect and/or common variation in genes of small effect (Korte and Farlow 2013). Multivariate linear mixed models (mvLMMs), which relate explanatory variables to multiple correlated outcome variables, have become increasingly important in GWA studies due to their effectiveness in accounting for relatedness among samples and population stratification and a growing appreciation of the power gains of mvLMMs over standard univariate analysis (Zhou and Stephens 2014). However, single-locus methods may fail to capture the true genetic model of complex traits, particularly due to the significant experimental errors inherent in field studies of crop genetics (Li et al. 2022; Y.M. Zhang et al. 2019). Additionally, inflation of the Type I error rate (false positives) is commonly observed with single-locus approaches and needs to be accounted for post model fitting (Kaler et al. 2020). Multi-locus GWA mapping methods have been recommended to overcome these problems, including the random-SNP-effect mixed linear model (mrMLM) (S. Wang, Feng, et al., 2016), fast MLM (FASTmrMLM) (Tamba and Zhang 2018), fast multi-locus random-SNP-effect efficient mixed model analysis (FASTmrEMMA) (Wen et al. 2018), polygenic-background-control–based least angle regression plus empirical Bayes (pLARmEB) (J. Zhang et al. 2017), integration of Kruskal–Wallis test with empirical Bayes (pKWmEB) (Ren et al. 2018), and iterative modified-sure independence screening expectation–maximization-Bayesian least absolute shrinkage and selection operator (ISIS EM-BLASSO) (Tamba et al. 2017).

The RDP1 is a rice panel representing a broad range of rice cultivars from over 70 countries (Eizenga et al. 2014). D. Wang, Agosto-Pérez, et al. (2018) imputed 5,231,435 SNPs for RDP1 by comparing the 700,000 SNPs (McCouch et al. 2016) with whole-genome sequence data of the 3000 sequenced rice cultivars (W. Wang, Mauleon, et al. 2018). Previously GWA mapping was conducted on the grain As phenotype characterized on approximately 370 rice cultivars from RDP1 grown in three countries with five different field locations/years (Norton et al. 2014); however, this QTL analysis used a single GWA method and a limited number of SNPs (44K). The objective of this study was to identify the QTLs and candidate genes associated with variance for grain As concentration within the RDP1 population by different GWA mapping methods using 5.2 million SNPs.

Methods and Materials

Material

A total of 253 rice accessions from the RDP1 (Eizenga et al. 2014) were used in this study (Table S1), with this subset from the larger RDP1 being defined as those RDP1 accessions for which there were data from all five field experiments (seeds of these accessions are publicly available from the Genetic Stocks-Oryza collection at ). This diverse population comprises genotypes from the five rice subpopulations: indica, aus, aromatic, tropical japonica, and temperate japonica. The publicly available 5.2 million–SNP database () (D. Wang, Agosto-Pérez, et al. 2018) of RDP1 was generated by imputing from the set of the intersection of 700k SNPs from the High Density Rice Array (HDRA; McCouch et al. 2016) and 18M SNPs from the resequencing data of the 3000 Rice Genome Project (W. Wang, Mauleon, et al. 2018), resulting in 4.8M high-quality HDRA-imputed SNPs (D. Wang, Agosto-Pérez, et al. 2018). SNPs with minor allele frequency (MAF) < 0.05 were filtered out and maximum per-SNP missing was set at 5% and a total of 3.5 million SNPs were kept after filtering.

Field Screen and Arsenic Analysis

The field screening for grain As was conducted in four locations: Bangladesh, China, and two US field locations (Arkansas and Texas). Soil chemistry and field management details for the production and harvest of the RDP1 grains analyzed were described in Norton et al. (2014) where single-locus GWA mapping using 44k SNP markers was conducted. Below is a summary of the field experiments. The field site in Bangladesh was located at Faridpur in 2009; the field site in China was located at Qiyang, Hunan province, in 2009. At both sites, the plants were transplanted in a randomized complete block design with four replicates and were maintained under continuous flooding until most of the cultivars had flowered, at which time the field was drained and allowed to dry before harvest of the grains. For the field site in Arkansas, the cultivars were grown in year 2006 and 2007; the field layout in both years was a randomized complete block design with two replicates. The flood was applied when the plants were at the five-leaf stage and drained 15–20 days after the majority of the cultivars had flowered to allow fields to dry before harvest. For the field site in Texas, three replicates of the RDP1 accessions were grown in 2009 using a randomized complete block design. The fields were flush-irrigated until plants reached an average height of 18 cm. At this stage, a 10-cm permanent flood was established and maintained until the seeds of all genotypes matured and were harvested manually. Total As concentration in grain samples were determined using inductively coupled plasma mass spectrometry (ICP-MS) following total sample digestion in acid, as described by Norton et al. (2014). Each unique combination of location and year was treated as an independent experiment, resulting in a total of five experiments conducted across four locations.

Single-Experiment GWA Mapping

GWA mapping was performed using the PIQUE pipeline, which utilizes EMMAX to conduct the analyses on each phenotype (Norton et al. 2018). The mean grain As value from the original data in each experiment was log₁₀ transformed to ensure a normal distribution and transformed values were used for subsequent GWA analysis. GWA mapping was conducted using 3.5M filtered SNPs and an efficient mixed model (EMMA) to control for population structure (PCA based on the first four principal components) and kinship (using IBS). A significance threshold of p < 0.0001 was used to determine significant SNPs as described in Norton et al. (2018). The false discovery rate (FDR) of detected associations was controlled for using Benjamini–Hochberg adjusted probabilities (Benjamini and Hochberg 1995). A significance threshold of 5% FDR was used to identify SNP associations (McCouch et al. 2016).

Multi-Experiment GWA Mapping

A multi-experiment analysis of GWA mapping for grain As concentration across the five field experiments in RDP1 was conducted using the mvLMM in GEMMA version 0.98 (Zhou and Stephens 2014). The log₁₀ mean grain As concentration for each field experiment was converted to a z-score by mean centering and dividing by the standard deviation. Trait data were analyzed for linkage against the 3.5M filtered SNPs. The mvLMM accounts for both population structure and relatedness among samples to control for potential confounding factors. The eigenvectors of the first four principal components were calculated using the smartpca program in EIGENSOFT (Patterson et al. 2006) and included in the model as fixed effects. The null hypothesis is that the effect of a single SNP across all experiments is zero, whereas the alternative hypothesis is a nonzero effect of at least one SNP tested using a Wald test (Ruang-areerate et al. 2020). p values of all association tests are presented as Manhattan plots and observed p values against expected p values are presented by Q-Q plots using the qqman package (Turner 2014). A significance threshold of p < 0.0001 was used to identify significant SNPs as described in Norton et al. (2018). SNPs were also tested to a 5% FDR based on the Benjamini–Hochberg procedure (Benjamini and Hochberg 1995).

Multi-Locus GWA Mapping

Arsenic concentrations from each experiment (log₁₀ transformed) were standardized using z-score by mean centering and dividing by the standard deviation. The mean values of five standardized As traits were used to calculate a single mean phenotype value for multi-locus GWA. To make sure the analysis was balanced, only the 253 genotypes that had As data in all five locations were used. The first four principal components were included in the GWA and the critical logarithm of odds (LOD) threshold was set at ≥ 3. Six multi-locus GWAS methods (mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, pKWmEB, and ISIS EM-BLASSO) were performed using the MrMLM 4.0 package (Y.W. Zhang et al. 2020) in R Version 4.1.2s. Significant quantitative trait nucleotides (QTNs) for grain As were reported in the multi-locus GWA mapping.

Clump Analysis of Significant SNPs

Following analysis using EMMAX and GEMMA, results were grouped across analyses based on empirical estimates of linkage disequilibrium (LD) between SNPs using CLUMP via the PLINK command “--clump-p1 0.0001 --clump-p2 0.0001 --clump-r2 0.5 --clump-kb 500” as described in Ruang-areerate et al. (2020). The --clump-kb option sets the physical distance for clumping and was set to the average LD-decay (500 kb) of all accessions in the RDP1 (K. Zhao et al. 2011). Clump QTLs with more than one significant SNP were reported in the single-locus GWA mapping. When the positions of peak SNPs for QTLs across different field sites or GWAS methods overlapped or were within 500 kb of each other, these QTLs were regarded as the same QTLs identified across multiple field sites or GWAS methods.

Haplotype Analysis

The notable QTLs that were associated with grain As concentration in multiple locations and/or detected from different GWA methods were further investigated. The SNPs within 500 kb of peak SNPs of QTLs were extracted via plink, and LD heatmaps of these SNPs were constructed using the R package “LDheatmap.” The SNPs that were significantly associated (p < 0.01) with grain As concentration were extracted using PLINK, and nonsynonymous and synonymous SNPs within the genes were identified based on the gene models of the annotated Nipponbare reference genome from RGAP using the ANNOVAR software (K. Wang et al. 2010). The haplotypes based on visualization of the SNPs in the exon of the candidate genes in these notable QTLs were determined, and the phenotypic response for the cultivars with each haplotype was analyzed. The pairwise Wilcox test was used to determine the significant difference between Haplotypes I and II of candidate gene OsPIP2;4. Tukey's honestly significant difference (HSD) test was used to assess the significant difference among haplotypes of candidate gene OsABCC1.

Statistical Analysis

Boxplots and Pearson correlation analyses of grain As concentrations in different experiments were performed using the R packages “ggplot2” (Wickham 2016) and “corrplot” (Wei et al. 2021). Plots showing locations of QTLs and candidate genes across the 12 rice chromosomes were generated using the R packages “plotrix” (Lemon 2006) and “shape” (Wickham 2007).

Results

Grain Arsenic Concentrations Under Different Field Experiments

The phenotypic data used in this study were collected and analyzed by Norton et al. (2014). The grain As concentrations of the RDP1 genotypes were significantly different across the four field sites. Additionally, at each site, the grain As concentrations differed significantly between the rice subpopulations (Norton et al. 2014). Grain As was also significantly correlated between the two Arkansas harvests (2006 and 2007), with r = 0.407, p < 0.001 (Norton et al. 2014). Here, we expand that analysis to include correlations across all five years and location experiments. Significant positive correlations between grain As from the three US field experiments at Arkansas in 2006 and 2007 and at Texas in 2009 were observed (p < 0.001, r = 0.42–0.51). The grain As from the experiment at Faridpur in Bangladesh correlated significantly with grain As in Arkansas in 2006 (p < 0.01, r = 0.19) and grain As in Texas (p < 0.001, r = 0.28). The grain As from the experiment at Qiyang in China significantly correlated with grain As in Arkansas in 2006 (p < 0.001, r = 0.22) and 2007 (p < 0.001, r = 0.31) (Figure 1).

[IMAGE OMITTED. SEE PDF]

Grain As QTLs Detected by Single-Experiment GWA Mapping

GWA mapping was conducted separately in five field experiments using EMMAX with 3.5M filtered SNPs (Figure 2). A total of 14, 31, 21, 10, and 14 QTLs were significantly associated with grain As concentration (−log₁₀ p > 4) in Arkansas (2006), Arkansas (2007), Texas, Faridpur, and Qiyang, respectively (Table S2 and Figure 2). Of these, 48 QTLs exhibited an effect size of 10% or greater (Table S2). The most statistically significant QTLs (−log₁₀ p > 5.95, p = p value of the peak SNP per QTL) were identified on Chromosome 5 at ~23.95 Mb (−log₁₀ p = 5.96) with a total of 534 significant SNPs, on Chromosome 5 at ~24.39 Mb (−log₁₀ p = 6.73) with a total of 135 significant SNPs, and on Chromosome 7 at ~15.56 Mb (−log₁₀ p = 5.99) with a total of 173 significant SNPs. As the effect sizes for the QTLs on Chromosome 5 at ~23.95 and ~24.39 Mb were opposite (17.1% and −15.6%, respectively) and the MAFs were different for each marker (11.9% and 24.7%, respectively), these QTLs were treated separately. The three QTLs with the highest -log₁₀ p were all identified from analysis of Arkansas 2007 data alone. A total of seven QTLs were detected across more than one field experiment (Table S2). A QTL was determined to be in the same location when the positions of peak SNPs for QTLs across different field sites overlapped or were within 500 kb of each other. One collocated QTL was detected on each of Chromosomes 2, 3, and 5, and two each on Chromosomes 1 and 7. It is notable that the QTL identified on Chromosome 1 ranging from 33.94 to 34.32 Mb was detected in three field experiments: Arkansas in year 2006, Arkansas in year 2007, and Texas in year 2009.

[IMAGE OMITTED. SEE PDF]

Grain As QTLs Detected by Multi-Experiment GWA Mapping

GWA mapping for grain As concentration (z-score after transformation) was performed using GEMMA on the 253 accessions that were common across all five field experiments (Figure 3). A total of 111 QTLs were significantly associated with grain As concentration (−log₁₀ p > 4), with at least one SNP in each QTL passing the 5% FDR threshold (Table S3, Figure 3). Six of these QTLs were notable based on statistical significance (−log₁₀ p > 8, p = p value of the peak SNP per QTL). They were located on Chromosome 5 at ~23.95 Mb with 1370 significant SNPs (−log₁₀ p of peak SNP = 9.17), on Chromosome 5 at ~24.38 Mb with 352 significant SNPs (−log₁₀ p = 8.41), on Chromosome 7 at ~6.73 Mb with two significant SNPs (−log₁₀ p = 8.67), on Chromosome 7 at ~16.58 Mb with 619 significant SNPs (−log₁₀ p = 8.24), on Chromosome 7 at ~18.71 Mb with 197 significant SNPs (−log₁₀ p = 8.16), and on Chromosome 7 at ~21.12 Mb with 55 significant SNPs (−log₁₀ p = 8.13).

[IMAGE OMITTED. SEE PDF]

Of the 111 QTLs detected by the multi-experiment GWA mapping, 40 QTLs were consistent with those identified based on single-experiment GWA mapping (Table S4 and Figures 2 and 3). For example, the QTL on Chromosome 1 ranging from 33.84 to 34.37 Mb was co-detected by three independent field experiments: Arkansas in 2006 and 2007 and Texas in 2009. The QTL on Chromosome 5 at ~23.95 Mb was also identified in the single experiment, Arkansas in year 2007. The QTL on Chromosome 7 at ~15.5 Mb with a total of 510 significant SNPs was co-detected by the single-experiment GWA mapping in Arkansas in the years 2006 and 2007. The QTL on Chromosome 9 at ~1.35 Mb with 1074 significant SNPs was also identified in the single-experiment GWA mapping in Qiyang with two significant SNPs.

A total of 71 QTLs for grain As that were identified using multi-experiment GWA analysis were not detected by single-experiment analysis (Table S4), for example, the QTL on Chromosome 1 at ~11.74 Mb with 181 significant SNPs, the QTL on Chromosome 5 at ~0.89 Mb with a total of 65 significant SNPs, and the QTL on Chromosome 12 at ~14.45 Mb with a total of 151 significant SNPs (Table S4 and Figures 2 and 3).

Grain As QTNs Detected by Multi-Locus GWAS Methods

When the mean value of the five standardized As phenotypic data was used to conduct the multi-locus GWA with the six methods in the mrMLM package in R, a combined total of 72 significant QTNs were identified (Table 1). The FASTmrMLM (23) methods identified the most QTNs (23), followed by pLARmEB (21), pKWmEB (19), mrMLM (12), and ISIS EM-BLASSO (11), while the FASTmrEMMA method identified the lowest number of QTNs (2). The five QTNs notable for having −log₁₀ p > 12 were SNP_1_30876860 (−log₁₀ p = 12.5) explaining 5.53% of phenotypic variation, SNP_1_31421621 (−log₁₀ p = 15.7) explaining 41.8% of phenotypic variation, SNP_3_24636405 (−log₁₀ p = 12.5) explaining 11.56% of phenotypic variation, SNP_7_25050720 (−log₁₀ p = 14.9) explaining 7.86% of phenotypic variation, and SNP_12_13859027 (−log₁₀ p = 12.5) explaining 3.93% of phenotypic variation (Table 1).

TABLE 1 Significant QTNs for grain As concentration detected by multi-locus methods.

SNP	Chr	Marker position (bp)	QTN effect	−log10(P)	r2 (%)a	Methodsb
SNP_1_2089929	1	2,089,929	−0.07	3.90	2.45	5
SNP_1_11876654	1	11,876,654	0.08	6.47	3.28	5
SNP_1_11955697	1	11,955,697	−0.10	6.48	2.37	2
SNP_1_13379811	1	13,379,811	0.13	4.72	2.58	4
SNP_1_13975014	1	13,975,014	0.13	5.07	2.33	3
SNP_1_28484917	1	28,484,917	0.09, 013	4.37, 6.92	1.37, 4.19	2, 6
SNP_1_28959591	1	28,959,591	−0.18	5.71	4.95	4
SNP_1_29463236	1	29,463,236	−0.14	8.87	6.48	5
SNP_1_30876860	1	30,876,860	0.20	12.55	5.53	2
SNP_1_30877239	1	30,877,239	0.20	3.72	3.00	4
SNP_1_31421621	1	31,421,621	0.52	15.67	41.80	4
SNP_1_31536782	1	31,536,782	0.22	10.26	11.44	2
SNP_1_34323683	1	34,323,683	−0.11	7.25	3.03	2
SNP_2_6906642	2	6,906,642	0.09	3.78	2.28	2
SNP_2_7398675	2	7,398,675	−0.13	5.06	2.15	3
SNP_2_9450368	2	9,450,368	0.14	5.68	0.76	6
SNP_2_13001295	2	13,001,295	0.17	4.05	2.30	5
SNP_2_19542765	2	19,542,765	−0.05	3.77	0.31	6
SNP_2_33277520	2	33,277,520	0.09, 0.10	4.27, 5.33	1.26, 3.15	2, 5
SNP_2_34564546	2	34,564,546	0.17	6.32	2.53	6
SNP_2_34591582	2	34,591,582	−0.33, −0.19	6.53, 11.94	3.79, 6.35	3, 5
SNP_2_34977775	2	34,977,775	0.17	4.00	2.82	4
SNP_3_4919439	3	4,919,439	−0.12	6.93	1.38	2
SNP_3_11299789	3	11,299,789	0.19	4.41	3.65	4
SNP_3_13056272	3	13,056,272	−0.15	6.30	3.41	3
SNP_3_13447204	3	13,447,204	−0.14	6.78	2.51	6
SNP_3_24636405	3	24,636,405	−0.18, −0.14	5.89, 12.45	2.54, 11.56	5, 6
SNP_3_29244009	3	29,244,009	0.33	3.83	4.32	1
SNP_4_6190448	4	6,190,448	−0.10	4.92	2.04	3
SNP_4_6488350	4	6,488,350	0.10, 0.13	5.58, 8.43	2.06, 2.46	2, 6
SNP_4_16168608	4	16,168,608	0.07–0.28	3.74–8.40	1.50–4.47	1, 2, 5
SNP_4_16170098	4	16,170,098	0.11	6.80	1.87	6
SNP_4_23433069	4	23,433,069	0.12	6.50	3.21	5
SNP_4_23981762	4	23,981,762	0.17	6.36	3.12	4
SNP_4_31803015	4	31,803,015	0.09	4.01	0.62	2
SNP_4_32081836	4	32,081,836	−0.05	4.27	1.54	5
SNP_4_33137021	4	33,137,021	−0.14	6.92	1.43	2
SNP_5_1982590	5	1,982,590	0.21	8.60	5.67	4
SNP_5_2207493	5	2,207,493	0.18	9.12	4.91	5
SNP_5_9134248	5	9,134,248	−0.14	5.66	2.26	5
SNP_5_14757146	5	14,757,146	−0.16	3.92	3.72	4
SNP_5_15728722	5	15,728,722	−0.22	7.40	2.30	4
SNP_6_1491837	6	1,491,837	0.13	6.49	2.87	3
SNP_6_4806530	6	4,806,530	0.11	4.27	2.15	5
SNP_6_23051425	6	23,051,425	−0.12	4.96	2.04	5
SNP_7_7814446	7	7,814,446	0.08, 0.10	3.90, 6.89	0.87, 2.08	2, 6
SNP_7_8168133	7	8,168,133	−0.06	4.13	1.70	5
SNP_7_10448596	7	10,448,596	0.05	4.26	0.32	6
SNP_7_15970751	7	15,970,751	−0.27	8.09	6.51	4
SNP_7_16105614	7	16,105,614	0.11, 0.11	4.75, 5.09	0.77, 1.03	2, 6
SNP_7_17198761	7	17,198,761	0.24, 0.25	7.04, 7.78	3.65, 5.80	3, 5
SNP_7_19392618	7	19,392,618	0.15	6.87	1.02	6
SNP_7_25050720	7	25,050,720	−0.22 to −0.08	3.72–14.89	0.50–7.86	2, 3, 4, 5, 6
SNP_7_26120188	7	26,120,188	0.15	6.85	2.18	2
SNP_7_27379415	7	27,379,415	0.16	7.95	2.00	2
SNP_8_15994461	8	15,994,461	0.19, 0.25	7.58, 7.75	2.86, 3.33	3, 5
SNP_8_21207866	8	21,207,866	0.15	5.13	1.58	3
SNP_8_26869785	8	26,869,785	0.19	5.35	1.92	3
SNP_9_12628970	9	12,628,970	0.19, 0.22	8.03, 11.93	2.56, 2.65	2, 6
SNP_9_14359329	9	14,359,329	0.17	8.82	3.65	5
SNP_10_1274853	10	1,274,853	−0.05	3.77	0.21	6
SNP_10_11757992	10	11,757,992	0.08	6.36	1.09	2
SNP_10_14875442	10	14,875,442	−0.15	9.86	3.61	6
SNP_11_5489458	11	5489458	0.11	6.00	1.73	6
SNP_11_18003112	11	18,003,112	0.08	3.89	0.45	2
SNP_11_21472458	11	21,472,458	0.11	8.15	1.54	2
SNP_11_22469803	11	22,469,803	0.18	10.54	4.66	6
SNP_11_22976960	11	22,976,960	−0.12	7.23	1.04	6
SNP_11_23730796	11	23730796	0.10	5.14	1.43	2
SNP_11_27538579	11	27538579	0.07	4.46	0.49	6
SNP_12_13859027	12	13859027	−0.17	12.54	3.93	6
SNP_12_14597524	12	14597524	0.17	6.92	1.75	2

A comparison of the QTNs detected by different multi-locus methods identified that 12 QTNs were co-detected by more than one method. These QTNs were located with one each on Chromosomes 1, 3, 8, and 9, two each on Chromosomes 2 and 4, and four each on Chromosome 7. The QTN SNP_7_25050720 was detected by five multi-locus methods (pLARmEB, ISIS EM-BLASSO, FASTmrMLM, mrMLM, and pKWmEB), and the QTN SNP_4_16168608 was detected by three multi-locus methods (FASTmrEMMA, FASTmrMLM, and pKWmEB) (Table 1).

A comparison between QTNs identified using multi-locus GWA methods and QTLs detected through single-locus GWA methods revealed that 21 of the 72 QTNs were colocated with QTLs (Table S4 and Figure 4). For example, three QTNs (SNP_2_34564546, SNP_2_34591582, and SNP_2_3497775) were found in the QTL on Chromosome 2 spanning 34.53–34.82 Mb, which was detected by single-locus GWAS methods (both EMMAX and GEMMA). The QTN SNP_7_15970751, detected via the multi-locus mrMLM method, was located within the QTL on Chromosome 7, ranging from 15.39 to 16.07 Mb, which was also detected by single-locus GWAS methods.

[IMAGE OMITTED. SEE PDF]

QTL on Chromosome 7 at ~15.55 Mb

The QTL on Chromosome 7 from 15.39 to 16.07 Mb was detected using both the single-experiment and multi-experiment GWA methods and also detected by multi-locus GWAS. The effect size of peak SNP (SNP_7_15551740) for this QTL identified by the single-experiment GWA is 23.7% (Table S2). A significant QTN (SNP_7_15970751) was co-detected in this QTL by multi-locus GWAS; this QTN effect is 27%, and this QTN explains 6.51% of phenotypic variation. Four genes annotated as aquaporin proteins are found in this QTL region, namely, OsPIP2;4 (LOC_Os07g26630), OsPIP2;9 (LOC_Os07g26640), OsPIP2;5 (LOC_Os07g26660), and OsPIP2;1 (LOC_Os07g26690). Due to their close proximity, these genes are linked within the studied population. OsPIP2;4 (positioned at 15,358,960–15,360,495 bp) was selected as the gene to report as a candidate for this QTL based on a priori knowledge about its function. The transcript level of OsPIP2;4 in rice roots and shoots has been reported to be significantly downregulated in response to As(III) treatment, while overexpression of OsPIP2;4 in Arabidopsis thaliana (Arabidopsis) yielded enhanced As(III) tolerance, and in long-term uptake experiments, these transgenic plants showed no significant accumulation of As (Mosa et al. 2012). The study of Mosa et al. (2012) did not explore the expression of the other three PIP genes in this location (OsPIP2;1, OsPIP2;5, and OsPIP2;9); however, in a previous analysis in rice roots exposed to arsenate, none of these genes where either upregulated or downregulated (Norton et al. 2008).

A total of 42 SNPs were identified in OsPIP2;4 in the 253 RDP1 used in this study and 12 SNPs were in the exons of this gene, three of which were nonsynonymous SNPs. The physical positions of these three nonsynonymous SNPs from the beginning of the chromosome were at 15,359,866 bp (T/C polymorphism), 15,359,874 bp (G/A polymorphism), and 15,360,078 bp (G/A polymorphism), which resulted in the following amino acid changes: V237A, V240I, and V270I, respectively. The cultivars in this study carried two haplotypes for these 12 SNPs: I (n = 80) and II (n = 167) (Table S5). The cultivars carrying Haplotype I (mean As concentration: 0.41, 0.36, 0.24, and 0.58 mg kg⁻¹ for plants growing in Faridpur, Arkansas in 2006 and 2007, and Texas, respectively) had significantly (p < 0.05) lower grain As concentration compared to those carrying Haplotype II (mean As concentration: 0.46, 0.40, 0.28, and 0.68 mg kg⁻¹ for plants growing in Faridpur, Arkansas in 2006 and 2007, and Texas, respectively) across four of the five field studies, resulting in a reduction in grain As of 11.78%–15.78% (Figure 5). Only the grains harvested from plants at the Qiyang location did not exhibit a statistical difference in As concentration between the two haplotypes. Among the rice subpopulations, the temperate japonicas have, as a group, the lowest grain As (Norton et al. 2014; Pinson et al. 2015), and Haplotype I, associated with reduced grain As, is the predominant haplotype among the temperate japonica cultivars (83%) (Table S5 and Figure S1a). While Haplotype I is completely absent among the indica RDP1 cultivars, all other subpopulations had five or more cultivars carrying Haplotype I that can be used by breeders as gene donors.

[IMAGE OMITTED. SEE PDF]

As already mentioned, the four OsPIP genes in this region are in close physical proximity to each other; therefore, while the haplotype analysis of OsPIP2;4 has been presented, it is important to note that the cultivars that have the low grain arsenic haplotype, Haplotype I, of OsPIP2;4 also have a low grain arsenic haplotype in the other three OsPIP genes in this region. Of the 81 cultivars identified as having the low grain haplotype, Haplotype I, for OsPIP2;4, 76, 77, and 68 of these cultivars also have a low grain arsenic haplotypes for OsPIP2;9, OsPIP2;5, and OsPIP2;1, respectively. These haplotypes are Haplotype I for OsPIP2;9, Haplotype V for OsPIP2;5, and Haplotype IV for OsPIP2;1 (Figures 5 and S2 and Table S5).

QTN on Chromosome 4 at ~31,803,015 bp

The QTN SNP_4_31803015 was detected by the multi-locus GWA method (FASTmrMLM), explaining 0.62% of phenotypic variation, and a QTL (on Chromosome 4 at ~31.64 Mb) close to this QTN was detected using a single-locus method (multi-experiment GWA by GEMMA). The gene OsABCC1 (positioned at 31,502,940–31,514,308 bp), which is 0.29 kb from this QTN, was considered as a candidate gene, as the ABC transporter OsABCC1 has been reported to sequester As into the vacuoles of phloem companion cells in stem nodes and thus reduce As accumulation in rice grain (Song et al. 2014). A total of 55 SNPs were found in this candidate gene from the 253 cultivars, of which seven were nonsynonymous SNPs. The physical positions of these seven SNPs were at 31,507,580 (C/G polymorphism), 31,513,040 (C/T polymorphism), 31,513,380 bp (G/A polymorphism), 31,514,034 bp (G/A polymorphism), 31,514,212 bp (T/A polymorphism), 31,514,242 bp (T/A polymorphism), and 31,514,248 bp (C/A polymorphism), which result in the following amino acid changes: V814L, R276Q, P206L, A92V, T33S, N23Y, and V21L, respectively. The cultivars in this study carried eight haplotypes (n ≥ 5) for these seven SNPs (Table S6): I (n = 44), II (n = 16), III (n = 14), IV (n = 84), V (n = 13), VI (n = 16), VII (n = 56), and VIII (n = 5) (Figure 6). The cultivars that had Haplotype VII (mean = 0.56 mg kg⁻¹) showed significantly lower grain As concentration (p < 0.01) than other haplotypes when the plants were grown in the fields at Qiyang, China (means ranged from 0.68 to 0.80 mg kg⁻¹ for Haplotypes I–VI and VIII), resulting in a reduction in grain arsenic concentration by 21.4%–44.8%. With cultivars of the indica subpopulation being preferred throughout most of China, it is important to note that none of the indica cultivars within RDP1 contain Haplotype VII, which was instead common among tropical japonica cultivars (Table S6 and Figure S1b). In contrast, the cultivars carrying Haplotypes I (mean = 0.21 mg kg⁻¹) and VII (mean = 0.22 mg kg⁻¹) had lower grain As concentrations than other haplotypes (mean ranged from 0.25 to 0.41 mg kg⁻¹) when grown in the fields at Arkansas in 2007 and Texas 2009 (Figure 6). Rice varieties produced commercially in the US Mid-South, including Arkansas and Texas, are predominantly of the tropical japonica subpopulation. While the more desirable Haplotype I was significantly more common among temperate japonica cultivars, two of the 67 tropical japonica cultivars in RDP1 do contain Haplotype I, including the historically grown US medium grain cultivar “Saturn” and the US long grain “C57-5043.” In fact, C57-5043 contains the haplotype associated with reduced grain As in the US Mid-South at both the OsABCC1 and OsPIP2;4 genes.

[IMAGE OMITTED. SEE PDF]

It is especially interesting to note that the haplotype common among tropical japonica cultivars was associated with reduced grain As in the experiment conducted in China, where indica cultivars prevail, but was not the haplotype with lowest grain As in the experiments conducted in the US Mid-South, where tropical japonica cultivars prevail. This might indicate that the various haplotypes/alleles have different impacts on inorganic versus organic As. The soil in Qiyang, China, is known to be relatively high in inorganic As (Norton et al. 2014), while the As in rice grains produced in the US Mid-South is generally 50% or more organic As (Williams et al. 2005; F. Zhao et al. 2013).

Discussion

Comparison of Different GWAS Methods

In this study, eight GWAS methods, EMMAX, GEMMA, mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, pKWmEB, and ISIS EM-BLASSO, were used to identify grain As QTLs and QTNs. This approach has led to the identification of a large number of QTLs and QTNs for grain arsenic. EMMAX and GEMMA are single-locus approaches while the others are multi-locus approaches. EMMA has been a gold standard method for GWAS because it corrects for population structure and genetic relatedness in model organism association mapping, which can reduce false positives (Type I errors), but this approach is computationally impractical (Kang et al. 2008). Kang et al. (2010) modified EMMA and developed EMMAX whereby the estimation of variance components was approximated, therefore removing the requirement to estimate this for each marker. EMMAX has reduced computation time (Kang et al. 2010). GEMMA has been developed from the EMMA algorithm (Kang et al. 2008). The GEMMA model can fit both univariate linear mixed models and mvLMMs for testing association with a single trait and multiple traits, respectively (Zhou and Stephens 2014). However, most quantitative traits are controlled both by a few genes with major effects and by numerous polygenes with minor effects, which the single-locus GWAS methods fail to detect. Multi-locus GWAS methods were developed to overcome this issue, such as the mrMLM, FASTmrMLM, FASTmrEMMA, pLARmEB, pKWmEB, and ISIS EM-BLASSO (X. Liu et al. 2016; Ren et al. 2018; Tamba et al. 2017; Tamba and Zhang 2018; J. Zhang et al. 2017). All six of these methods involve two-step algorithms. During the first step, a single-locus GWAS method is applied to scan the entire genome, and putative QTNs are detected according to a less stringent critical value, such as p < 0.005. During the second step, all selected putative QTNs are examined by a multi-locus GWAS model to detect true QTNs (X. Liu et al. 2016; Ren et al. 2018; Tamba et al. 2017; Tamba and Zhang 2018; J. Zhang et al. 2017). Each of these methods has its own advantages and limitations. For example, FASTmrEMMA is more powerful in QTN detection and model fit, exhibits less bias in QTN effect estimation, and requires less running time, though its two-step process complicates the determination of suitable thresholds for marker selection in the first step (Wen et al. 2018; Y. Xu et al. 2018). The LASSO method performs well for detecting large-effect QTL but poorly for small-effect QTL and fails to handle a large number of markers (Y. Xu et al. 2018). ISIS EM-BLASSO reduced the number of markers via correlation learning to a moderate number and was then employed to select variables in the reduced model (Tamba et al. 2017). The results of this study were consistent with these reports, with the FASTmrEMMA detecting the most QTNs (23) and ISIS EM-BLASSO detecting fewer QTNs (11). To improve the power of GWAS, several studies used different GWAS methods to mine target QTLs for complex agronomic traits (Cui et al. 2018; Y. Xu et al. 2018; Y.K. Zhang et al. 2015; Zhong et al. 2021), with the logic that QTLs detected by multiple methods may be more reliable.

Colocalization of QTLs for Arsenic Concentration Across Studies

This study employed two single-locus GWAS methods (EMMAX and GEMMA) and six multi-locus GWAS methods to identified QTLs and QTNs associated with grain As concentration using data from five field experiments (Figures 2 and 3 and Table 1). The phenotypic data used here had previously been analyzed for GWAS using EMMA, but only using a 44K SNP array (Norton et al. 2014). In that earlier study, significant associations were identified across the entire RDP1 population as well as within each subpopulation (indica, aus, temperate japonica, and tropical japonica). From the whole RDP1, a total of 15 QTLs were identified for grain As concentration when based on the criterion of the QTLs having two or more significant SNPs. In the current study, which utilized 3.5M filtered SNPs, 90 QTLs were detected using EMMAX, with nine QTLs overlapping between the two studies. The absence of six QTLs identified in the 2014 study but not detected here could be attributed to several factors. First, while the same phenotypic data were used, the number of genotypes analyzed differed; Norton et al. (2014) included all genotypes with data available for individual experiment, whereas the current study was limited to 253 genotypes that have As data across all field sites, potentially affecting QTL detection. Second, the current study employed data transformation, which is recommended for highly skewed datasets to reduce false positives. This transformation might have influenced the detection of certain QTLs near significance threshold. It is important to note that the most statistically significant QTLs remained unaffected by the data transformation.

QTLs and QTNs (± 5 kb, the average LD decay of RDP1) identified in this study were compared with those reported in the previously published studies. This comparison that 85 QTLs (and QTNs) identified here colocated with QTLs previously reported from the studies of various other populations (Dong et al. 2023; X. Liu et al. 2019; H. Liu et al. 2021; Masato et al. 2013; Norton et al. 2019; Pinson et al. 2022; Murugaiyan et al. 2019; J. Zhang et al. 2008) (Figure 4). Several colocalized QTLs and QTNs were particularly noteworthy due to their detection by multiple methods in this study and their previous identification in other studies. For example, the QTL located on Chromosome 3 at ~24.5 Mb, detected here by both single-locus and multi-locus GWA methods, was previously reported in two studies (Norton et al. 2019; J. Zhang et al. 2008), though it was not detected by prior analysis of this same phenotypic data (Norton et al. 2014). The QTL on Chromosome 9 at ~20.8 Mb was detected here by single-locus (both EMMAX and GEMMA) GWAS and had been reported in three previous studies (H. Liu et al. 2021; Norton et al. 2019; Murugaiyan et al. 2019). The QTN located on Chromosome 4 at ~16.16 Mb, identified in this study by three multi-locus methods (FASTmrEMMA, FASTmrMLM, and pKWmEB), was previously reported by Norton et al. (2019). These results indicate that there are environmentally stable QTLs for grain As and the underlying genes could be important for breeding low grain As cultivars.

Candidate Genes for As

In this study, 10 candidate genes previously reported to be involved with As stress in rice were located within 500 kb of peak SNPs and/or QTNs identified (Table S7 and Figure 4). The candidate gene OsABCC1, a rice C-type ATP-binding cassette (ABC) transporter, was identified within the QTL on Chromosome 4 at ~31.64 Mb, close to the QTN SNP_4_31803015. OsABCC1 has been reported to be involved in the detoxification and reduction of As in rice grains by sequestering As in the vacuoles of the phloem companion cells of the nodes (Song et al. 2014).

Two candidate genes, OsNIP1;1 and OsNIP2;2 (OsLsi6), encode NIPs. OsNIP1;1 was found in the QTL on Chromosome 2 at ~7.1 Mb and a QTN (SNP_2_7398675), which colocated with the QTL spanning 7.02–7.28 Mb, by Murugaiyan et al. (2019). OsNIP2;2 (OsLsi6) was found in the QTL on Chromosome 6 from 6.54 to 6.61 Mb, which was also reported by Masato et al. (2013). Overexpression of OsNIP1;1 can restrict As(III) loading into the xylem, thereby reducing As accumulation in rice grain (Sun et al. 2018). OsNIP2;2 (Lsi6) is a silicon transporter in rice, directing silicon to either the seed husk or the flag leaf (Yamaji and Jian 2009); it mediates As(III) transport and may contribute to As loading into rice grain (Lindsay and Maathuis 2017).

A PIP, OsPIP1;3, was also identified. OsPIP1;3 was found in the QTN on Chromosome 2 at ~34.97 Mb, detected by the multi-locus GWA method, which colocated with the QTL identified in the BAAP population (Norton et al. 2019). This gene has been demonstrated to be significantly downregulated in rice root and shoot under As(III) stress (Mosa et al. 2012). In addition to OsPIP1;3, a cluster of PIP genes on Chromosome 7 was identified in the QTL on Chromosome 7 between 15.39 and 16.01 Mb; this QTL was detected by both single-locus (EMMAX and GEMMA) and multi-locus (mrMLM) GWAS methods. This region is interesting as it contains four PIP genes in close proximity. Haplotype analysis of the PIP genes in this region showed that there was a low grain arsenic haplotype (made up of predominantly japonica cultivars) that was present across all four genes. Therefore, it is not possible to separate out which of these is the most likely candidate gene based on allelic variation. However, it has been demonstrated that OsPIP2;4 was significantly downregulated in rice root and shoot under As(III) stress and that overexpression of OsPIP2;4 in Arabidopsis resulted in enhanced As(III) tolerance, with transgenic plants showed no significant accumulation of As (Mosa et al. 2012). The Mosa et al. (2012) study did not report on the other PIPs in this cluster, and it is noteworthy that the transcriptomic study of rice root response to arsenate by Norton et al. (2008) did not reveal any of these aquaporins as arsenic responsive. Based on the information above, any or all of these PIPs should be considered as potentially contributing to arsenic uptake in rice.

OsPT4 (OsPht1;4), a member of the Pht1 family of phosphate transporters, was identified as the candidate gene for a QTL on Chromosome 4 at ~6.19 Mb. The Pht1 family of genes encodes inorganic phosphate (Pi) transporters localized to the plasma membrane. These transports have been implicated in As uptake and transport in rice (Kamiya et al. 2013). Two candidate genes, OsAIR1 and OsAIR2, are Oryza sativa arsenic-induced RING E3 ligase genes. OsAIR1 was the candidate gene for the QTL on Chromosome 5 ranging from 24.00 to 24.78 Mb, detected via EMMAX. OsAIR2 was the candidate gene of the QTL (QTN) on Chromosome 11 at ~21.47 Mb, identified by both single-locus and multi-locus GWA methods. Both OsAIR1 and OsAIR2 were induced under As exposure (Hwang et al. 2017, 2016). Overexpression of OsAIR1 in rice influenced seeding root length and cotyledon expansion under As(V) expression (Hwang et al. 2016). Heterogeneous overexpression of OsAIR2 in Arabidopsis enhanced seed germination and increased root length under As(V) stress conditions (Hwang et al. 2017).

Conclusion

In this study, two single-locus GWAS methods (EMMA and GEMMA) and six multi-locus GWAS methods (FASTmrEMMA, FASTmrMLM, ISIS EM-BLASSO, mrMLM, pKWmEB, and pLARmEB) were used to identify the loci associated with grain As accumulation across five field experiments. There were 90 and 111 QTLs for grain As detected by using the EMMAX and GEMMA, respectively. A total of 72 QTNs for grain As were detected from six multi-locus GWAS methods. Of these, 21 QTLs (QTNs) were co-detected by both single-locus and multi-locus GWAS methods. A total of 10 candidate genes were identified. The rice ABC transporter OsABCC1 was detected by GEMMA and FASTmrMLM GWAS methods. This approach of using an array of statistical methods and focusing on loci that appear in more than one study has the advantage of not being reliant on a single set of computational criteria. These QTLs, QTNs, and candidate genes represent promising targets for further genetic and breeding efforts aimed at reducing As accumulation in rice grain.

Author Contributions

C.C., A.H.P., G.J.N.: conceptualization; C.C., P.R., A.J.T., A.D.: methodology; C.C., A.D., S.R.M.P., G.C.E., A.H.P., G.J.N.: formal analysis; D.E.S., S.R.M.P., G.C.E.: resources; C.C.: data curation; C.C.: writing – original draft; C.C., P.R., A.D., A.J.T., D.E.S., S.R.M.P., G.C.E., A.H.P., G.J.N.: writing – review & editing; C.C., G.C.E.: visualization.

Acknowledgments

We acknowledge Professor Mahendra Modi from Assam Agricultural University, India, for providing access to the high-performance computer Aryabhatta for analysis of the data.

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

The phenotype data used to conduct the GWA mapping are available in the Supporting Information (Table S1).

References

Benjamini, Y., and Y. Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society, Series B 57: 289–300. https://doi.org/10.1111/j.2517‐6161.1995.tb02031.x.

Chen, C., G. J. Norton, and A. H. Price. 2020. “Genome‐Wide Association Mapping for Salt Tolerance of Rice Seedlings Grown in Hydroponic and Soil Systems Using the Bengal and Assam Aus Panel.” Frontiers in Plant Science 11: 1–18. https://doi.org/10.3389/fpls.2020.576479.

Chen, C., A. J. Travis, M. Hossain, R. Islam, A. H. Price, and G. J. Norton. 2021. “Genome‐Wide Association Mapping of Sodium and Potassium Concentration in Rice Grains and Shoots Under Alternate Wetting and Drying and Continuously Flooded Irrigation.” Theoretical and Applied Genetics 134: 2315–2334. https://doi.org/10.1007/s00122‐021‐03828‐9.

Codex Alimentarius. 2017. “General Standard for Contaminants and Toxins in Food and Feed.” Codex Alimentarius, Rome. https://www.fao.org/fao‐who‐codexalimentarius.

Cui, Y., F. Zhang, and Y. Zhou. 2018. “The Application of Multi‐Locus GWAS for the Detection of Salt‐Tolerance Loci in Rice.” Frontiers in Plant Science 9: 1–9. https://doi.org/10.3389/fpls.2018.01464.

Dong, Z., L. Guo, X. Li, et al. 2023. “Genome‐Wide Association Study of Arsenic Accumulation in Polished Rice.” Genes (Basel) 14: 2186. https://doi.org/10.3390/genes14122186.

Eizenga, G. C., M. L. Ali, R. J. Bryant, K. M. Yeater, A. M. McClung, and S. R. McCouch. 2014. “Registration of the Rice Diversity Panel 1 for Genomewide Association Studies.” Journal of Plant Registrations 8: 109–116. https://doi.org/10.3198/jpr2013.03.0013crmp.

Etesami, H., B. R. Jeong, F. J. M. Maathuis, and J. Schaller. 2023. “Exploring the Potential: Can Arsenic (As) Resistant Silicate‐Solubilizing Bacteria Manage the Dual Effects of Silicon on As Accumulation in Rice?” Science of the Total Environment 903: 166870. https://doi.org/10.1016/j.scitotenv.2023.166870.

Huang, X., X. Wei, T. Sang, et al. 2010. “Genome‐Wide Association Studies of 14 Agronomic Traits in Rice Landraces.” Nature Genetics 42: 961–967. https://doi.org/10.1038/ng.695.

Huang, X., Y. Zhao, X. Wei, et al. 2012. “Genome‐Wide Association Study of Flowering Time and Grain Yield Traits in a Worldwide Collection of Rice Germplasm.” Nature Genetics 44: 32–39. https://doi.org/10.1038/ng.1018.

Hwang, S., S. Chapagain, A. Han, et al. 2017. “Molecular Characterization of Rice Arsenic‐Induced RING Finger E3 Ligase 2 (OsAIR2) and Its Heterogeneous Overexpression in Arabidopsis thaliana.” Physiologia Plantarum 161: 372–384. https://doi.org/10.1111/ppl.12607.

Hwang, S. G., H. M. Park, A. R. Han, and C. S. Jang. 2016. “Molecular Characterization of Oryza sativa Arsenic‐Induced RING E3 Ligase 1 (OsAIR1): Expression Patterns, Localization, Functional Interaction, and Heterogeneous Overexpression.” Journal of Plant Physiology 191: 140–148. https://doi.org/10.1016/j.jplph.2015.12.010.

Jian, F. M., N. Yamaji, N. Mitani, et al. 2008. “Transporters of Arsenite in Rice and Their Role in Arsenic Accumulation in Rice Grain.” Proceedings of the National Academy of Sciences of the United States of America 105: 9931–9935. https://doi.org/10.1073/pnas.0802361105.

Kaler, A. S., J. D. Gillman, T. Beissinger, and L. C. Purcell. 2020. “Comparing Different Statistical Models and Multiple Testing Corrections for Association Mapping in Soybean and Maize.” Frontiers in Plant Science 10: 1–13. https://doi.org/10.3389/fpls.2019.01794.

Kamiya, T., M. R. Islam, G. Duan, S. Uraguchi, and T. Fujiwara. 2013. “Phosphate Deficiency Signaling Pathway Is a Target of Arsenate and Phosphate Transporter OsPT1 Is Involved in As Accumulation in Shoots of Rice.” Soil Science & Plant Nutrition 59: 580–590. https://doi.org/10.1080/00380768.2013.804390.

Kang, H. M., J. H. Sul, Service SK, et al. 2010. “Variance Component Model to Account for Sample Structure in Genome‐Wide Association Studies.” Nature Genetics 42: 348–354. https://doi.org/10.1038/ng.548.

Kang, H. M., N. A. Zaitlen, C. M. Wade, et al. 2008. “Efficient Control of Population Structure in Model Organism Association Mapping.” Genetics 178: 1709–1723. https://doi.org/10.1534/genetics.107.080101.

Karagas, M. R., T. Punshon, M. Davis, et al. 2019. “Rice Intake and Emerging Concerns on Arsenic in Rice: A Review of the Human Evidence and Methodologic Challenges.” Current Environmental Health Reports 6: 361–372. https://doi.org/10.1007/s40572‐019‐00249‐1.

Korte, A., and A. Farlow. 2013. “The Advantages and Limitations of Trait Analysis With GWAS: A Review Self‐Fertilization Makes Arabidopsis Particularly Well Suited to GWAS.” Plant Methods 9: 29.

Lemon, J. 2006. “Plotrix: A Package in the Red Light District of R.” R‐News 6: 8–12.

Li, L., X. Wu, J. Chen, et al. 2022. “Genetic Dissection of Epistatic Interactions Contributing Yield‐Related Agronomic Traits in Rice Using the Compressed Mixed Model.” Plants 11: 2504. https://doi.org/10.3390/plants11192504.

Lindsay, E. R., and F. J. M. Maathuis. 2017. “New Molecular Mechanisms to Reduce Arsenic in Crops.” Trends in Plant Science 22: 1016–1026. https://doi.org/10.1016/j.tplants.2017.09.015.

Liu, H., S. X. Long, S. R. M. Pinson, et al. 2021. “Univariate and Multivariate QTL Analyses Reveal Covariance Among Mineral Elements in the Rice Ionome.” Frontiers in Genetics 12: 638555. https://doi.org/10.3389/fgene.2021.638555.

Liu, X., S. Chen, M. Chen, et al. 2019. “Association Study Reveals Genetic Loci Responsible for Arsenic, Cadmium and Lead Accumulation in Rice Grain in Contaminated Farmlands.” Frontiers in Plant Science 10: 1–16. https://doi.org/10.3389/fpls.2019.00061.

Liu, X., M. Huang, B. Fan, E. S. Buckler, and Z. Zhang. 2016. “Iterative Usage of Fixed and Random Effect Models for Powerful and Efficient Genome‐Wide Association Studies.” PLoS Genetics 12: e1005767. https://doi.org/10.1371/journal.pgen.1005767.

Liu, Z., H. Li, Z. Gou, et al. 2020. “Genome‐Wide Association Study of Soybean Seed Germination Under Drought Stress.” Molecular Genetics and Genomics 295: 661–673. https://doi.org/10.1007/s00438‐020‐01646‐0.

Masato, K., T. Abe, A. Kawasaki, et al. 2013. “Correction to: Genetic Diversity of Arsenic Accumulation in Rice and QTL Analysis of Methylated Arsenic in Rice Grains (Rice, (2013), 6).” Rice 11, no. 1: 29. https://doi.org/10.1186/s12284‐018‐0221‐6.

McCouch, S. R., M. H. Wright, C. W. Tung, et al. 2016. “Open Access Resources for Genome‐Wide Association Mapping in Rice.” Nature Communications 7: 10532. https://doi.org/10.1038/ncomms10532.

McNally, K. L., K. L. Childs, R. Bohnert, et al. 2009. “Genomewide SNP Variation Reveals Relationships Among Landraces and Modern Varieties of Rice.” Proceedings of the National Academy of Sciences of the United States of America 106: 12273–12278. https://doi.org/10.1073/pnas.0900992106.

Meharg, A. A., and M. Rahman. 2003. “Arsenic Contamination of Bangladesh Paddy Field Soils: Implications for Rice Contribution to Arsenic Consumption.” Environmental Science & Technology 37: 229–234. https://doi.org/10.1021/es0259842.

Mitani, N., N. Yamaji, and J. F. Ma. 2008. “Characterization of Substrate Specificity of a Rice Silicon Transporter, Lsi1.” Pflügers Archiv: European Journal of Physiology 456: 679–686. https://doi.org/10.1007/s00424‐007‐0408‐y.

Mitani‐Ueno, N., N. Yamaji, F.‐J. Zhao, and J. F. Ma. 2011. “The Aromatic/Arginine Selectivity Filter of NIP Aquaporins Plays a Critical Role in Substrate Selectivity for Silicon, Boron, and Arsenic.” Journal of Experimental Botany 62: 4391–4398. https://doi.org/10.1093/jxb/err158.

Mosa, K. A., K. Kumar, S. Chhikara, et al. 2012. “Members of Rice Plasma Membrane Intrinsic Proteins Subfamily Are Involved in Arsenite Permeability and Tolerance in Plants.” Transgenic Research 21: 1265–1277. https://doi.org/10.1007/s11248‐012‐9600‐8.

Murugaiyan, V., J. Ali, A. Mahender, et al. 2019. “Mapping of Genomic Regions Associated With Arsenic Toxicity Stress (Oryza sativa L.).” Rice 12: 61. https://doi.org/10.1186/s12284‐019‐0321‐y.

Norton, G. J., A. Douglas, B. Lahner, et al. 2014. “Genome Wide Association Mapping of Grain Arsenic, Copper, Molybdenum and Zinc in Rice (Oryza sativa L.) Grown at Four International Field Sites.” PLoS ONE 9: e89685. https://doi.org/10.1371/journal.pone.0089685.

Norton, G. J., D. E. Lou‐Hing, A. A. Meharg, and A. H. Price. 2008. “Rice–Arsenate Interactions in Hydroponics: Whole Genome Transcriptional Analysis.” Journal of Experimental Botany 59: 2267–2276. https://doi.org/10.1093/jxb/ern097.

Norton, G. J., A. J. Travis, A. Douglas, et al. 2018. “Genome Wide Association Mapping of Grain and Straw Biomass Traits in the Rice Bengal and Assam Aus Panel (BAAP) Grown Under Alternate Wetting and Drying and Permanently Flooded Irrigation.” Frontiers in Plant Science 9: 1–18. https://doi.org/10.3389/fpls.2018.01223.

Norton, G. J., A. J. Travis, P. Talukdar, et al. 2019. “Genetic Loci Regulating Arsenic Content in Rice Grains When Grown Flooded or Under Alternative Wetting and Drying Irrigation.” Rice 12: 54. https://doi.org/10.1186/s12284‐019‐0307‐9.

Patterson, N., A. L. Price, and D. Reich. 2006. “Population Structure and Eigenanalysis.” PLoS Genetics 2: 2074–2093. https://doi.org/10.1371/journal.pgen.0020190.

Pinson, S. R. M., D. J. Heuschele, J. D. Edwards, A. K. Jackson, S. Sharma, and J. Y. Barnaby. 2022. “Relationships Among Arsenic‐Related Traits, Including Rice Grain Arsenic Concentration and Straighthead Resistance, as Revealed by Genome‐Wide Association.” Frontiers in Genetics 12: 787767. https://doi.org/10.3389/fgene.2021.787767.

Pinson, S. R. M., L. Tarpley, W. Yan, et al. 2015. “World‐Wide Genetic Diversity for Mineral Element Concentrations in Rice Grain.” Crop Science 55: 294–311. https://doi.org/10.2135/cropsci2013.10.0656.

Rahman, M. S., A. H. M. S. Reza, M. A. B. Siddique, M. A. Akbor, and M. Hasan. 2023. “Accumulation of Arsenic and Other Metals in Soil and Human Consumable Foods of Meherpur District, Southwestern Bangladesh, and Associated Health Risk Assessment.” Environmental Sciences Europe 35: 47. https://doi.org/10.1186/s12302‐023‐00751‐2.

Ren, W. L., Y. J. Wen, J. M. Dunwell, and Y. M. Zhang. 2018. “PKWmEB: Integration of Kruskal–Wallis Test With Empirical Bayes Under Polygenic Background Control for Multi‐Locus Genome‐Wide Association Study.” Heredity 120: 208–218. https://doi.org/10.1038/s41437‐017‐0007‐4.

Ruang‐areerate, P., A. J. Travis, S. R. M. Pinson, et al. 2020. “Genome‐Wide Association Mapping for Grain Manganese in Rice (Oryza sativa L.) Using a Multi‐Experiment Approach.” Heredity 126: 505–520. https://doi.org/10.1038/s41437‐020‐00390‐w.

Shi, S., T. Wang, Z. Chen, et al. 2016. “OsHAC1;1 and OsHAC1;2 Function as Arsenate Reductases and Regulate Arsenic Accumulation.” Plant Physiology 172: 1708–1719. https://doi.org/10.1104/pp.16.01332.

Song, W. Y., T. Yamaki, N. Yamaji, et al. 2014. “A Rice ABC Transporter, OsABCC1, Reduces Arsenic Accumulation in the Grain.” Proceedings of the National Academy of Sciences of the United States of America 111: 15699–15704. https://doi.org/10.1073/pnas.1414968111.

Sun, S. K., Y. Chen, J. Che, et al. 2018. “Decreasing Arsenic Accumulation in Rice by Overexpressing OsNIP1;1 and OsNIP3;3 Through Disrupting Arsenite Radial Transport in Roots.” New Phytologist 219: 641–653. https://doi.org/10.1111/nph.15190.

Tamba, C. L., Y. L. Ni, and Y. M. Zhang. 2017. “Iterative Sure Independence Screening EM‐Bayesian LASSO Algorithm for Multi‐Locus Genome‐Wide Association Studies.” PLoS Computational Biology 13: 1–20. https://doi.org/10.1371/journal.pcbi.1005357.

Tamba, C. L., and Y. M. Zhang. 2018. “A Fast mrMLM Algorithm for Multi‐Locus Genome‐Wide Association Studies.” bioRxiv. https://doi.org/10.1101/341784.

Turner, S. D. 2014. “qqman: An R Package for Visualizing GWAS Results Using QQ and Manhattan Plots.” bioRxiv. https://doi.org/10.1101/005165.

Wang, D. R., F. J. Agosto‐Pérez, D. Chebotarov, et al. 2018. “An Imputation Platform to Enhance Integration of Rice Genetic Resources.” Nature Communications 9: 1–10. https://doi.org/10.1038/s41467‐018‐05538‐1.

Wang, K., M. Li, and H. Hakonarson. 2010. “ANNOVAR: Functional Annotation of Genetic Variants From High‐Throughput Sequencing Data.” Nucleic Acids Research 38: 1–7. https://doi.org/10.1093/nar/gkq603.

Wang, M., J. Yan, J. Zhao, et al. 2012. “Genome‐Wide Association Study (GWAS) of Resistance to Head Smut in Maize.” Plant Science 196: 125–131. https://doi.org/10.1016/j.plantsci.2012.08.004.

Wang, P., W. Zhang, C. Mao, G. Xu, and F. J. Zhao. 2016. “The Role of OsPT8 in Arsenate Uptake and Varietal Difference in Arsenate Tolerance in Rice.” Journal of Experimental Botany 67: 6051–6059. https://doi.org/10.1093/jxb/erw362.

Wang, S. B., J. Y. Feng, W. L. Ren, et al. 2016. “Improving Power and Accuracy of Genome‐Wide Association Studies via a Multi‐Locus Mixed Linear Model Methodology.” Scientific Reports 6: 1–10. https://doi.org/10.1038/srep19444.

Wang, W., R. Mauleon, Z. Hu, et al. 2018. “Genomic Variation in 3,010 Diverse Accessions of Asian Cultivated Rice.” Nature 557: 43–49. https://doi.org/10.1038/s41586‐018‐0063‐9.

Wei, T., V. R. Simko, M. Levy, Y. Xie, Y. Jin, and J. Zemla. 2021. “Package “Corrplot”: Visualization of a Correlation Matrix.” Version 0.84.

Wen, Y. J., H. Zhang, Y. L. Ni, et al. 2018. “Methodological Implementation of Mixed Linear Models in Multi‐Locus Genome‐Wide Association Studies.” Briefings in Bioinformatics 19: 700–712. https://doi.org/10.1093/bib/bbw145.

Wickham, H. 2007. “Reshaping Data With the Reshape Package.” Journal of Statistical Software 21: 1–20.

Wickham, H. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer‐Verlag.

Williams, P. N., M. R. Islam, E. E. Adomako, et al. 2006. “Increase in Rice Grain Arsenic for Regions of Bangladesh Irrigating Paddies With Elevated Arsenic in Groundwaters.” Environmental Science & Technology 40: 4903–4908. https://doi.org/10.1021/es060222i.

Williams, P. N., A. H. Price, A. Raab, S. A. Hossain, J. Feldmann, and A. A. Meharg. 2005. “Variation in Arsenic Speciation and Concentration in Paddy Rice Related to Dietary Exposure.” Environmental Science & Technology 39: 5531–5540. https://doi.org/10.1021/es0502324.

Williams, P. N., A. Villada, C. Deacon, et al. 2007. “Greatly Enhanced Arsenic Shoot Assimilation in Rice Leads to Elevated Grain Levels Compared to Wheat and Barley.” Environmental Science & Technology 41, no. 19: 6854–6859.

Xu, J., S. Shi, L. Wang, et al. 2017. “OsHAC4 is Critical for Arsenate Tolerance and Regulates Arsenic Accumulation in Rice.” New Phytologist 215: 1090–1101. https://doi.org/10.1111/nph.14572.

Xu, Y., T. Yang, Y. Zhou, et al. 2018. “Genome‐Wide Association Mapping of Starch Pasting Properties in Maize Using Single‐Locus and Multi‐Locus Models.” Frontiers in Plant Science 9: 1–10. https://doi.org/10.3389/fpls.2018.01311.

Yamaji, N., and F. M. Jian. 2009. “A Transporter at the Node Responsible for Intervascular Transfer of Silicon in Rice.” Plant Cell 21: 2878–2883. https://doi.org/10.1105/tpc.109.069831.

Yang, M., K. Lu, F. J. Zhao, et al. 2018. “Genome‐Wide Association Studies Reveal the Genetic Basis of Ionomic Variation in Rice.” Plant Cell 30: 2720–2740. https://doi.org/10.1105/tpc.18.00375.

Ye, Y., P. Li, T. Xu, et al. 2017. “Ospt4 Contributes to Arsenate Uptake and Transport in Rice.” Frontiers in Plant Science 8: 1–12. https://doi.org/10.3389/fpls.2017.02197.

Zhang, J., J. Y. Feng, Y. L. Ni, et al. 2017. “PLARmEB: Integration of Least Angle Regression With Empirical Bayes for Multilocus Genome‐Wide Association Studies.” Heredity (Edinb) 118: 517–524. https://doi.org/10.1038/hdy.2017.8.

Zhang, J., Y. G. Zhu, D. L. Zeng, W. D. Cheng, Q. Qian, and G. L. Duan. 2008. “Mapping Quantitative Trait Loci Associated With Arsenic Accumulation in Rice (Oryza sativa) (New Phytologist (2008) 177, (350–355)).” New Phytologist 178: 688. https://doi.org/10.1111/j.1469‐8137.2008.02402.x.

Zhang, Y. K., D. F. Zhu, Y. P. Zhang, H. Z. Chen, J. Xiang, and X. Q. Lin. 2015. “Low pH‐Induced Changes of Antioxidant Enzyme and ATPase Activities in the Roots of Rice (Oryza sativa L.) Seedlings.” PLoS ONE 10: 1–12. https://doi.org/10.1371/journal.pone.0116971.

Zhang, Y. M., Z. Jia, and J. M. Dunwell. 2019. “Editorial: The Applications of New Multi‐Locus GWAS Methodologies in the Genetic Dissection of Complex Traits.” Frontiers in Plant Science 10: 1–6. https://doi.org/10.3389/fpls.2019.00100.

Zhang, Y. W., C. L. Tamba, Y. J. Wen, et al. 2020. “mrMLM v4.0.2: An R Platform for Multi‐Locus Genome‐Wide Association Studies.” Genomics, Proteomics & Bioinformatics 18: 481–487. https://doi.org/10.1016/j.gpb.2020.06.006.

Zhao, F. J., S. P. McGrath, and A. A. Meharg. 2010. “Arsenic as a Food Chain Contaminant: Mechanisms of Plant Uptake and Metabolism and Mitigation Strategies.” Annual Review of Plant Biology 61: 535–559. https://doi.org/10.1146/annurev‐arplant‐042809‐112152.

Zhao, F. J., Y. G. Zhu, and A. A. Meharg. 2013. “Methylated Arsenic Species in Rice: Geographical Variation, Origin, and Uptake Mechanisms.” Environmental Science & Technology 47: 3957–3966. https://doi.org/10.1021/es304295n.

Zhao, K., C. W. Tung, G. C. Eizenga, et al. 2011. “Genome‐Wide Association Mapping Reveals a Rich Genetic Architecture of Complex Traits in Oryza sativa.” Nature Communications 2: 1–10. https://doi.org/10.1038/ncomms1467.

Zhong, H., S. Liu, T. Sun, et al. 2021. “Multi‐Locus Genome‐Wide Association Studies for Five Yield‐Related Traits in Rice.” BMC Plant Biology 21: 1–12. https://doi.org/10.1186/s12870‐021‐03146‐8.

Zhou, X., and M. Stephens. 2014. “Efficient Multivariate Linear Mixed Model Algorithms for Genome‐Wide Association Studies.” Nature Methods 11: 407–409. https://doi.org/10.1038/nmeth.2848.

Zhu, Y. G., G. X. Sun, M. Lei, et al. 2008. “High Percentage Inorganic Arsenic Content of Mining Impacted and Nonimpacted Chinese Rice.” Environmental Science & Technology 42: 5008–5013. https://doi.org/10.1021/es8001103.