1,135 ionomes reveal the global pattern of leaf and seed mineral nutrient and trace element diversity in Arabidopsis thaliana

Soil is a heterogeneous reservoir of essential elements needed for plant growth and development. Plants have evolved mechanisms to balance their nutritional needs based on availability of nutrients. This has led to genetically based variation in the elemental composition, the ‘ionome’, of plants, both within and between species. We explore this natural variation using a panel of wild-collected, geographically wide-spread Arabidopsis thaliana accessions from the 1001 Genomes Project including over 1,135 accessions, and the 19 parental accessions of the Multi-parent Advanced Generation Inter-Cross (MAGIC) panel, all with full-genome sequences available. We present an experimental design pipeline for high-throughput ionomic screenings and analyses with improved normalisation procedures to account for errors and variability in conditions often encountered in large-scale, high-throughput data collection. We report quantiﬁcation of the complete leaf and seed ionome of the entire collection using this pipeline and a digital tool, Ion Explorer, to interact with the dataset. We describe the pattern of natural ionomic variation across the A. thaliana species and identify several accessions with extreme ionomic proﬁles. It forms a valuable resource for exploratory genetic mapping studies to identify genes underlying natural variation in leaf and seed ionome and genetic adaptation of plants to soil conditions.


INTRODUCTION
The ionome represents the sum of all mineral nutrient and trace elements of biological or environmental importance in an organism (Lahner et al., 2003). The coupling of the analyses of the ionome with genome-enabled genetics is termed ionomics (Salt, Baxter and Lahner, 2008). Ionomics has been successfully applied in plants and other organisms such as yeast (Eide et al., 2005;Yu et al., 2011Yu et al., , 2012, milk (Hadsell et al., 2018), human serum (Konz et al., 2017), mammalian organs (Ma et al., 2015) and human cell lines (Malinouski et al., 2014). Plants, being sessile organisms, have evolved intricate regulatory mechanisms to balance uptake and distribution of mineral nutrients and trace elements in response to physical and chemical variation of the soil. Arabidopsis thaliana is a non-crop, genetic model plant species extensively used in plant biology due to its small genome size and extensive genetic resources. Ionomic studies in this species have led to the identification of genes involved in numerous processes, such as Casparian strip formation (Hosmani et al., 2013;Kamiya et al., 2015), sphingolipid biosynthesis (Chao et al., 2011), phloem transport (Tian et al., 2010), vesicle trafficking (Gao et al., 2017), one-carbon metabolism  and iron homeostasis (Hindt et al., 2017). Further, there is an increasing body of organ-specific ionomic data that provide insights into the underlying functions (Baxter et al., 2014;McDowell et al., 2013;Pauli et al., 2018). These make ionomics on plants particularly interesting as it can reflect variation in soil mineral composition, altered function of ion transporters, the genes that control the concentration of essential or toxic elements and processes that indirectly impact their transport (Huang and Salt, 2016).
One of the first plant genome-wide association studies (GWAS) in A. thaliana explored the allelic variation in a selection of 96 wild-collected accessions and investigated 107 different traits, including the leaf ionome (Atwell et al., 2010). The study identified the gene HKT1;1 as being associated with phenotypic variation in leaf Na concentration, validating what had previously been observed in a bi-parental mapping population (Rus et al., 2006). An extension of this approach using the A. thaliana HapMap population of 349 accessions showed improved detection, identifying three causal genes underlying variation in the concentration of elements in leaves-MOT1 for Mo, HMA3 for Cd and HAC1 for As Chao et al., 2014b;Chao et al., 2012;Forsberg et al., 2015)-and reidentifying HKT1;1. In addition, the intraspecific ionomic variation of accessions has been used with success in crops like rice (Oryza spp.) (Pinson et al., 2015;Ricachenevsky et al., 2018;Yang et al., 2018), soybean (Glycine max)  and Brassica rapa (Bus et al., 2014;Thomas et al., 2016). It provides a promising basis to uncover the genetic and physiological mechanisms that regulate the elemental concentrations in plants, both between (Watanabe et al., 2007) and within species (Huang and Salt, 2016;Neugebauer et al., 2018), and the potentially adaptive significance of this variation (Busoms et al., 2018;Poormohammad Kiani et al., 2012). Further, it highlights the power of screening large collections of accessions for the identification of rare alleles associated with differences in the ionome. The 1001 Genomes Project made available detailed genome sequence information for a large diverse collection of around 1,135 wild-collected A. thaliana accessions from Europe, Central Asia, North Africa and North America (Alonso-Blanco et al., 2016;Cao et al., 2011;Weigel and Mott, 2009). With a threefold higher number of accessions and a genetic diversity of more than 10 million bi-allelic SNPs, we decided to experimentally characterise the full natural ionomic variation of this population to help uncover the genetic architecture underlying the natural ionomic variation within A. thaliana in both leaves and seeds. The A. thaliana seed ionome would be a unique dataset to verify if the genetic basis of the seed ionome in A. thaliana and important crops such as rice, maize (Zea Mays) and soybean are the same. Further, it forms a unique dataset to explore and better understand the relatedness between the leaf and seed ionomes.
Due to the size of the 1001 collection, large-scale ionomic experimental workflows were developed where plants were grown and analysed on a continuous basis over the course of several months. Previous ionomic studies on collections of natural accessions faced the challenge of suboptimal normalisation and experimental design that may have contributed to underpowered results (Atwell et al., 2010;Baxter, 2010;Chao et al., 2014a;Chao et al., 2012). This highlights the importance of a robust experimental design and sound normalisation of stochastic variation in experimental conditions. Normalisation methods can be divided into control-and sample-based methods (Birmingham et al., 2009). In most high-throughput ionomic experiments with A. thaliana, control-based normalisation techniques have been used to filter out differences between experimental blocks (plant growth trays) (Baxter, 2010;Chao et al., 2014b;Chao et al., 2012). Methodologically, high-throughput ionomic experiments on plants resemble high-throughput microplate assessments. For these kinds of experiments, a plethora of normalisation methods have been developed and reviewed to compensate for introduced noise into these kinds of experiments (Baryshnikova et al., 2010;Caraus et al., 2015;Dragiev et al., 2011;Malo et al., 2006;Murie et al., 2015;Wiles et al., 2008;Yu et al., 2011). Spatial bias can be compensated for by statistical techniques such as spatial smoothing, in which spatial bias is modelled based on samples and compensated for with generalised additive models (GAMs) (Hastie and Tibshirani, 1986) or LOESS curves (Baryshnikova et al., 2010;Murie et al., 2015). In this study, we propose an experimental design and a two-step normalisation procedure for high-throughput ionomic studies, which extend previously described normalisation procedures such as restricted maximum likelihood (REML) normalisation (Broadley et al., 2010), by using a samplebased approach to normalise for spatial variation between and within experimental blocks and a control-based approach to evaluate the normalisation.
With this work, we aimed to provide a unique referential dataset of this large natural population that forms a valuable phenotypic resource to enable large-scale GWAS to identify the genetic basis of ionomic traits. The study also explores the relationship between the leaf and seed ionome and identifies new extreme accessions that can be used to generate bi-and multi-parental populations for quantitative trait locus (QTL) analyses. Finally, the wide geographic range of this population, along with dense sampling at specific locations, can be a starting point for new local adaptation studies.

Validation of experimental design
We developed and optimised a high-throughput analysis pipeline for elemental profiling of 22 elements-Li, B, Na, Mg, P, S, K, Ca, Cr, Mn, Fe, Co, Ni, Cu, Zn, As, Se, Rb, Sr, Mo, Cd and Pb-in leaves and seeds and present data for 1,135 accessions of A. thaliana from the 1001 Genomes Project, which includes 180 Swedish accessions and the 19 parental lines of the MAGIC population (Alonso-Blanco et al., 2016;Kover et al., 2009;Long et al., 2013;Weigel and Mott, 2009).
We included four normalisation lines-Col-0, Cvi-0, Fab-2 and Ts-1-in our experimental design to correct for experimental variation in the ionomic profile between plant growth trays. We confirmed that the ionomic profiles of these selected accessions are distinct ( Figure S1(a)). We observed that Cvi-0 showed an extreme phenotype for several elements-high Mo, As and Cd, low Sr, Ca and Mgand was the most distinct accession among the four selected. To account for variation between and within trays we included four replicates per normalisation line, which were distributed evenly in growth trays on a row and column basis ( Figure S1(b,c)). To account for the impact of instrumental noise caused by drift in the sensitivity of the Inductively Coupled Plasma-Mass Spectrometer (ICP-MS) on the analysis, we applied a linear model that incorporates the effect of signal drift within a given run as a covariate to compensate for variation between ICP-MS runs. Our results indicate significant drift of the analytical signal within the ICP-MS runs for all elements except Na and Zn, and the drift patterns varied between ICP-MS runs. There were also significant differences in the sensitivity of the instrument across all elements. Once the analytical signal for each element was corrected for drift, we accounted for the limit of detection (LOD) for each element on the ICP-MS by excluding quantifications of a given element in a sample that had values below 3 standard deviations of the mean value of the blanks (Table 1). This resulted in the exclusion of a large number of values for elements, such as B (81.6%), Cr (70.9%), Ni (66.5%) and Pb (95.6%), which were therefore considered unreliable and excluded from further analyses.

Normalisation strategies employed and validated across the pipeline
To increase the robustness of the ionomic dataset, we utilised four common check-up lines-Col-0, Cvi-0, Fab-2 and Ts-1, rather than just Col-0. This provided a better fit for the data with a goodness of fit for the calculated versus measured weights of 92.3% for leaf samples and 83.6% for the seed samples. We found significant differences in elemental concentrations for all elements between the different trays used to grow the plants, as was reflected in similar patterns of variation between trays of the check-up lines ( Figures S2 and S3, Tables S1 and S2). The application of a tray-specific smoother using a GAM to account for the systematic noise successfully compensated for the differences within trays for all elements with no further within tray differences for either leaf or seed ionomes, except for Zn for the seed ionome (Tables S1 and S2).

Distribution of elemental data and identification of new accessions with extreme ionome profiles
We examined the leaf and seed elemental distribution across the accessions in the study. We found high levels of variation in the leaf elemental concentrations across the accessions, which included some new and already described extreme accessions for most elements appearing at the tail of the distribution Chao et al., 2013;Huang and Salt, 2016;Morrissey et al., 2009b;Rus et al., 2006) (Figure 1a). We were also able to detect several accessions with extreme leaf ionome phenotypes never described before by selecting the 10 top and bottom accessions for each leaf element concentration (Table S3). Among the 10 accessions with highest and lowest concentrations of each element some occurred multiple times, indicating a multi-element extreme ionomic profile. For some elements, these were chemical analogues, such as S and Se, K and Rb or Ca and Sr. However, even when not considering the chemical analogues, there were 10 accessions that stood out with more than three elements among these extremes, namely Stilo-1, Cvi-0, Pva-1, CS75436, Ham 1, Ala-0, Castelfed-2-201, UKID13, Ha-SB and Reg-0 ( Figure 2a). The check-up lines Cvi-0 and Fab-2 and the accession Ler-0 from the MAGIC parents appear in this list. High variation across accessions was observed for most elements quantified in the seeds, with new and already described extreme accessions appearing in the tails of the distribution (Figure 1b). Several accessions with extreme seed ionome phenotypes were identified by selecting the 10 top and bottom accessions for seed elemental concentration, with several being described for the first time (Table S4). Not considering chemical analogues, nine accessions stood out for having more than three elements among these extremes-Lip-0, Pa-1, Bar-1, Aln-30, Lam-0, Urd-1, Iso-4, Ven-0 and Bela-3 ( Figure 2b). There were 41 accessions that had two elements among these extremes for the leaf and seed data (Tables S3 and S4). The MAGIC parents Hi-0, Tsu-0 and Zu-0 appear in this list. We then calculated the broadsense heritability to estimate to what extent the variation in the leaf and seed ionomes across the accessions were due to genetic factors (Table 2). We see the heritability values being lower than reported in previous studies for the leaf ionome (Atwell et al., 2010;. The broad-sense heritability of elemental concentration in leaves across accessions ranged from 4.32% for Zn to 77.14% for Mo. For the seed ionome, heritability values ranged from 2.07% for Ni to 36.31% for Se (Table 2). These findings provide evidence of the power of using larger populations, such as the 1001 Genomes Project collection, to capture a more complete picture of the existing natural variation for traits such as the leaf and seed ionomes in A. thaliana and its potential to be replicated in other plant species.

Analysis of the leaf and seed ionomes
Elemental correlation analysis. To understand the relationship between the elemental concentrations, we performed a correlation analysis to cluster the elements based on unassisted hierarchical clustering and defined the groups based on K-means clustering ( Figure 3a). For the leaf ionome, the elements grouped into 10 clusters with a stronger association among elements within a cluster than between clusters. As expected, the chemical analogues S and Se, K and Rb or Ca and Sr had the highest correlation and grouped in the same clusters. The elements Li, Fe, Ca and Sr clustered together and were positively correlated between themselves, but negatively correlated with K, Rb, S and Se. Although not in the same cluster, some of these elements also showed a positive correlation with Mg. Other clusters of elements positively correlated were Cu, Zn, Mn and Cd. The elements Mg, Co, Mo and Na did not   The parameters shown are: mean concentration (ppm), standard deviation (SD, ppm), relative standard deviation (RSD, %) and broad-sense heritability (%). Principal component analysis. Given that different soils and environments may drive differences in the ionome, one might expect the ionomes of accessions to reflect their geographical origins. In order to explore the existence of these patterns, we performed a principal component analysis (PCA) on the leaf and seed ionomes as it highlights extreme phenotypes taking into account the ionome as a whole. First, the accessions were grouped based on similarities and differences in the ionome. Then, for signs of regional adaptation, they were grouped based on the geographical groupings described in a previous study on these accessions (Alonso-Blanco et al., 2016). For the leaf ionome, the first two principal components (PCs) together explained only 34.5% of the total variation in the leaf elements across the studied accessions ( Figure S4a). The PCA largely supports the clustering of the different elements found in the elemental correlation analysis. The chemical analogues S and Se, K and Rb or Ca and Sr and elements Cd, Zn and Cu are grouped together, likely due to low specificity of some uptake and transport proteins for these elements. The first PC explains 21% of the variation, mainly driven by variation in the concentrations of Cd, Mn, Zn, Cu and P. The second PC explains 13.5% of the variation in the leaf ionome, largely determined by the chemical analogues Ca and Sr, Se and S or K and Rb and Fe, Mg and Li. Some elements had opposing loadings in the PCA plot, which could be indicative of a crosstalk between their homeostases. The notable ones were P versus Mn and Mg versus Co. However, no strong correlation was found between these pairs of elements ( Figure 3a). Five of the 10 multi-element extreme accessions identified stand out as extremes in the PCA plot ( Figure S4a). However, the PCA was unable to identify strong signs of regional identity based on the geographical groups. Some extreme accessions for the same element that stand out in the PCA were from the same geographical group: Men-2 and Vas-0 from Spain with high Se and Panik-1 and Borsk-2 from Asia with high Sr and Ca. Interestingly, although the accessions Castelfed come from the same region and sampled very close to each other, they do not have a similar leaf ionomic profile, and do not cluster together in the PCA. For the seed ionome, we found that the first two PCs together explained 31.3% of the total variation in the seed ionome across the accessions ( Figure S4b). The first PC explained 18.4% of the variation, mainly driven by variation in the concentrations of Ca, Mg, Cu and Cd. The second PC explained 12.9% of the variation, largely determined by the elements Sr and Mn and the chemical analogues K and Rb. Indications of crosstalk were observed for Fe and Mo versus Cu and Mg or Zn and Co versus Li and As versus Cd. However, only As had a clear negative correlation with Cd. The accessions Lam-0 and Bar-1, previously identified as extremes for multiple elements (Huang and Salt, 2016), likewise stand out as extremes in the seed ionome PCA.

Clustering of accessions based on the ionome and evidence of local adaptation
To examine possible signs of regional adaptation, first we grouped the accessions based on similarities and differences in the ionome against the above-mentioned geographical groupings (Alonso-Blanco et al., 2016). In general, neither the leaf nor the seed ionome showed distinct clustering based on the geographical groups either looking at the total set of accessions or the 10 top and bottom extreme accessions (Figures S5a,b and S6a,b).
Geographic distribution of extreme accessions per element. To further explore the possible relationship between geographic regions and the leaf and seed ionomes, we looked at the geographic distribution of the 10 highest and 10 lowest extreme accessions for each leaf and seed element concentration. These were plotted onto the geographic location of the original collection sites (Figure 4a,b). We could confirm two clusters of accessions from North and South Sweden with a common extreme leaf ionome phenotype that also clustered in the heatmap ( Figure S6a). Accessions with extremely low K concentrations occurred together in the northern region of Sweden, and accessions with extremely low S concentrations occurred close to each other in the southern region of Sweden. The clusters identified with this approach consisted of accessions that occurred together in Azerbaijan at the Caspian Sea with extremely high S concentrations in leaves, at Castelfed with high Cd and Sr and in Spain with high K and Rb (Figure 4a,b). While using the same approach for the seed ionome we also identified new clusters of accessions not found with the heatmaps. A cluster of accessions from Spain, with extremely high K and low Sr concentrations in seeds, and several accessions with extremely high Se and low Zn concentrations in seeds distributed across Spain (Figure 4c,d). Accessions with extremely low Ni concentrations in seeds occurred together in the United Kingdom and accessions with extremely low Cd concentrations occurred in the Southeast of Sweden and in Azerbaijan at the Caspian Sea (Figure 4a,b). Finally, we could also confirm the cluster of accessions from the South of Sweden with extremely low Se concentrations in seeds found in the heatmaps (Figures S5b and S6b). We conclude, based on these observations, that the highthroughput ionomics approach using a large-scale worldwide collection was able to identify novel geographical patterns of variation based on analysis of the leaf and seed elemental concentrations.

Comparative analysis between different tissues-the leaf and seed ionomes
Correlation analysis of the leaf and seed elemental profiles. We explored the relationship between the concentrations of different elements in different tissues of the plant across all accessions and found the elements to group into 13 clusters ( Figure 5). Clusters identified when considering the seed and leaf ionomic data together largely support the clusters also identified when analysing the leaf and seed tissues independently. The exceptions, however, were S, Mo and Na, whose seed and leaf concentrations are in the same cluster and are positively correlated.
PCA analysis of the leaf and seed ionomes combined. To examine the relationship between the leaf and seed ionomes, we performed a PCA on the element concentrations of the two tissues. The first two PCs together explained 20.6% of the total variation observed and gave an indication of accessions with extreme ionomic phenotypes for both tissues besides supporting some of the clusters of the different elements previously described in this paper ( Figure 6). The first PC explains 12% of the variation and is mainly driven by Mn, Cd, Sr and K concentrations in seeds, while the second PC explains 8.6% of the variation, largely driven by S, Se, Ca and Sr concentrations in leaves. Even though the PCA analyses take the ionomes for both tissues as a whole, elements with opposing loadings, which indicate possible crosstalk for homeostasis, were only found for the same tissue. Several elemental concentrations in both seeds and leaves had extremely low PCA loadings. Most of these were the ones that correlated weakly with other elements, when considering the seed and leaf ionomes either together or independently.
Identification of accessions with extreme profiles in both leaf and seed ionomes. Although not all elements in the ionome showed a strong relationship between the leaf and seed ionomes, we were able to identify several accessions with extreme leaf and seed ionome profiles. The extreme accessions identified were separated in two groups: (i) accessions with an extreme phenotype for the same element or chemical analogues in both seed and leaf (Figure 7a) and (ii) accessions which are among the 10 top and bottom accessions for each element concentration in both leaves and seeds across the accessions in the study (Figure 7b). Four out of these accessions also stood out as seed and leaf extremes in the PCA plot: T980, Sim1, Lam-0 and Castelfed-2-201.  to compare the Z-scores of all elements of select accessions, correlation and PCA that considers the filtered data as a whole. One can interact with the leaf and seed data independently and in parallel. The complete and filtered dataset can be downloaded as a .csv file, while the maps, correlograms and radar and PCA plots generated can be downloaded as image files. Further, the user can also upload data to the tool to perform comparative analyses against the dataset. The user functionality of the tool is constantly being updated.

DISCUSSION
Mineral nutrients are critical for a plant's normal growth and development. These nutrients are primarily acquired from the soil, which is a highly heterogeneous mixture that varies in chemical and physical properties both spatially and temporally. Much attention has been given to understanding the mechanisms that underlie the plastic responses that allow plants to respond to fluctuations in mineral nutrients over time, during their lifetimes (e.g. de Bang et al., 2020). Much less attention has been given to investigating the molecular basis of evolutionary adaptive genetic change driven by selection to a particular edaphic property (physical, chemical or biological property of soil) over many generations. This question goes back to the beginnings of evolutionary biology, with both Wallace and Darwin recognising that edaphic discontinuities such as boundaries between soil types produce strong natural selection (Darwin and Bynum, 1859;Wallace and Darwin, 1858). Plant adaptation to unique soils such as serpentine  and mine sites is well described (Brady, Kruckeberg and Bradshaw, 2005;Macnair, 1993), but a molecular genetic understanding of edaphic adaptation lags behind, though progress is being made (e.g. Arnold et al., 2016;Busoms et al., 2018;Selby and Willis, 2018). Genetically determined natural variation in the A. thaliana ionome is known to exist, and at least in certain cases has been suggested to be adaptive (Baxter and Dilkes, 2012;Poormohammad Kiani et al., 2012). The possibility therefore stands that natural variation in the ionome reflects the evolutionary adaptation of A. thaliana to the varied edaphic conditions of its dispersed habitats, across much of the northern hemisphere. Mapping this natural ionomic variation to the landscape provides a pathway for the identification of both the adaptive loci upon which selection is acting and the specific nature of the selective agents.
To help achieve this we present the largest available dataset to date of the leaf and seed ionomes of over 1,135 accessions of A. thaliana from a diverse geographical range, and that were a part of the 1001 Genomes Project

Experimental design and data analyses of high-throughput ionomic profiling experiments
There have been ongoing efforts to optimise the efficiency of ionomics workflows and the quality of the phenotypic data produced. In this study, we present an improved method for large-scale ionomic studies, including experimental design and data analysis. Large-scale ionomic experimental workflows where plants are grown and analysed on a continuous basis often span several months. Given the length of time required to grow and analyse thousands of plants, even with best efforts to maintain standard conditions, variation in environmental conditions (temperature, watering, etc.) between experiments is unavoidable. Therefore, it is essential to design experiments to include multiple replicates of each accession spread across multiple plant growth trays to take into account the interaction between tray and accessions, which is what we have done in this study. Our observations emphasise the need for randomisation and robust normalisation strategies to reduce type I (false positive) and II (false negative) errors. Spatial bias introduced during plant growth can be compensated by statistical techniques such as spatial smoothing in which spatial bias is modelled and compensated with GAMs (Hastie and Tibshirani, 1986) or LOESS curves (Baryshnikova et al., 2010;Murie et al., 2015). Previous studies on data normalisation have found that spatial smoothing techniques were able to properly model this type of noise (Baryshnikova et al., 2010;Murie et al., 2015). In this study, we successfully expanded the approach to plant ionomic research and were able to better compensate for experimental noise than previous high-throughput ionomic screenings applying a control-based normalisation approach without modelling within tray differences both for A. thaliana (Atwell et al., 2010;Baxter et al., 2008Baxter et al., , 20102014b;Chao et al., 2012;Forsberg et al., 2015) and for other crops (Pinson et al., 2015;Ricachenevsky et al., 2018;Yang et al., 2018).

Heritable ionomic variation
In genetic studies, the broad-sense heritability values indicate the relative importance of genes and environment to the phenotypic variation of traits within and across populations (Visscher, Hill and Wray, 2008). Here, the heritability of variation in the leaf ionome in general is lower than previously reported for most elements except Mo. Although lower, we found correlations (0.62-0.86) between our observed heritability values and previous studies with 96 accessions, or the 360 HapMap population (Atwell et al., 2010;Baxter, 2010;Chao et al., 2014b) and the study with two recombinant inbred line populations (Buescher et al., 2010). For the seed ionome, even though we observed high variability, our heritability values were in general lower than those reported in previous studies with A. thaliana  and other crops (Baxter et al., 2014;Pinson et al., 2015). This highlights the consistency in the general trend of elements with high and low heritability values. The possible cause for the lower heritability values in our studies could be the randomisation of replicate accessions across the growth trays implemented in our experimental design but not in previous studies (Atwell et al., 2010;Baxter et al., 2010;Buescher et al., 2010;Chao et al., 2014b).

Extreme ionomic accessions-known and new
Using a larger population from geographically dispersed locations we were able to identify several new singleand multi-element extreme accessions not reported previously (Tables S3 and S4, Figures 2 and 7). This demonstrates the ability of the dataset to capture a more complete picture of the existing natural variation of the ionome. Previously identified A. thaliana accessions with an extreme ionomic phenotype have been used in studies to generate bi-parental mapping populations and identify key genes involved in the homeostasis of Cd (Chao et al., 2012), Co (Morrissey et al., 2009a), Na Rus et al., 2006), Mo , As (Chao et al., 2014b), S and Se (Chao et al., 2014a) and Zn (Chen et al., 2018) in plants. New accessions with extreme ionomic phenotypes identified in this study form an important starting point for the development of new bi-and multiparent experimental populations for further genetic studies. Their unique ionomic phenotypes can aid in the study of multi-element nutrient transport pathways, their genetic bases and how plants respond to different growth conditions and nutritional stresses. Biological insights from multi-element species-wide ionomic data The relationship and behaviour of multiple elements in the plant can be determined by three potential processeschemical analogues, elemental transporters or common biological mechanisms-that regulate the concentration of several elements at the same time (Briat et al., 2015;Salt, Baxter and Lahner, 2008;Vreugdenhil et al., 2004). Strong correlations and relationships were observed in the leaf ionome of K and Rb, Ca and Sr or S and Se and Li, Fe and Mg among the accessions. K and Rb, Ca and Sr or S and Se are chemical analogues, and the strong correlations can be attributed to a robust normalisation strategy applied and in the biological context a lack in selectivity to take up and transport these elements in the plant (Broadley et al., 2010;White, 2001). On the other hand, some elements, such as Na, Mo, Co and As, were weakly or not correlated with other elements, and had extremely low PCA loadings, indicative that their homeostasis is highly element-specific. Leaf concentrations of Mn, Cd, Cu and Zn across the studied accessions showed strong associations between these elements which have been previously described and considered to be likely a result of shared transporter proteins. For example, transporters of the NRAMP family show affinity to the elements Mn, Zn, Cd and Cu (Cailliatte et al., 2010;Oomen et al., 2009). IRT1, primarily identified as an Fe transporter, can also take up Mn, Zn, Co and Cd (Connolly, Fett and Guerinot, 2002;Korshunova et al., 1999); ZIP4 can take up Zn and Cu (Wintz et al., 2003); HMA3 can take up Cd, Zn and Co (Hussain et al., no date;Morel et al., 2008;Chao et al., 2012) and HMA4 can take up Cd and Zn (Hussain et al., no date;Mills et al., 2003Mills et al., , 2005Verret et al., 2004).
Other ionomic patterns have also previously been identified. A genome-wide scan of the effects of all genes on the ionome of the budding yeast Saccharomyces cerevisiae (yeast) revealed 26 possible ionomic clusters when all genes were evaluated by knockout (Yu et al., 2012). This number of clusters is small compared to the total number of possibilities, and the authors conclude that 'the ionome only has a limited number of sets points that are consistent with cell viability'. These fundamental biological trade-offs in the ionome limit the potential natural ionomic variation that can occur. Major ionomic gene clusters from this study in yeast reveal the importance of vesicular sorting (Cd, Mn, Mo and Na), vascular acidification (P, Mg and Co) and mitochondrial function (Mn and K) (Yu et al., 2012). In plants, various studies have started to reveal the molecular basis of observed ionomic patterns. Perturbation in vesicular trafficking causes a distinctive ionomic phenotype involving Na, Mn, Fe, Zn and Mo and is due to protein sorting defects disrupting plasmodesmata and ion-transporter localisation (Gao et al., 2017). The regulation of the key nutrients P, S, Zn and Fe has also been suggested to be coordinated by the MYB-like transcription factor Phosphate Starvation Response 1 (PHR1) (Briat et al., 2015). The Ca antiporters CAX1 and CAX3 appear to play a redundant function regulating a distinct ionomic pattern of P, Ca, Mn, Zn and Mg via a PHO2-independent impact also on phosphate homeostasis (Cheng et al., 2005;Liu et al., 2011).
The seed ionome dataset could form the basis to uncover mechanisms regulating mineral nutrient and trace element concentrations in the seed. The lack of knowledge in this field is one of the major constraints to the development of biofortified grains with high concentrations of essential mineral nutrients such as Fe, Zn, Ca and Se and avoid the sometimes consequent accumulation of trace elements such as As and Cd (Bouis and Saltzman, 2017;McDowell et al., 2013;Mendoza-C ozatl et al., 2011;Waters and Sankaran, 2011;Yoneyama, Ishikawa and Fujimaki, 2015). Although within the same tissue several elements were correlated, across the leaf and seed only S, Mo and Na were positively correlated, consistent with observations from a previous study using 96 accessions . This suggests that although the seeds depend on vegetative tissue as a source of nutrients, the regulatory control mechanisms over the seed and leaf ionomes are distinct.

Evidence of species-wide local adaptation
Although several genes have been identified to regulate the ionome in A. thaliana, far less is known about their role in local adaptation (Huang and Salt, 2016). Some of the detailed studies on the influence of soil or climate on the A. thaliana ionome have been related to saline soils, where the ability of these plants to accumulate Na in leaves has been shown based on their distance from the coast/salinity level Busoms et al., 2015;Rus et al., 2006). Another study demonstrates the influence of Mo levels in native soils of West Asian accessions on Mo homeostasis (Poormohammad Kiani et al., 2012). Our attempts to identify associations between the leaf and seed ionomes and geographic regions did not lead to any clear trends. Previous local adaptation studies have highlighted the heterogeneous nature of soils and have suggested local adaptation is likely to be found at a small to very small geographic scale (Busoms et al., 2015;Macel et al., 2007). Also, the ecology of A. thaliana as a ruderal species often occurring in urban environments with a large heterogeneity in local soil composition may contribute to that (Guilland et al., 2018;Hovick and Whitney, 2019;Takou et al., 2019). The most prominent signals were found in the south of Sweden and at the Caspian Sea shore, with accessions from these areas having extreme low leaf S and Se levels. The presence of a large regional collection of 180 accessions from Sweden classified into regional groups as North and South Sweden may have contributed to our observation of an apparent association between location of origin and the leaf ionome of some Swedish accessions. New accessions identified in this study with not only high, but also low elemental concentrations in leaves and seeds hold potential for further study to identify novel alleles or genes involved in element homeostasis and local adaptation.
Taken together, this study presents a robust experimental design pipeline that allows the study of thousands of accessions with a normalisation strategy that deals with systematic noise within and between plant growth trays and accessions. The Ion Explorer tool holds potential to become a comprehensive resource for visualisation and interpretation of ionomic data. We envisage addition of further datasets and views to perform comparative analysis with other experimental conditions and links to other databases to increase functionality and output. The power of the presented ionomic dataset lies in the ability to undertake large-scale genetic studies, including GWAS and linkage mapping using recombinant populations, to identify new genes and alleles that regulate the plant ionome. The use of the newly identified extreme elemental accessions to generate bi-and multi-parent populations will help pave the way to uncover nutrient transport pathways in the plant and its response to different growth conditions and nutritional stresses. Finally, we present indications that local adaptation of plants to soil conditions is reflected to a limited degree in their ionome, mostly at small geographical scales. The identification of promising regions would be an interesting focus for future local adaptation studies. Such information can help improve our understanding of how plants efficiently acquire mineral nutrients from the soil and adapt to their local soil habitat.

Plant materials
The plants used for leaf and seed ionomic analysis included 1,135 accessions of A. thaliana from the 1001 Genomes Project, which also includes 180 Swedish accessions and the 19 parental lines of the MAGIC population (Alonso-Blanco et al., 2016;Kover et al., 2009;Long et al., 2013;Weigel and Mott, 2009

Experimental design and growth conditions
Leaf ionome. Plants were grown in 32-mm 7C Jiffy-7 â peat pellets (http://www.jiffypot.com/) pre-soaked in a solution (104 pellets in 2.6 L) containing 0.142 µM Na 2 HAsO 4 •7H 2 O, 5.935 µM Cd (NO 3 ) 2 •4H 2 O, 0.534 µM Co(NO 3 ) 2 •6H 2 O, 0.053 µM LiNO 3 , 0.534 µM Ni(NO 3 ) 2 •6H 2 O, 0.427 µM K 2 SeO 4 , 0.106 µM Sr(NO 3 ) 2 (98%), 0.0774 µM RbNO 3 to provide these elements in trace amounts to enable detection in the plants. The final concentrations of these elements in dry soil were 5.6 ppm As, 0.2 ppm Cd, 1.1 ppm Co, 1.3 ppm Li, 1.1 ppm Ni, 1.9 ppm Se, 8.7 ppm Sr and 11.7 ppm Rb. The Jiffy-7 â peat pellets were evenly distributed in trays with 104 places, and the trays lined with capillary matting soaked in 1.5 L of 18 MO deionised water. Three seeds of each A. thaliana accession were sown in the pre-soaked Jiffy-7 â pellets and stratified at 4°C in the dark for 1 week to help synchronise germination. Germinated seedlings were thinned to one after 1 week. A growth pipeline was established with one tray sown per day, and after 4 weeks, one tray harvested per day. Each tray had 104 plants distributed randomly of which 88 were a single replicate of different accessions and four normalisation lines replicated four times each ( Figure S1b). The accessions used for normalisation lines were Col-0, Cvi-0, Fab-2 and Ts-1. Six replicates of each accession were grown, each in a separate tray comprising 120 plant growth trays grown over a 37-week period. Plants were cultivated in a controlled growth room with 10-h days (approximately 100 µmol photons m À2 sec À1 )/14-h nights and temperatures of approximately 25°C (day) and 17°C ( Seed ionome. Plants were grown in 32-mm 7C Jiffy-7 â peat pellets prepared as above. The experiment was divided into four phases: Phase 1: Sowing-three seeds of each accession were sown in pre-soaked Jiffy-7 â pellets and stratified at 4°C in the dark for one week. Germinated seedlings were thinned to one after 1 week. Phase 2: Initial growth-Plants were grown for 2 weeks in a greenhouse with supplemental lighting set at 12 h (day)/12 h (night) at 26°C (day)/16°C (night). Plants were bottom-watered once a week with 0.259 Hoagland solution in which Fe was supplied as 10 µM Fe-HBED. Phase 3: Vernalisation-Plants were placed in a cold room for 8 weeks with 10-h days (50 µmol photons m À2 sec À1 )/14-h nights and temperatures of 4°C (day)/2°C (night). They were bottom-watered twice with tap water. Phase 4: Final growth-Vernalised plants were returned to the greenhouse and placed randomly in 61 new trays distributed across four benches ( Figure S1c). Plants were bottom-watered with 0.259 Hoagland solution four times in the first 2 weeks and kept in the greenhouse for approximately 2 months until all seed pods were dry. Dry seeds were harvested and stored in paper envelopes. In total six replicates of each accession were grown and seeds were harvested.

Ionomic profiling by inductively coupled plasma mass spectroscopy analysis
Leaf and seed ionomes. The samples were analysed using the method described by Danku et al., (2013) with some modifications. For the leaf ionome, two or more leaves (2-5 mg dry weight) were harvested per plant, washed with 18 MO deionised water and placed in Pyrex digestion tubes. The leaf samples were dried in the Pyrex tubes for 20 h at 88°C and digested with 1 ml of concentrated nitric acid (HNO 3 ; Baker Instra-Analyzed) in dry block heaters at 115°C for 4 h. For the seed ionome 2-5 mg of seeds from each plant was placed in Pyrex tubes and digested with 1.5 ml of concentrated HNO 3 in dry block heaters at 115°C for 5 h. After digestion, samples were diluted to 10 ml with 18 MO Milli-Q filtered water and elemental analysis was performed using an Inductively Coupled Plasma -Mass Spectrometer (ICP-MS) (Perkin Elmer, NexION 300D, Shropshire, UK) in the standard mode and monitored for 22 elements (Li, B, Na, Mg, P, S, K, Ca, Cr, Mn, Fe, Co, Ni, Cu, Zn, As, Se, Rb, Sr, Mo, Cd and Pb). In order to minimise the effect of analytical drift within and between ICP-MS runs, two liquid reference solutions composed of pooled digested leaf or seed samples of the first 10 trays of the leaf and seed experiments, respectively, were prepared and analysed every 10 samples for all sample sets run on the ICP-MS. To compensate for the effect of variation within and between the ICP-MS runs, we used a linear model that incorporated the effect of within run drift as a covariate and between run differences and their interaction (Holmberg and Artursson, 2004). To calculate the LOD of the instrument per element, 40 process blanks were included in each ICP-MS sample set (eight blanks per plant growth tray). The blanks consisted of empty tubes that went through all the sample preparation processes alongside the samples. The detection limit was calculated based on the blanks as follows: LOD = (mean value of a given element in the blank) + 3 9 (standard deviation over all 40 blanks) (Lever et al., 2014) and samples below the detection limit were set as empty. Calibration standards (with Indium internal standard) were prepared from single-element standards solutions.

Data analysis
Weight calculation. Leaf ionome: For every plant growth tray, the dry weights of seven reference samples were measured and used to estimate the weight and final elemental concentration of the remaining samples based on a heuristic algorithm that uses the best quantified elements in a sample to calculate the weight and final elemental concentrations of the samples (Lahner et al., 2003). The approach was extended to include all four of the checkup lines Col-0, Cvi-0, Fab-2 and Ts-1, rather than just Col-0, for improved normalisation. Based on this, the dry weight of seven samples of each of the check-up lines was recorded per tray. For the calculations, only the elements with a relative standard deviation of 50 or less were included. The quality of the calculated weights was verified by performing a linear regression analysis between the calculated weights of the weighed samples and their measured weight.
Seed ionome: ICP-MS data processing for drift correction and calculation of the detection limits and sample weights were performed as described above for leaf ionome data. The weight of 14 samples in each tray was measured and used for the weight calculation heuristic algorithm method.
Normalisation. In order to reduce the effect of external factors and to account for variation in growth conditions, we used a normalisation approach that compensates for variation both between and within trays. First, per element, the outlier samples that were 4.5 times higher than the interquartile range across all trays were excluded from the analysis. This threshold was chosen over the more standard 1.5 times so as to not exclude accessions with an extreme phenotype from the analysis. In addition, there were idiosyncrasies observed in some trays for some elements due to, for instance, issues with the ICP-MS instrument, which caused a large number of samples in a tray to drop below the LOD. These trays were excluded based on elemental concentrations of the check-up lines, that is, if both the upper and lower quartile across all four check-up lines combined were zero. Next, to test for spatial patterns we applied a linear model with elemental concentrations as 'response variable' and row, column, tray and their interactions as 'explanatory variables' (Dragiev, Nadon and Makarenkov, 2011). Next, we modelled per element per tray with a spatial smoother using a GAM in R (Version 3.2.2; R Core Team, 2015) using the package mgcv (Version 3.2.2; Wood, 2011). We normalised for spatial patterns by subtracting the GAM predicted values for a tray from the un-normalised values of the tray and subsequently adding the overall mean predicted values. To evaluate the effectiveness of the normalisation we first tested if spatial patterns were still observable by repeating the linear model on the normalised values. Then, we analysed whether the normalisation reduced the variance among replicates of the accessions by constructing a linear model that compared per element the 'variance for all accessions of the unnormalised values', 'ionomic values normalised based on tray means' and 'spatial smoother normalised values'. After normalisation, the elemental values of the ionome were extracted per accession by calculating the Best Linear Unbiased Predictors (BLUPs) using an REML mixed model with only accessions as random effect and intercept as fixed effect (Table 1). The BLUPs allowed phenotypes to be predicted from genotypes, which is important for fitting functional structural plant models. For the seed ionome, to compensate for both between and within tray variation we used the same normalisation approach described above for leaf ionomic data.
Ion Explorer. The interactive tool was developed using R shiny, all plotting was performed using a combination of leaflet (to render the map), ggplot2 (to plot the heatmaps, histograms and PCA results) and plotly (to add interactive functionality and draw the radar plots). The user interface should work across most devices and browsers but has been optimised for devices with a minimum display width of 800 pixels. The code base is publicly available at https://bitbucket.org/ADAC_UoN/dr000081-web-service-ionomeseed-and-leaf-map/.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article. Figure S1. (a) Radar plot representing the elemental leaf ionome profile of the Arabidopsis thaliana accessions Col-0, Ts-1, Cvi-0 and Fab-2, which were used to normalise the experimental data. Axes display Z-scores calculated per element. (b) Distribution of normalisation lines over trays in the leaf experiments. (c) Randomisation map of normalisation lines used for the seed ionome experiment. Figure S2.    Figure S5. Heatmap that shows the hierarchical clustering of analysed accessions based on their leaf (a) and seed (b) ionomes. The top row represents the regional groups based on genetic similarities and geographic distance (Alonso-Blanco et al., 2016). Figure S6. Heatmap showing the hierarchical clustering of the accessions studied based on 10 top and bottom extreme accessions for the leaf (a) and seed (b) ionomes. The top row represents the regional groups based on genetic similarities and geographic distances (Alonso-Blanco et al., 2016). Table S1. Output of the linear model on row, column, tray and their interaction on the leaf ionome of the accessions studied. Table S2. Output of the linear model on row, column, tray and their interaction on the seed ionome of the accessions studied. Table S3. Overview of the 10 accessions with the highest and lowest concentrations of each element measured in leaves across the studied A. thaliana accessions. Table S4. Overview of the 10 accessions with the highest and lowest concentrations of each element measured in seeds across analysed A. thaliana accessions.

OPEN RESEARCH BADGES
This article has earned an Open Data Badge for making publicly available the digitally shareable data necessary to reproduce the reported results.