Trait-to-trait relationships and QTL were identified using phenotypic and genotypic data on the USDA-ARS Rice Minicore Collection (RMC), a set of 202 pure line rice (Oryza sativa) accessions collected from 14 global rice-growing regions (Agrama et al., 2009). The RMC was previously used for GWAS using a marker map consisting of 156 PCR-based molecular markers distributed across the twelve rice chromosomes, including identification of hull-Si QTL (Bryant et al., 2011). Subsequent to that study, 173 RMC accessions were resequenced (Wang et al., 2016), and a marker map containing 3,200,320 (3.2 million) filtered, reliable SNPs was created and used for GWA mapping (Huggins et al., 2019). The present study accumulated phenotypic data on the subset of 167 RMC accessions having a high-density marker map listed in Supplementary Table S1 which includes accessions from six genetic subpopulations in rice, specifically 30 accessions classified as being from the aus (AUS) subpopulation, 54 being indica (IND), 29 temperate japonica (TEJ), 28 tropical japonica (TRJ), and 6 aromatic (ARO), with the remaining 20 accessions being admixtures of two or more subpopulation groups. With the AUS, IND, and AUS/IND admixtures being considered members of the indica subspecies (INDAUS), there were 89 INDAUS accessions. Similarly, the study set includes 66 accessions of the japonica subspecies, identified as TEJ, TRJ, or TEJ/TRJ admixtures.

Bryant et al. (2011) evaluated the hull silica concentrations of RMC accessions using rice harvested from 2 replications × 2 locations (Stuttgart, AR, and Beaumont, TX) then used least square means (LSMeans) along with a low-density marker map (164 DNA markers) to identify 12 putative hull-Si QTL. Raw sample data were obtained and reanalyzed for the present study. As a subset of the USDA Rice Core Collection, the RMC was evaluated by Pinson et al. (2015) for grain concentrations of 16 elements, including As, P, S, Ca, and Cu, using grains harvested from 2 replications × 2 years × 1 location (Beaumont, TX) from both flooded and unflooded fields. For the present GWA analyses, raw sample data from flooded plots were obtained for grain concentrations (mg kg⁻¹) of As, P, S, Ca, and Cu, hereafter called grain-As, grain-P, and so forth. Details of sample production, separation of hulls and grains, sample digestions, and laboratory equipment and methods for determining element concentrations are in Bryant et al. (2011) and Pinson et al. (2015).

While some of the RMC accessions had been included in previous StHD QTL mapping studies (Li et al., 2016; Li et al., 2017), numerous RMC accessions had not yet been characterized for StHD. The RMC was evaluated for resistance to MSMA-induced StHD in the USDA-ARS Straighthead Nursery at Stuttgart, AR over 2 years, using two replications planted on May 5, 2015, and four replications planted on May 30, 2016. Growing rice accessions with red- or purple-colored bran in the Stuttgart, AR Straighthead Nursery is prohibited to prevent pollen outcrossing into nearby seed production fields. We therefore evaluated StHD on the 156 brown bran RMC accessions using 2-row plots drill-seeded with 3 g seed per 1.5 m row to ensure a dense plant stand, with 0.3 m between rows and plots. Plots were arranged in a randomized complete block design, with plots of a StHD-susceptible (“Cocodrie,” PI 606331, Linscombe et al., 2000) and resistant (“Zhe733,” PI 629016) cultivar inserted approximately every 10 plots in order to assure consistent StHD symptoms throughout the study area. The same RMC were similarly planted concurrently in an adjacent field area not treated with “MSMA,” hereafter called Native soil plots. Monosodium methanearsonate (MSMA) was applied on the day of planting at a rate of 6.7 kg ha⁻¹ using Target^® 6.6 (Luxembourg-Pamol, Inc., Houston, TX, United States) which corresponds with 1.6 kg As ha⁻¹ and incorporated into the upper 15 cm of soil before planting. Fertilization of fields began a month prior to planting with incorporation of P (triple super phosphate, 20 kg P ha⁻¹) and potassium (K) (muriate of potash, 56 kg K ha⁻¹). Nitrogen fertilizer (112 kg N ha⁻¹ as dry urea) was applied to the soil surface just before the permanent flood was established at the five-leaf stage. The flood was maintained throughout the growing season, until maturity of all plots, to ensure ideal conditions for StHD development. Weeds were controlled with 9.3 l ha⁻¹ of propanil (3′,4′-dichloropropionanilide) mixed with 0.4 kg ha⁻¹ of Quinclorac (3,7-dichloroquinoline-8-carboxylic acid; Facet^®, BASF, Ludwigshafen, Germany) when the rice was at the four-leaf stage. The soil at the site is Dewitt silt loam (fine, smectitic, thermic, Typic Albaqualfs) (5% sand, 78% silt, 17% clay), and all fields in this study had been managed in a rice-soybean rotation for more than 20 years. With the Straighthead nursery area receiving a fresh application of MSMA prior to planting rice every other year, the total As (iAs + oAs) in the soil has increased over time to 13 ± 6 kg ha⁻¹ (Maguffin et al., 2020) compared to 4 ± 2 kg As ha⁻¹ in the Native soil field area (Norton et al., 2012b). Regardless of the increased soil As content, a fresh application of MSMA is required to ensure development of StHD disorder in the nursery plots.

The MSMA-treated and Native soil plots were rated for StHD at maturity, approximately 35–40 days after heading. Heading was defined as 50% of the plants per plot having at least 1 panicle at anthesis. Days to heading (DHD) is the number of days between planting and heading per plot, with DHDms denoting data from MSMA-treated plots, and DHDnt for data from Native soil plots. Straighthead severity was rated on a 0 to 9 scale based on a visual observation of floret sterility and panicle emergence from the flag leaf sheath, as described by Yan et al. (2005). The plant height (PHT) of three random plants per plot was measured also at maturity and averaged per plot. The PHT for this study was defined as the distance (cm) from the soil surface to the tip of the tallest leaf to accommodate the fact that panicles are often not fully exerted on MSMA-treated plants. As for DHD, PHTms denotes data from MSMA-treated plots, and PHTnt denotes data from the Native soil plots. The 2015 trial was planted using RMC accessions obtained from the USDA Genetic Stocks Oryza (GSOR) Collection (www.ars.usda.gov/GSOR); 2016 plantings used seed harvested from 2015 Native soil plots.

Analyses of variance (ANOVAs) were conducted using the generalized linear mixed model (GLIMMIX) procedure in SAS version 9.4 (SAS Institute Inc., 2012). Least square means (LSmeans) were calculated using the GLIMMIX procedure, with replication considered as a random effect and the RMC accessions as fixed effects. For calculating LSmeans combined across environments (either locations or years), the statistical model considered environment and replications nested within the environment as random effects.

Trait distributions, summary statistics, and mean comparisons across the seven population panels were calculated in JMP 14 (SAS Institute Inc. 2018) using trait LSmeans. Best linear unbiased predictions (BLUPs) computed per trait from raw replication data were used for the Pearson correlations and Bayesian network analyses. After a review of the data distributions, it was decided to transform the grain element data to reduce the strong skewing of these trait datasets prior to GWA analyses. Raw element concentration data were log-transformed and new BLUPs calculated for GWAS.

Bayesian network learning is a multivariate probabilistic modeling which computes, through an iterative learning process, the relationship between random variables represented by “nodes” and probabilistic dependencies represented by “arrows” in a directed acyclic graph (DAG) (Pearl, 1988; Scutari, 2010). The DAG was computed using the R package “bnlearn” (Scutari, 2010), as described in more detail in Sharma et al. (2021). The Bayesian network (BN) model was trained with 10-fold cross validation to compute the DAG (Scutari, 2010). Each validation used a random 17 accessions as the validation set (VS = 17) with the rest of the RMC accessions serving as the training set (TS = 150). Trait BLUPs were computed anew using the lme4 R package (Bates et al., 2011), this time accounting for population structure by using a mixed-effect model with a simple nested family structure (Piepho and Williams, 2006) that considered subgroups in the RMC as blocks. The growth stage progression of the traits in the order of flowering (DHDnt), plant height (PHTnt), grain element compositions, and lastly StHD disease ratings (StHDms) was encoded in the BN model using blacklisted arcs to compute the DAG (Scutari 2010).

The “bnlearn” package learns by iteratively selecting and estimating the model. The procedure first uses feature selection to find “parents” and “children” relationships among traits in the Markov blanket by the Semi-Interleaved HITON PC algorithm. In our model, the dependence was analyzed using Student’s t-test for Pearson’s correlations with alpha at 0.1. The large values of alpha were used because it allows Markov blankets to initially involve traits that are weakly associated, so that they are not initially discarded. Next, the procedure uses structure learning to find the DAG by computing the conditional independence present in data with a score-based algorithm using a heuristic optimization technique. In each iteration, the candidate DAG is assigned a network score reflecting its goodness of fit, which the algorithm then attempts to maximize. Then across all iterations, the analysis selects the DAG that maximizes the Bayesian information criterion.

GWA was conducted using 7 different “panels” derived from the RMC data, including the full RMC population (“All”); four subpopulations AUS, IND, TEJ, and TRJ; and the two subspecies indica (INDAUS = IND + AUS + IND/AUS admixtures) and japonica (JAP = TEJ + TRJ + TEJ/TRJ admixtures). For each filtered panel, Tassel V.5 (Bradbury et al., 2007) was used to generate a centered identity by state (IBS) kinship matrix and to conduct principal component analysis (PCA). The number of principal components (PCs) included in the mixed linear model analysis of each panel was adjusted in the following manner to account for population structure within each panel. For GWAS of the full RMC, the first three PCs were used as covariates in the mixed linear model; two PCs were used with the INDAUS subspecies; 1 PC with the JAP subspecies and the IND subpopulation; and no PCs used for TEJ, TRJ, and AUS. Mixed linear models were performed in Tassel version 5 (Zhang et al., 2010) with the variance components estimated for each marker and no compression options.

As described in Huggins et al. (2019), a script was used that identified associated chromosome regions from individual SNPs or groups of physically linked SNPs. Chromosome regions included 50 kb in both directions around each individual significant SNP and were extended to include nearby significant SNPs occurring within 200 kb. A “Peak SNP” was designated for each region, which corresponded to the SNP with the most significant p-value. The observed frequencies of the alleles present at the peak SNP as well as the allele effect value were outputted (Huggins et al., 2019). Quantile–quantile (QQ) and Manhattan plots were created using the R package qqman (Turner, 2014). SNPs with −log10(p) > 5 were considered significant, and increased stringency was applied during data interpretation as needed, such as when only one SNP in a region met the log10(p) = 5 threshold, or regions where all associated SNPs had rare alleles (in <6 accessions), in which cases QTL were not claimed in these results. Running GWA on multiple panels provides more opportunity to identify chromosomal regions containing QTL but also prevents identification of a single “most strongly associated” SNP per QTL since it is common for SNPs to not be polymorphic in all panels, causing panels to identify QTL via association with different SNPs. Following the long-held precedent of claiming the fewest number of genes fitting a genetic model, we declared a single QTL region rather than multiple QTL per trait when different panels identified different but closely linked SNPs, merging linked SNPs into a single QTL (or not) per the following criteria. When considering SNP peaks across panels, a single QTL was declared if the string of SNPs significantly associated with a trait in one or more panels and/or study environments did not have a gap >800 kb and the additive effects were consistent across the SNPs and panels (example, qStHD1-1 in Supplementary Table S2). A chromosomal region containing multiple associated SNPs was declared as two QTL when a gap >800 kb existed between associated SNPs (e.g., qStHD1-1 to qStHD1-2), or if there were no SNP gaps >800 kb, but the additive effect of the predominant allele changed from positive to negative within a string of associated SNPs (e.g., qSi10-2 to qSi10-3). QTL for different traits were considered to overlap if the SNPs significantly associated with the two traits were interspersed with each other or if the ends of the QTL regions were ≤500 kb apart.

The multiple SNPs identified by GWA analyses for each QTL region were further evaluated using pivot tables in Microsoft Excel to characterize their allele effects on additional traits, e.g., to evaluate the effect of Si QTL on StHD or grain-As. Most of the QTL identified in a subpopulation panel were also significant in the GWA of all RMC, allowing Table 2, discussion text, and the pivot analysis to focus on SNP peaks as identified in the GWAS of “All” RMC. Calcium was unusual among the traits with several QTL significant in one subpopulation and not also significant in the related subspecies or “All” panels.

Annotated functions of genes within the QTL regions were evaluated to identify candidate genes underlying the QTL for traits central to this study. Gene annotations from multiple public databases were merged into one file for this effort (Huggins et al., 2019) and included the candidate genes for general biotic and abiotic stress response genes identified in Cohen and Leach (2019), merged with all gene annotations in the Os-Nipponbare-Reference-IRGSP-1.0 assembly (Ouyang et al., 2006; Kawahara et al., 2013), the Rice Annotation Project (RAP1; http://rapdb.dna.affrc.go.jp/, accessed 26 Sept. 2019) (Sakai et al., 2013), and Oryzabase (OrzbaseGeneListEn_20190424010057; https://shigen.nig.ac.jp/rice/oryzabase/download/gene, accessed 26 Sept. 2019), and in Cohen and Leach (2019), which examined general biotic and abiotic stress response genes in rice.

The StHD ratings for the susceptible Cocodrie checks averaged 6.0 ± 0.9 in 2015 and 5.3 ± 0.7 in 2016 on the 0–9 rating scale (Yan et al., 2005). Zhe733 had notably less severe StHD symptoms with rating averages of 2.3 ± 0.5 in 2015, 2.1 ± 0.3 in 2016, and plot ratings ranging 2 to 3 in both years. Check plot ratings verified the production of MSMA-induced symptoms across the experimental areas in both years. Among the RMC, the average StHD 2-year LSmean in MSMA-treated plots was 6.7 ± 1.3 and ranged from 3.4 to 9 (Figure 1), with none of the accessions being as StHD resistant as Zhe733, the resistant check. The majority of the RMC (101 of 156 accessions) were as or more susceptible than Cocodrie, with LSmeans >6.0. None of the plots planted in the Native soil treatment received StHD ratings >3.

Summaries of all traits observed across the entire RMC (“All”) as well as divided among the subspecies (indica and japonica) and among the subpopulations (ARO, AUS, IND, TEJ, TRJ) are presented in Figure 1. Table 1 presents Pearson correlation coefficients between the trait BLUPs for the whole RMC (“All”) and for the indica and japonica subspecies. Subpopulation correlations are not presented or discussed further because their smaller panels identified relatively few significant correlations and their correlation patterns were similar to those of the indica and japonica subspecies panels. On average, the As-induced stress in the MSMA plots delayed DHD by 2 days and decreased height by 13 cm compared to the Native soil plots, but both DHD and PHT were highly correlated between the MSMA and Native soil treatments (r = 0.90 and 0.79, respectively, Table 1). While hull-Si was lower in Beaumont, TX (avg. 188 mg kg⁻¹), compared to Stuttgart, AR (avg. 215 mg kg⁻¹), the hull-Si calculated across the two locations was highly correlated with the Si data from each individual location (r = 0.86 for both Beaumont, TX, and Stuttgart, AR, Table 1). A positive correlation was observed between PHT and DHD under both Native and MSMA conditions among all panels in this RMC study.

Among the subpopulations, the TEJ were notably more resistant to StHD, earlier, and shorter (Figure 1). The TEJ were also lower in hull-Si, grain-As, grain-S, and grain-Ca, but higher in grain-Cu. This exemplifies the need to account for population structure when conducting GWA analyses across multiple subpopulations, accomplished in this study by using PCs as covariates.

StHD and grain-As were always positively correlated, although the correlations were significant in the entire population (r = 0.23, p = 0.011) and among the indica accessions (r = 0.36, p = 0.003), but not significant among the japonica accessions (r = 0.50, p = 0.10) (Table 1). Contrary to our working hypothesis, Si, P, S, and Ca were not associated with either StHD or grain-As. Copper was negatively correlated with StHD, as predicted if Cu reduces StHD severity by increasing ROS scavenging, but Cu was also negatively correlated with grain-As. In fact, the strongest correlations with StHD are PHT from Native soil (PHTnt) (r = 0.39, p < 0.0001) and Cu (r = −0.31, p < 0.001), followed by grain-As (r = 0.23, p < 0.05), and the traits most strongly correlated with grain-As are Cu (r = −0.41, p < 0.0001) and PHTnt (r = 0.36, p < 0.0001), suggesting that correlations with PHT might be underlying the correlations observed between Cu with either StHD or grain-As.

BN analysis is useful for predicting the likelihood that any one of several possible causes was the contributing factor to an observed event and presents the modeled predictions graphically using uni- or bidirectional arrows in the resulting DAG. Because of the high correlations (r > 0.75) observed between PHT and DHD in both Native and MSMA soils, and between hull-Si calculated across both locations with the hull-Si per individual location (r > 0.80), the BN analysis included only one version of these traits, namely, DHDnt, PHTnt, and hull-Si combined by BLUP across locations. Figure 2 presents the DAG determined by BN to best fit the phenotypic data from “All” RMC. The relationship between DHD and PHT was considered bidirectional (no arrow points on the connecting line), while all other relationships included in the BN model are unidirectional cause–effect relationships. The DAG indicates that the Pearson correlations observed between StHD, grain-As, and grain-Cu are likely caused by all three traits being affected by DHD, either directly or indirectly through PHT. The DAG also indicates no direct relationships between Si, S, Ca, and Cu with neither StHD nor grain-As. The DAG indicates a directional effect of As on P, a relationship which was positive but not significant in the Pearson correlations.

Figure 3 presents the 195 GWA-QTL identified among the RMC for nine traits. Tables 2–4 present the QTL for Si (33 QTL), StHD (23), and As (15), respectively, while Supplementary Table S3 presents the QTL identified for the remaining As-related traits, P (11), S (18), Ca (19), Cu (9), DHD (39), and PHT (28). Manhattan and QQ plots are in Supplementary Figure S1. Because identification of a QTL in more than one population or environment not only validates that QTL but also demonstrates the reliability of the methods used to identify QTL in a particular study, when an RMC QTL was co-located with a previously reported QTL for the same or similar trait, this is noted in Tables 2–4 and Supplementary Table S3 by citing the previous QTL study or studies. Because the present study reevaluated the same data previously used to identify 12 hull-Si QTL (Bryant et al., 2011), we expected finding many of the same loci, and eight of the 33 Si RMC GWA-QTL did coincide with a QTL identified by Bryant et al. (2011) (Table 2). Seven additional Si RMC QTL coincided with QTL for hull, stem, or root Si concentrations reported by Dai et al. (2005) or Wu et al. (2006), and one encompassed the Lsi1 gene. These 16 validated Si QTL are now mapped more precisely in the present high-density GWA study. Of the four known Lsi transporter genes, only Lis1 on chr2 coincides with a RMC Si QTL (Figure 3; Table 2). Among the 23 StHD GWA-QTL (Table 3), seven were close (≤1.2 Mb distant) from a QTL previously reported for StHD (Agrama and Yan 2009; Pan et al., 2012; Li et al., 2016; Li et al., 2017), arsenite toxicity (Murugaiyan et al., 2019), or the chemical mimic, germanium toxicity (Talukdar et al., 2015), increasing confidence also in the 15 novel StHD RMC GWA-QTL. The most commonly identified StHD QTL was qStHD8-2, which co-located with QTL previously identified in four different populations (Pan et al., 2012; Li et al., 2016; Murugaiyan et al., 2019). When the RMC QTL for grain-As were compared with previously reported As QTL determined in various biparental mapping or GWA populations (Zhang et al., 2008; Norton et al., 2012a; Norton et al., 2014; Zhang et al., 2014; Norton et al., 2019; Liu et al., 2021), nine of the 15 RMC grain-As GWA-QTL coincided with a previously reported grain-As locus (Table 4). Four of the 11 RMC QTL for P, six of 18 for S, one of 19 for Ca, and 2 of 9 for Cu were validated by known transporter genes or previous QTL studies (Zhang et al., 2014; Liu et al., 2021) (Supplementary Table S3). Thirty-two percent (41/128) of the StHD and elemental QTL were validated by other studies, leaving 84 of the RMC QTL for As, Si, StHD, P, S, Ca, and Cu as novel. With the literature containing more than 200 reports of QTL affecting rice DHD and PHT, encompassing much of the rice genome, the present DHD and PHT GWA-QTL were compared only with genes confirmed to affect DHD or PHT. Seven of the 39 DHD QTL found among the RMC and 13 of the 28 PHT QTL encompassed a known DHD or PHT gene (Supplementary Table S3).

Among the 23 StHD QTL and 15 As QTL, four were co-located (Tables 3, 4; Figure 3). Eleven of the StHD QTL and four of the As QTL coincided with Si. Although less common, P, S, Cu, and Ca QTL were found co-located with some StHD and As QTL. In all four instances of coincidence between StHD and As QTL, a pivot table evaluation of allele effects showed that the allele that decreased StHD at qStHD1-3, qStHD1-4, and qStHD8-1 also decreased As and Si. Similarly, the qStHD9-1 allele that decreased StHD severity also decreased both As and P. Thus, we did find some QTL overlaps lending support to the hypothesis that decreased As-uptake, as evidenced by decreased grain-As, hull-Si, or grain-P concentrations, results in decreased StHD as well. Because the DAG (Figure 2) suggested that overlap of StHD and As QTL might instead be caused by mutual dependence on DHD, it must be noted that all three instances of coincidence between StHD, As, and Si QTL also coincide with DHD. The one instance of overlap between P, StHD, and As (at qStHD9-1) does appear to be independent of DHD and PHT, as well as all other traits in this study (Table 3; Figure 3).

In agreement with the DAG (Figure 2) showing a strong influence of DHD on StHD via PHT, 10 of the 23 StHD QTL co-located with a PHT QTL, 8 with DHD, for a total of 12 that coincide with DHD, PHT, or both (Figure 3; Table 3). The DAG did not indicate a close relationship between StHD and Si, nor between Si and DHD. Even so, co-location among StHD and Si QTL was high, with 12 of the 23 StHD loci overlapping or being within 1 Mb of a Si QTL. This suggests that Si concentrations may in fact affect StHD severity in a manner not detected in this data by BN. For example, Si might affect StHD but less strongly than DHD and PHT. Similarly, four As QTL (qAs2, qAs3, qAs7-1, and qAs8-1), and three StHD QTL (qStHD1-1, qStHD3, and qStHD8-1) were on or near S QTL (Tables 3, 4; Figure 3), suggesting that increased As-chelation and/or increased ROS scavenging by increased S might be affecting grain-As or StHD. Overlaps between StHD and Ca or Cu QTL were notably uncommon and always involved another trait such as DHD or Si (qStHD1-4 and qStHD7 with Cu QTL, and qStHD5-1, qStHD6-1, and qStHD10 with Ca QTL), indicating that increased ROS scavenging from increased Ca or Cu is not contributing significantly to the observed variance for StHD resistance.

There are four instances of co-location between As and P QTL (chr 2, 6, 7, and 9; Table 4 and Figure 3), with pivot table analyses showing the allele that increased As also increased P. The instance of overlap on chr9 was discussed above as independent from DHD, but the other three do coincide with a DHD QTL. While the DAG does not show a direct connection between DHD and P, it does show a DHD effect on P channeling through As (Figure 2), and DHD-P correlations were significant and positive among the entire RMC and indica subspecies, but the DHD-P correlations were near zero among japonica accessions (Table 1).

To explore candidate genes in the 33 Si RMC GWA-QTL (Table 2), we searched annotated genes to find those involved with transmembrane transport of silicon transport or metals and metalloids. The four known Low silicon transporters (Lsi1, Lsi2, Lsi3, and Lsi6) are of two types, with Lsi2 and Lsi3 being silicon efflux transporters (SIETs), and Lsi1 and Lsi6 being aquaporins of the Nodulin-26 like intrinsic protein (NIP) type. Only one of the four Lsi genes, specifically Lsi1, coincided with an RMC Si QTL (Table 2; Figure 3). NIPs are a subclass of membrane intrinsic transporters (MIP) which include also plasma membrane intrinsic proteins (PIPs) and tonoplast intrinsic proteins (TIPs). Overall, across the 33 Si RMC GWA-QTL regions (Table 2), we identified as candidate genes two SIETs, seven NIPs, three PIPs, five (TIPs), and 13 multidrug and toxic compound extrusion (MATE) genes. All these genes were located within 1 Mb of 17 of the 33 Si RMC GWA-QTL.

To explore candidate genes in the 15 As RMC GWA-QTL (Table 4), we again looked in QTL regions for components of silicon and heavy metal uptake or transport and added to the search genes for signaling, and/or tolerance mechanisms (Abedi and Mojiri 2020) as well as chelation and detoxifying enzymes reported in the literature as affecting metalloid metabolism (Chen et al., 2017; Abedi and Mojiri, 2020). Overall, across the 15 As RMC GWA-QTL regions, we identified three NIPs, one SIET/Lsi3, four MATE genes, two ABC family proteins, one PIP, one TIP, one phosphate transporter, two transmembrane proteins (i.e., transmembrane receptor protein and EnS-10, endosperm-specific gene with metal ion transmembrane transporter activity), three metabolizing enzymes (phosphatase, arsenate reductase (HAC1;2), and one similar to phytochelatin synthetase), and one heat shock protein as noted in Table 4. These genes were located within 1 Mb of 10 of the 15 As RMC GWA-QTL.

The candidate gene identified by the overlapping qSi4-2/qAs4 QTL, PIP2;6, was initially reported to be involved in influx and efflux of boron transport (Mosa et al., 2016), but the heterologous expression of OsPIP2;6 in Xenopus laevis oocytes was also found to increase the uptake of arsenite (Mosa et al., 2012). Several NIP proteins have been shown to transport As, including the well-known Lsi1 and Lsi6 aquaporin transporters, NIP1;1 (Kamiya et al., 2009), and Arabidopsis NIP3;1 (Xu et al., 2015). The NIP3;2 protein (candidate for qSi8-1) has also been shown permeable to arsenite (Bienert et al., 2008; Ma et al., 2008). However, while this Si QTL overlapped with qStHD8-1, it did not overlap with the nearby qAs8. MATE genes play an important role in cellular detoxification processes. Using a transcriptomic study, Seth et al. (2020) identified nine MATE genes that were upregulated upon exposure to As. The constitutive expression of OsMATE2 in transgenic tobacco plants decreased root-to-shoot As transfer, and transgenic rice plants wherein RNAi was used to silence OsMATE2 produced grains with less grain-As than wild-type rice plants (Das et al., 2018). While MATE2 was not among the grain-As candidate genes, it is a candidate gene for qStHD5-3.

Because StHD has been attributed to As toxicity (Tang et al., 2020a), to explore candidate genes underlying the 23 StHD RMC GWA-QTL (Table 4), we included in the search the same components of Si and heavy metal uptake, transport, signaling, and/or tolerance mechanisms, chelating and detoxifying enzymes included in the As candidate gene search, and added also stress-protective agents (e.g., redox homeostasis and antioxidant defense). Across the 23 StHD RMC GWA-QTL regions (Table 4), we identified 10 ABC family transporter proteins, six NIPs, two MATEs, one Si transporter, two membrane proteins, eight enzymes (i.e., vacuolar H + -ATPase subunit C antioxidant protein 1 (ATX1), glycolate oxidase 6 (GLO6), flavonol synthase, oxidoreductase, glutathione S transferase 5 (GST5)), two purple acid phosphatases (PAP3B, PAP20A), and one plant cadmium resistance (PCR1) gene within 1 Mb of 12 of the 23 StHD RMC GWA-QTL. Vacuolar ATPase is a proton pump protein whose role in acidification of the vacuolar compartment activates the uptake and release of ions and metabolites (Sze et al., 1992; Lüttge and Ratajczak, 1997). Furthermore, ATX1 is known to interact with heavy metal P1B-ATPases, suggesting its role in delivering Cu to heavy metal P1B-ATPases for Cu trafficking and distribution in order to maintain Cu homeostasis in rice (Zhang et al., 2018). RLK1 encodes a plasma membrane-localized protein that acts upstream of mitogen-activated protein kinase (MPK) cascades and positively regulates defense-related MPKs and WRKY transcription factors to respond differentially to external stimuli (Meng and Zhang, 2013; Macho and Zipfel, 2014). The ubiquitin-proteasome system (UPS) is one of major protein regulation pathways that enhance the adaptation and survival of plants under various environmental stresses such as toxic metalloid exposure as well as other environmental stresses (drought, salinity, and cold). Recent study showed active E3 ligase activity of OsAIR2 through in vitro ubiquitination assay and further examined that overexpression lines of OsAIR2 in Arabidopsis improved the seed germination and increased the root length under arsenate stress conditions, suggesting its role as a positive regulator of As stress tolerance (Hwang et al., 2017). Germin-like protein (GLP) is a plant glycoprotein associated with the plant cell wall, and its various proposed roles in plant development and defense are known (Dunwell et al., 2008; Caliskan, 2009).

The region of Chr8 encompassing qAS8 and qStHD8-1 (Chr8 0.45-0.82 MB) contains a series of BTB-domain-containing proteins (MB17&18, MBTB16, 17, 18, 19; BTBN17) that are also of interest. A BTB domain-containing protein was reported to be involved in enhancing iron homeostasis in apple (Zhao et al., 2016). Recent study also showed that rice plants overexpressing a BTB domain protein (OsTAZ24) promoted plant resistance against various heavy metals (Shalmani et al., 2021), demonstrating the role of the BTB-domain-containing protein in heavy metal homeostasis, making these genes candidates for both the grain-As and StHD QTL. We therefore list this cluster of BTB-domain proteins as candidate genes for both qAs8 (Table 4) and qStHD8-1 (Table 3).

There were two RMC GWA-QTL regions that were associated with Si along with As and StHD (qAs1-2/qStHD1-3/qSi1-2, and qAs1-3/qStHD1-4/qSi1-3). Interestingly, these two regions, one on Chr1 at 23.34–24.45 Mb and the other one on Chr1 at 31.72–34.22 Mb, contain eight and 11 F-box domain-containing proteins, respectively. Recent study showed that 12 F-box domain-containing proteins were located in the As QTL regions from a japonica diversity panel of 228 accessions and from 95 advanced breeding lines with japonica genetic backgrounds phenotyped for As concentration in the flag leaf as well as in dehulled grain (Frouin et al., 2019). One of the F-box proteins identified as a candidate gene in our study, OsBOX028, was also identified by Frouin et al. (2019), suggesting that F-box proteins may have a role in affecting As tolerance. F-box proteins are known to define the specific substrates of the SCF complexes that are part of the ubiquitin/26S proteasome pathway for degrading unwanted or misfolded proteins (Zheng et al., 2002; Lechner et al., 2006). Thus, these F-box gene complexes are listed as candidate genes for qStHD1-3 and qStHD1-4 QTL (Table 3), but not for the coincident As and Si QTL (Tables 2, 4).

Another interesting aspect of the QTL overlaps was a region on chr9 containing QTL for StHD and As but not for Si. A QTL for P (qP9) was also located with qStHD9-1/qAs9. When we further examined the genes in this overlapping region, we found that more than 50% of genes in this region were retro/transposon proteins, forty retrotransposon proteins and five transposon proteins out of 99 total genes in this overlapping QTL region of qAs9, qStHF9-1, and qP9 (Chr. 9, 3,356,592–3,937,895). Transposons and retrotransposons play a role in post-transcription regulation through silencing of transposable elements (TEs), which can elicit gene variation and functional changes (Gao et al., 2012). TEs are also a source of small RNAs. A role for small RNAs in regulation of P starvation has been reported (Fang et al., 2009; Hsieh et al., 2009; Chiou and Lin 2011; O’Rourke et al., 2013).

Eleven of the 33 StHD QTL were co-located with Si QTL in the RMC, only three of which were also co-located with As QTL. This suggests a benefit from Si on StHD that is not directly associated with As and shared Si–As transporters. Numerous recent studies have confirmed Si beneficial effects on a variety of plant species growing under a wide range of environmental conditions. Si not only has a role in ameliorating biotic and abiotic stress but also can alleviate nutrient deficiency and toxicities (Pavlovic et al., 2021). Most of the studies regarding Si nutrient interactions were studies conducted on Si accumulators like rice. A general model for Si interactions has been proposed based on those studies for silicon-mediated response to mineral deficiency and toxicity (Ali et al., 2020). Under nutritional stress, uptake of Si increases through Lsi transporters (Yamaji and Ma, 2011) or aquaporins (Deshmukh et al., 2013). The increase in Si induces genes encoding the uptake and translocation of minerals in roots and facilitates the increase or decrease uptake of the corresponding nutrients such as N (Wu et al., 2017) and P (Owino-Gerroh and Bascho, 2005). In shoots, Si interacts with phytohormones and amino acids, and metabolites lower oxidative stress by further modulating the antioxidant enzymes (Liang et al., 2003; Li et al., 2012; Chalmardi et al., 2014) and simultaneously maintaining or even inducing higher photosynthetic efficiency resulting in plant growth and development (Detmann et al., 2012). On the other hand, for mineral or heavy metal toxicity, Si forms a complex with some heavy metals (Cd) or minerals (Zn, Fe, Cu) in the root cell wall, therefore reducing uptake and translocation to the shoots (Liu et al., 2013; Kopittke et al., 2017; Nozawa et al., 2018; Bosnic et al., 2019). In shoots, Si accumulation further mitigates metal toxicity by sequestration of minerals (Cu, Zn) and heavy metals (As, Cd) into the vacuoles of leaf cells (Keller et al., 2015; Kopittke et al., 2017; Nozawa et al., 2018). Overall, Si has the potential of enabling the plant to react adaptively against nutritional stress and promotes tolerance.

The initial hypothesis was that both StHD and grain-As would be regulated strongly enough by the concentration of free (unbound, non-sequestered) As in the plant that they would share numerous QTL. We further hypothesized that overlap of the StHD and As QTL with Si, P, or S would provide further insight as to whether the shared StHD and As QTL were being driven by As-uptake or post-uptake metabolism. Because data existed to add Ca and Cu to the study, and because these elements have been shown to increase ROS scavenging, which could in turn reduce StHD severity, they were included as well. Correlations were positive between StHD, As, and Si, as would be expected if our first two hypotheses were true. However, finding a stronger StHD-Cu correlation than seen for StHD-As was unexpected. The BN DAG results cautioned to not overlook possible confounding effects from DHD in our data interpretations. When QTL were identified, and overlaps considered, the number (4 total) and percentage of StHD (0.12) and As (0.27) QTL that were co-located were lower than hypothesized. We also predicted that Si would impact StHD via its influence on As uptake. In contrast, we found that significantly more StHD QTL overlapped with Si QTL (11/33) than with As QTL (4/33). The results indicate that StHD and grain-As are both more complex than initially hypothesized and indicate that Si QTL may be affecting StHD response in ways not connected by As concentration, such as by increasing ROS scavenging, or improving overall plant health and nutrient balance.

Even though the hypothesis based on an StHD-As association was not well supported, the StHD and Si QTL identified among the RMC accessions in this study are the first to be identified using high-density mapping and, thus, are mapped to smaller QTL regions than previously available for these traits. Furthermore, the high proportion of QTL identified among the RMC accessions that were validated from those in literature indicates that the QTL-identification methods used were robust, increasing confidence in both the validated and novel QTL we report. The novel Si QTL, in particular, are worth deeper investigation since many of them were not explainable by currently known Si transporters and because Si did appear to have a stronger impact on StHD than could be explained by As alone. The region of chr8 that has been commonly reported as associated with StHD (qStHD8-2) also warrants further investigation. This qStHD8-2 QTL region was not associated with either As or Si but was found to contain multiple genes of interest including genes for antioxidants and redox factors which could be reducing ROS damage, suggesting that differences in ROS damage should be considered in future research on StHD resistance mechanisms. The region of chr9 containing qStHD9-1, qAs9, and qP9 which was independent from DHD would also be a genomic segment warranting further research. Other genes of interest for future study, because they were found in multiple grain-As and StHD regions, include several multidrug and toxic compound extrusion (MATE) genes, F-box genes, BTB-domain proteins, and NIPs not documented to date to transport As. The candidate genes we identified could be a basis for gene editing studies to determine the genes underlying the StHD, As, and Si QTL identified in the RMC.