The Interaction of Genotype and Environment Determines Variation in the Maize Kernel Ionome

Plants obtain soil-resident elements that support growth and metabolism from the water-flow facilitated by transpiration and active transport processes. The availability of elements in the environment interacts with the genetic capacity of organisms to modulate element uptake through plastic adaptive responses, such as homeostasis. These interactions should cause the elemental contents of plants to vary such that the effects of genetic polymorphisms will be dramatically dependent on the environment in which the plant is grown. To investigate genotype by environment interactions underlying elemental accumulation, we analyzed levels of elements in maize kernels of the Intermated B73 × Mo17 (IBM) recombinant inbred population grown in 10 different environments, spanning a total of six locations and five different years. In analyses conducted separately for each environment, we identified a total of 79 quantitative trait loci (QTL) controlling seed elemental accumulation. While a set of these QTL was found in multiple environments, the majority were specific to a single environment, suggesting the presence of genetic by environment interactions. To specifically identify and quantify QTL by environment interactions (QEIs), we implemented two methods: linear modeling with environmental covariates, and QTL analysis on trait differences between growouts. With these approaches, we found several instances of QEI, indicating that elemental profiles are highly heritable, interrelated, and responsive to the environment.

Quantitative genetics using structured recombinant inbred populations is a powerful tool for dissecting the factors underlying elemental accumulation and relationships. By breaking up linkage blocks through recombination and then fixing these new haplotypes of diverse loci into mosaic sets of lines, these populations allow similar sets of alleles to be repeatedly tested in diverse environments (Koornneef et al. 1997). A variety of quantitative statistical approaches can then be used to identify quantitative trait loci (QTL) by environment interactions (QEI).
Here, we used elemental profiling of a maize recombinant inbred population grown in multiple environments to identify QTL and QEI underlying elemental accumulation. We sought both environmental and genetic determinants by implementing single-environment QTL mapping and analyses of combined data from multiple environments. Overall, we detected 79 loci controlling elemental accumulation, many of which were environment-specific, and identified loci exhibiting significant QEI.  (87 lines). In all but three environments, NY05, NC06, and MO06, one replicate was sampled per line. In NY05, three replicates of 199 lines, two replicates of 50 lines, and one replicate of seven lines were sampled. A replicate is considered pooled ears from a row. Several ears were harvested and kernels were subsampled from pooled ears from the row. After harvesting, seeds were stored in local temperature and humidity controlled seed storage rooms. Subsequently they were shipped to the ionomics laboratory, where they were stored in temperature-controlled conditions. Because each batch of seed was treated identically, any losses in weight or increases in weight due to differing hydration should not affect the relative, weight-adjusted concentrations used for analysis. We do not expect any changes in ion composition due to storage. Supplemental Material, Table S1 includes planting dates and line numbers after outlier removal and genotype matching. After outlier removal, 199 of the 233 unique lines in the experiment were present in three or more of the 10 environments. 106 lines were present in seven or more of the environments.

Elemental profile analysis
Elemental profile analysis is conducted as a standardized pipeline in the Baxter laboratory. The methods used for elemental profile analysis are as described in Ziegler et al. (2013). Descriptions taken directly are denoted by quotation marks.
Sample preparation and digestion: Lines from the IBM population from each environment were analyzed for the concentrations of 20 elements. "Seeds were sorted into 48-well tissue culture plates, one seed per well. A weight for each individual seed was determined using a custom built weighing robot. The weighing robot holds six 48-well plates and maneuvers each well of the plates over a hole which opens onto a three-place balance. After recording the weight, each seed was deposited using pressurized air into a 16 · 110 mm borosilicate glass test tube for digestion. The weighing robot can automatically weigh 288 seeds in 1.5 hr with little user intervention." "Seeds were digested in 2.5 ml concentrated nitric acid (AR Select Grade, VWR) with internal standard added (20 ppb In, BDH Aristar Plus). Seeds were soaked at room temperature overnight, then heated to 105°for 2 hr. After cooling, the samples were diluted to 10 ml using ultrapure 18.2 MV water (UPW) from a Milli-Q system (Millipore). Samples were stirred with a custom-built stirring rod assembly, which uses plastic stirring rods to stir 60 test tubes at a time. Between uses, the stirring rod assembly was soaked in a 10% HNO 3 solution. A second dilution of 0.9 ml of the 1st dilution and 4.1 ml UPW was prepared in a second set of test tubes. After stirring, 1.2 ml of the second dilution was loaded into 96 well autosampler trays." Ion coupled plasma mass spectrometry analysis: Elemental concentrations of B, Na, Mg, Al, P, S, K, Ca, Mn, Fe, Co, Ni, Cu, Zn, As, Se, Rb, Sr, Mo, and Cd "were measured using an Elan 6000 DRC-e mass spectrometer (Perkin-Elmer SCIEX) connected to a PFA microflow nebulizer (Elemental Scientific) and Apex HF desolvator (Elemental Scientific). Samples were introduced using a SC-FAST sample introduction system and SC4-DX autosampler (Elemental Scientific) that holds six 96-well trays (576 samples)." Measurements were taken with dynamic reaction cell (DRC) collision mode off. "Before each run, the lens voltage and nebulizer gas flow rate of the ICP-MS were optimized for maximum Indium signal intensity (.25,000 counts/sec) while also maintaining low CeO+/Ce+ (,0.008) and Ba++/Ba+ (,0.1) ratios. This ensures a strong signal while also reducing the interferences caused by polyatomic and double-charged species. Before each run a calibration curve was obtained by analyzing six dilutions of a multi-element stock solution made from a mixture of single-element stock standards (Ultra Scientific). In addition, to correct for machine drift both during a single run and between runs, a control solution was run every tenth sample. The control solution is a bulk mixture of the remaining sample from the second dilution. Using bulked samples ensured that our controls were perfectly matrix matched and contained the same elemental concentrations as our samples, so that any drift due to the sample matrix would be reflected in drift in our controls. The same control mixture was used for every ICP-MS run in the project so that run-to-run variation could be corrected. A run of 576 samples took 33 hr with no user intervention. The time required for cleaning of the instrument and sample tubes as well as the digestions and transfers necessary to set up the run limit the throughput to three 576 sample runs per week."

Computational analysis
Drift correction and analytical outlier removal: Analytical outliers were removed from single-seed measurements using a method described by Davies and Gather (1993). Briefly, values were considered an outlier and removed from further analysis if the median absolute deviation (MAD), calculated based on the line and location where the seed was grown, was .6.2.
Normalization for seed weight by simply dividing each seed's solution concentration by sample weight resulted in a bias where smaller seeds often exhibited a higher apparent elemental concentration, especially for elements whose concentration is at or near the method detection limit. This bias is likely either a result of contamination during sample processing, a systematic over or under reporting of elemental concentrations by the ICP-MS, or a violation of the underlying assumption that elemental concentration in seeds scales linearly with seed weight. Instead, we developed a method taking residuals from the following linear model: where Y is the nonweight normalized measure of elemental concentration for each seed after digestion, b 0 is the population mean, X 1 is the seed weight, X 2 is the analytical experiment the seed was run in (to further correct for run-to-run variation between analytical experiments), and e is the residual (error) term. The residuals in this linear model represent how far each data point departs from our assumption that analyte concentration will scale linearly with sample weight. If all samples have the same analyte concentration then the linear model will be able to perfectly predict analyte concentration from weight and the residuals will all equal zero. However, if a sample has a higher or lower concentration of an analyte than the general population being measured, then it will have a residual whose value represents the estimated concentration difference from the population mean. For this reason, we have termed this value the estimated concentration difference from the mean (ECDM).
Heritability calculation: Broad-sense heritability was calculated for seed weight and 20 elements across environments and within three environments for which we had substantial replicate data. To estimate the broad-sense heritability across 10 environments, the total phenotypic variance was partitioned into genetic and environmental variance, with the broad-sense heritability being the fraction of phenotypic variance that is genetic. This was done using an unbalanced, type II ANOVA in order to account for the unbalanced common line combinations across environments. Two models were fit using the lmfit function in R. The first model included genetic variance as the first term, and environmental variance as the second. The second model had the opposite form. The variances for genetic or environmental components were obtained using the ANOVA function on the model in which that component was the second term. Broad-sense heritability was calculated by dividing the genetic variance by the total (genetic plus environmental plus residuals) variance. Heritability was calculated within environments for NY05, NC06, and MO06. Data with outliers designated as NA was used for each environment. For each element within an environment, lines with NA were removed, and lines with only one replicate were removed, leaving only lines with two or more replicates. The heritability was then calculated for seed weight and each element using the lmfit and anova functions to obtain the variances for the genetic component and the residuals. Broad-sense heritability was calculated as the proportion of total variance (genetic plus residuals) explained by the genetic component.

QTL mapping and elemental traits:
The R package R/qtl was used for QTL mapping. For each of the 10 environments, elemental trait line averages and genotypes for all lines, 4217 biallelic single nucleotide polymorphisms (SNPs) distributed across all 10 maize chromosomes, were formatted into an R/qtl cross object. The stepwiseqtl function was used to implement the stepwise method of QTL model selection for 21 phenotypes (seed weight, 20 elements). The maximum number of QTL allowed for each trait was set at 10 and the penalty for addition of QTL was set as the 95th percentile LOD score from 1000 scanone permutations, with imputation as the selected model for scanone. A solely additive model was used; epistatic and interaction effects were not considered and thus heavy and light interaction penalties were set at 0. QTL positions were optimized using refineqtl, which considers each QTL one at a time, in random order, iteratively scanning in order to move the QTL to the highest likelihood position. QTL models for each trait in each environment were obtained using this procedure. QTL within 5 cM of each other were designated as the same QTL.
QTL by environment analysis: Linear model comparison: Linear modeling was used to determine instances and strength of QEI using all data from two different years within three locations (FL, IN, and NY). The specific growouts analyzed together were FL05, FL06, IN09, IN10, NY05, and NY12. FL, IN, and NY were then used as covariates in QTL analysis. Two QTL models, one with location as an additive and interactive covariate and one with location as only an additive covariate, were fit for each phenotype (sample weight, 20 elements) using the scanone function in R/qtl, where y i is the phenotype of individual i, g i is the genotype of individual i, and x i is the location of individual i. B g and B x are additive effects of genotype and environment, respectively, and g is the effect of genotype by environmental interaction. LOD scores for each marker using model (2) were subtracted from LOD scores for each marker using model (1) to the isolate genetic by location effect. QTL by location interaction was determined as QTL with a significant LOD score after subtraction. The significance threshold was calculated from 1000 permutations of the three-step procedure (fitting the two models and then subtracting LOD scores), and taking the 95th percentile of the highest LOD score.
Mapping on within-location differences: QTL were mapped on phenotypic differences between common lines grown over 2 yr at a single location. This procedure was used to compare FL05 and FL06, IN09 and IN10, and NY05 and NY12 by calculating the differences for each trait value between common lines in location pairs (FL05-FL06, IN09-IN10, and NY05-NY12) and using these differences for analysis using the previously described stepwiseqtl mapping and permutation procedure.

Data availability
All data and scripts are available on Ionomics Hub (iHUB) in the Maize Database at www.ionomicshub.org.

Genetic regulation of elemental traits
The data used for this study is comprised of 20 elements measured in the seeds from the Zea mays L. intermated B73 · Mo17 recombinant inbred line (IBM) population grown in 10 different location/year settings. The IBM population is a widely studied maize population of 302 intermated recombinant inbred lines, each of which have been genotyped with a set of 4217 biallelic SNP genetic markers (Lee et al. 2002). The four rounds of intermating and subsequent inbreeding generated increased recombination and a longer genetic map for the IBM than for typical biparental recombinant inbred line populations. The number of individuals, marker density, and greater recombination facilitates more precise QTL localization than a standard RIL population ( (Table S1). While very few of the 233 unique IBM lines in the experiment were grown in all environments, 106 of the 233 lines were grown in seven or more environments and 199 were grown in three or more environments. Within each growout, all samples were treated identically: seeds from all environments were stored in temperature and humidity controlled storage rooms after harvest and then shipped to the ionomics lab. We do not expect any change in ion composition from storage within a growout; however, we cannot rule out that some of the differences between growouts might be due to slightly different moisture content. These differences are not likely to account for the genetic by environment interactions we observe as they should have similar effects on all lines. Single seeds were profiled for the quantities of 20 elements using ICP-MS. These measurements were normalized to seed weight and technical sources of variation using a linear model, with the resulting values used as the elemental traits for all analyses (Shakoor et al. 2016). After outlier removal, seed element phenotypes were derived by averaging line replicates (kernels subsampled out of pooled ears from a row) within an environment. Variation in the elemental traits was affected by both environment and genotype. Elemental traits generally exhibited lower heritability among genotypes grown across multiple environments than among genotype replicates within a single environment (Table 1). The broadsense heritability (H 2 ) of seed weight, 15 of 21 elements in NY05, 13 of 21 elements in NC06, and 13 of 21 elements in MO06 exceeded 0.60. Elements exhibiting low heritability within environments corresponded to the elements that are prone to analytical artifacts or present near the limits of detection by our methods, such as B, Al, and As. Seven elements had a broad sense heritability of at least 0.45 in a single environment (NY05, NC06, and NY06) but ,0.1 across all environments. This decrease in heritability across the experiment, which was particularly striking for Mg, P, S, and Ni, is consistent with strong genotype by environment interactions governing the accumulation of these elements.
A stepwise algorithm, implemented via stepwiseqtl in the R package R/qtl (Broman and Speed 2002), was used to map QTL for seed weight and 20 seed elemental phenotypes. The stepwise algorithm iterates through the genome and tests for significant allelic effects of each marker on a phenotype. Forward and backward regression generates the final genome-wide QTL models for each trait. This QTL mapping procedure on 21 traits was completed as a separate analysis for each subset of lines from the IBM populations grown in each of the 10 environments. For the sake of completeness and to comprehensively investigate all of the traits to which we had access, all elemental traits in each environment were tested, even in cases where heritability for a given element was low in an environment. QTL significance were determined using the 95th percentile threshold from 1000 scanone permutations as a penalty score for adding QTL to the stepwise model (Churchill and Doerge 1994). We examined the relationship between the heritability of an element in a given environment and number of QTL identified in that environment ( Figure S1). As expected, elements with very low heritability had few to no QTL, while larger numbers of QTL were identified for higher heritability elements.
The environmental dependence on QTL detection was first estimated by identifying QTL common to multiple environments. If QTL detected in two or more growouts affected the same element and localized within 5 cM of each other, they were considered to be the same locus. Across the 10 environments, a total of 79 QTL was identified for seed weight, and 18 of the 20 elemental traits tested (none for Al or Co) (Figure 1, B and C). Of these QTL, 63 were detected in a single n   (IN09 and IN10), and New York (NY05 and NY12) growouts. All lines within each location were included, with covariates designated based on location. QTL mapping output using model with location as an additive covariate is shown as dotted line. QTL mapping output from model with location as both an additive and interactive covariate is shown as dashed line. Subtracted LOD score profile from the two models (QTL by location interactive effect only) is shown as solid line. Horizontal line corresponds to significance threshold for QTL by location interaction effect, derived from 1000 iterations of the three step procedure using randomized data: scanone QTL mapping with the additive model, scanone QTL mapping with the additive and interactive model, and subtraction of the two models.
environment and 16 were detected in multiple environments. The 16 QTL found in multiple environments included QTL detected in nearly all of the environments and QTL detected in only two. One QTL for Mo accumulation, on chromosome 1 in the genetic region containing the maize ortholog of the Arabidopsis molybdenum transporter MOT1 (Baxter et al. 2008), was found in nine environments ( Figure 1A). Another QTL affecting Cd accumulation, on chromosome 2 and without a clear candidate gene, was found in eight environments. Other QTL were present only in a smaller set of environments, such as the QTL for Ni accumulation on chromosome 9, which was found in five environments ( Figure 1D). The strength of association and percent variance explained showed strong differences between environments even for these QTL that were detected in multiple environments (Table S2).
As seen in the full-genome view of all QTL colored by environment ( Figure 1C), there is a high incidence of QTL found in single locations. Three hypotheses could explain the large proportion of QTL found only in a single location: (1) strong QTL by environment interaction effects; (2) false positive detection of a QTL in an individual location; and (3) false negative assessment of QTL absence due to genetic action, but statistical assessment below the permutation threshold in other environments. To reduce the risk of false positives in a single environment's QTL set, the significance threshold was raised to the 99th percentile, where 31 of the 63 environment-specific QTL remained significant. Despite the large number of trait/environment combinations tested (20 traits in 10 environments), the number of QTL detected is much larger than the null expectation derived from a Bonferroni correction: 10 QTL (95th percentile threshold) and two QTL (99th percentile threshold). To account for false negatives, we scanned for QTL using a more permissive 75th percentile cutoff. Of the 63 singleenvironment QTL, only nine had QTL in other environments by this more permissive threshold. Thus, the majority of the 63 singleenvironment QTL most likely result from environmentally contingent genetic effects on the ionome.

QTL by environment interactions
That QTL detection was so strongly affected by environment suggested the effects of allelic variation on element concentration were heavily dependent on environmental variables. These results, however, did not specifically test for QEI. Comparison between environments with our data are additionally complicated because different subsamples of the IBM population were grown at these multiple locations and years. While there are many different approaches to identifying QEI described in the literature (summarized in El-Soda et al. 2014), we focused on two previously implemented methods. The first considered location (but not year) by comparing the goodness-of-fit for linear models with and without an interactive covariate (Nichols et al. 2007;Leinonen et al. 2013;Bhatia et al. 2014). The second method takes advantage of the ability to grow the same RILs in multiple years. Trait values measured in the same IBM line for the same element at the same site but in different years were subtracted from each other and the difference between years was assigned as the trait value for that RIL genotype for QTL detection (Ungerer et al. 2003;Tétard-Jones et al. 2011).

Linear model estimation of QTL by location effects:
The most common approach to analyze QEI is to fit a linear model with environment as both a cofactor and an interactive covariate and compare results to a model with environment as an additive covariate (Des Marais et al. 2013). This method is most amenable when data are available for the same lines grown in every environment, which was not the case across all of our dataset. Data from the three locations with two replicate years each (FL, IN, NY) were analyzed to reduce the number of covariates and increase the power to detect variation from the environment. The data for both years in each location were combined (FL05 and FL06; IN09 and IN10; and NY05, NY06, and NY12), designating covariates based on location.
Two linear QTL models introduced in the Materials and Methods section were fit to the combined data using the FL, IN, and NY locations as covariates. These models reflect the dependence of phenotype on genotype, environment, and genotype-by-environment interactions. Stating these models again: The first equation fit (1) is the full model considering the phenotype of individual i (y i ) as controlled by genotype (g i ), location (x i ), and genotype by location interaction (g i x i ), while the reduced model (2) estimates phenotype without considering genotype by location interaction, using genotype and location as purely additive factors. B g and B x represent the additive effects of genotype and environment, respectively, while g represents the effect of the genotype by environment interaction. By using likelihood ratio tests on full and reduced models, we can test the hypothesis that genotype by environment interactions significantly improve the fit of the model to the data and estimate the effects of genetic by environment interactions. The program R/qtl was used to fit QTL using both the full and reduced models for sample weight and 20 elements, with three locations encoded as covariates in the environment term. For each marker, LOD scores resulting from the reduced QTL model were subtracted from LOD scores determined by the full model, leaving a LOD score for each marker representing solely the significance of the genetic by location n component. The significance threshold for the subtracted LOD scores was calculated by using 1000 permutations of the three step procedure (fitting the two models with randomized data and then subtracting LOD scores). Even with this underpowered dataset, 10 QTL by location interactions exceeded the threshold (a = 0.05, Table 2). Interactions between QTL and location are likely to be due to a combination of soil and weather differences across different locations. In the case of Ni, our initial single-element QTL mapping conducted separately on data from each environment identified differences in QTL presence or strength between FL, IN, and NY locations for a QTL located at the beginning of chromosome 9 (Figure 2). This QTL corresponds to a locus found to have a significant QTL by location effect (Table 2). Remarkably, all elemental QTL by location interactions detected by this approach affected trace element accumulation. These elements are both low in concentration in the grain and often variable among soils (White and Zasoski 1999). Cd, an element for which we found significant QEI, has detrimental effects on both human and plant health (Godt et al. 2006) and is toxic in food at levels as low as 0.05 ppm. (United States Department of Agriculture 2014). The locus with the strongest QEI for Cd does not follow location averages of Cd content in the grain (Table S3) and therefore is unlikely to be affected by crossing a detection threshold driven by higher Cd in the soils at those locations. The lack of direct correlation between QTL significance and grain content also occurs for the loci with strong by-location effects for Mo and Ni. This demonstrates that reduced cadmium or enhanced micronutrient contents in grain require plant breeding selections that consider complex genetic by environment interactions rather than genotypes assessed in a single soil environment.
QTL for trait differences within location: The previous method identified genotypes with interactions with location but not with year.
Year to year variation will also have effects due to differences in rainfall, temperature and management practices. To examine variation that occurs within a location over different years, we examined intralocation QEI in the three previously used locations with 2 yr samples (FL05 and FL06,IN09 and IN10,and NY05 and NY12). QTL were mapped using the stepwise algorithm on trait differences for sample weight and 20 elements between common lines among the two different years from a location. This approach identified loci affecting phenotypic differences between the same lines grown on the same farm but in different years. Six QTL were found for FL05-FL06 differences, one QTL for IN09-IN10 differences, and two QTL for NY05-NY12 differences (Table 3). These trait-difference QTL included loci identified in our single element/environment QTL experiment where a locus was present in one year, but not the other, or the QTL was found in both years with differing strength (Figure 3). Six of the difference QTL were detected at loci where no QTL were detected when the years were mapped separately, revealing novel gene by environment interactions not obvious from the single year data. These significant effects of year-to-year environmental variation within the same location indicated that factors beyond location are both influencing the ionome and determining the consequences of genetic variation. Figure 3 Comparison of QTL mapped on traits in single environments and trait differences between environments. Examples from stepwise QTL mapping on trait differences between two years at one location, calculated between IBM lines common to both years. Scanone QTL mapping output is also plotted for the same trait from each year separately. LOD score is shown on the y-axis and centimorgans posi-

DISCUSSION
The results described here demonstrate that the concentrations of elements in the kernels of maize are strongly affected by the interaction of genetics with growth environment. The majority of elements exhibited higher heritability within each environment and a dramatic drop in heritability across multiple environments. Combined with the presence of a large number of single-environment QTL, these data support the hypothesis that environment has a significant impact on genetic factors influencing the ionome. By changing the stringency of the statistical tests, we are able to discount the likelihood that that these single environment QTL are the result of a large number of false positives or false negatives. The structure of our data, with few lines measured across all locations and years, limited our ability to test for direct QTL by environment interactions. As a result, we have likely underestimated the extent of QEI. Future studies with uniform lines across environments will allow for inclusion of data from all environments and lines and increase power to detect additional genetic by environment interactions. Nevertheless, we were able identify QEI over different locations and QEI at a single location over different years. We identified a strong nickel QTL on chromosome 9 that was found in Indiana and New York with single-environment QTL mapping, but not in Florida. This same locus also was found to be a significant location-interacting QTL when using a model that included Indiana, New York, and Florida as covariates. One possible cause for this, and other location specific QTL, might be differences in element availability between local soil environments. Interestingly, the presence/absence of the QTL does not seem to correlate with the mean levels of the elements in the grains sampled from that location, suggesting that QEI are not being driven solely by altered availability of the elements in the soil. Local soil differences are less likely to be driving the QTL found for pairwise differences between two years at one location. Soil content should remain relatively similar from year to year at the same farm, suggesting that the loci identified by comparison between years and within location will encode components of elemental regulatory processes responsive to precipitation, temperature, or other weather changes. Experiments with more extensive weather and soil data, or carefully manipulated environmental contrasts, are needed to create models with additional covariates and precisely model environmental impacts.
Although the mapping intervals do not provide gene-level resolution, several QTL overlap with known elemental regulation genes, such as the QTL on chromosome 1 at 378 cM which coincides with ZEAMMB73_045160, an ortholog of the Arabidopsis molybdenum transporter, MOT1. We observe strong effects and replication of this QTL across nearly all environments, suggesting that this MOT1 ortholog plays a role in a variety of environments. Other large effect QTL found in several environments merit further investigation, as they may recapitulate important element-associated genes that have yet to be identified. Identification of the genes underlying these QTL and the gene/environmental variable pairs underlying the QEIs, will improve our understanding of the factors controlling plant elemental uptake and productivity. Given the high levels of variability that the interaction between genotype and environmental factors can induce in these traits, conventional breeding approaches that look for common responses across many different environments for a single trait may fail to improve the overall elemental content, necessitating rational approaches that include both genetic and environmental factors.

Conclusions
Here we have shown that the maize kernel ionome is determined by genetic and environmental factors, with a large number of genetic by environment interactions. Elemental profiling of the IBM population across 10 environments allowed us to capture environmentally driven variation in the ionome. Our QTL analysis on elements found mainly single-environment QTL, indicative of substantial genetic by environment interaction in establishment of the elemental composition of the maize grain. This approach, along with identification of QEI occurring both within a single location over different years and QEI between different locations, demonstrated that gene by environment interactions underlie elemental accumulation in maize kernels.