Genomic analysis of ionome-related QTLs in Arabidopsis thaliana

Ionome contributes to maintain cell integrity and acts as cofactors for catalyzing regulatory pathways. Identifying ionome contributing genomic regions provides a practical framework to dissect the genetic architecture of ionomic traits for use in biofortification. Meta-QTL (MQTL) analysis is a robust method to discover stable genomic regions for traits regardless of the genetic background. This study used information of 483 QTLs for ionomic traits identified from 12 populations for MQTL analysis in Arabidopsis thaliana. The selected QTLs were projected onto the newly constructed genetic consensus map and 33 MQTLs distributed on A. thaliana chromosomes were identified. The average confidence interval (CI) of the drafted MQTLs was 1.30 cM, reduced eight folds from a mean CI of 10.88 cM for the original QTLs. Four MQTLs were considered as stable MQTLs over different genetic backgrounds and environments. In parallel to the gene density over the A. thaliana genome, the genomic distribution of MQTLs over the genetic and physical maps indicated the highest density at non- and sub-telomeric chromosomal regions, respectively. Several candidate genes identified in the MQTLs intervals were associated with ion transportation, tolerance, and homeostasis. The genomic context of the identified MQTLs suggested nine chromosomal regions for Zn, Mn, and Fe control. The QTLs for potassium (K) and phosphorus (P) were the most frequently co-located with Zn (78.3%), Mn (76.2%), and Fe (88.2% and 70.6%) QTLs. The current MQTL analysis demonstrates that meta-QTL analysis is cheaper than, and as informative as genome-wide association study (GWAS) in refining the known QTLs.

www.nature.com/scientificreports/ proteomics, or metabolomics) levels and highlights critical concepts for assessment of the ionome 12,13 . Thus, ionomics should depict the functional status of a complex biological organism in both a quantitative and qualitative pattern of elements in various components of the organism 10 . The ionomics method refers to quantitative analysis of ionomic traits and their changes in living organisms in response to developmental state, physiological stimuli and genetic modifications. In addition to ionomics, the application of practical genetic analysis such as the genetic mapping approach leads to an understanding of the genetic architecture of ionomic traits. The use of molecularly well-characterized and diverse germplasms of model species like Arabidopsis thaliana expands the knowledge about ionomic variation and their associated genomic regions. Various QTLs have been identified using linkage mapping approaches for numerous important ionomic traits in plants [14][15][16][17][18][19] . Traditional QTL mapping has been successful in identifying genomic regions controlling the variation of ionomic traits. In traditional QTL mapping models, the genetic effects of QTL identified in an individual population may not be validated in other genetic backgrounds or environments and more often, only a limited number of traits can be tested in a single study 20 . However, meta-quantitative trait loci (Meta-QTL or MQTL) overcomes the limitation posed by traditional QTL analysis 21,22 . In the meta-QTL model, QTL data from independent studies and populations are used to detect consensus QTL regions known as Meta-QTL 23 . Further, the information on QTLs derived from various population structure, origin and sizes can be used to refine the known QTLs regardless of the genetic backgrounds, marker density and phenotypic variations 21,24,25 . In this way, a meta-QTL analysis helps to narrow down the confidence interval (CI) of the initial QTLs identified in independent populations and lays the foundation for a better understanding of traits underlying a QTL region than what is possible in independent QTL mapping studies 26 . Moreover, the MQTL analysis helps to validate the genetic association of loci identified by genome-wide association study (GWAS) approach [27][28][29] .
Despite the wide use of meta-QTL analysis in plant species [30][31][32][33][34][35][36][37][38][39] , little is known about ionomic meta-QTL analysis in A. thaliana. This study aims to perform a meta-QTL analysis to uncover genomic control of ionome related traits in A. thaliana. Due to the importance of essential elements in biofortification studies, the emphasis of the current study is placed on the assessment of genomic regions controlling trace elements, especially Zn, Fe, and Mn. In addition, the ability of MQTL analysis to validate the loci associated with micronutrients identified by GWAS 40 was assessed in our study.

Materials and methods
Development of QTL database. A database containing information on 483 QTLs for 38 different traits ( Fig. 1) derived from 12 segregating populations was developed. The number of QTL per study ranged from 5 to 83 QTLs that were used to find the most stable genomic regions controlling ionomic quantitative traits in A. thaliana (Table 1). Except for those lacking proper genetic map and QTL-related information, QTL mapping studies including important quantitative traits with different markers in A. thaliana were used for the identification of MQTLs. The information of populations used in the meta-QTL analysis is summarized in Table 1. The number of QTLs for each ionomic trait and their distribution on A. thaliana chromosomes are presented in Fig. 1 and Supplementary Figure S1. The assessed ionomic QTLs were distributed across different chromosomes (χ (4) 2 = 17.40, P = 0.0016) with chromosome 1 showing the highest number of QTLs (119 QTLs), followed by chromosome 5 (117 QTLs). Chromosomes 4 and 2 had the lowest number of QTLs (Fig. 1). The chromosomes 1 and 5 with 9 and 7 QTLs possessed the highest QTL numbers for Mn and Fe, respectively (Fig. 1). The highest number of QTLs for Zn were located on chromosome 5 with 14 QTLs. Details regarding the position, the proportion of phenotypic variance (% Expl.), and the log of odds ratios (LOD score) for QTLs were used for the analysis of meta-QTLs.
Constructing consensus genetic map and QTL projection. A consensus genetic map was constructed by integrating 12 map files. The constructed map file for each population included information about the type of cross, population size, map function, map units and the position of DNA markers used in different linkage groups. For map comparisons, the genome lengths from the individual and consensus genetic maps were calculated according to Hubert and Hedgecock 45 . First, the average spacing (χ) between markers was calculated by dividing the total length of all linkage groups by the number of intervals (number of markers minus the number of linkage groups). Then, the genome lengths for individual and consensus linkage maps were calculated by adding 2χ to the length of each linkage group to account for terminal chromosome regions 46 . After the development of the final consensus genetic map, the individual QTLs obtained from each population were projected onto the consensus genetic map in BioMercator v 4.2 24,25 . Meta-QTL analysis. The BioMercator v4.2 was 24,25 used to perform the meta-QTL analysis (available at https:// urgi. versa illes. inra. fr/ Tools/ BioMe rcator-V4). The Veyrieras et al. 22 approach was used for meta-QTL analysis where the number of QTLs in each chromosome was higher than 10 and the best QTL model was selected based on Akaike information criterion (AIC) 22 , corrected AIC criterion (AICc) 22 , AIC 3 candidate (AIC3) 22 , Bayesian information criterion (BIC) 22 , and Approximate Weight of Evidence (AWE) 22 . The best QTL model shows lower values in at least three of these five criterions. Consensus QTLs from the optimum model was regarded as meta-QTLs (MQTLs). The MQTL position and distribution in each linkage group were presented as a heatmap using the RIdeogram R package 47 (Fig. 2).
The variation of QTLs of traits over genetic consensus map and the MQTLs density in both genetic and physical maps regions were estimated following the study by Martinez et al. 32 . The QTL and MQTL densities on the genetic map were evaluated by counting the number of QTLs and MQTLs in every 50 cM interval across the A. thaliana genome, starting from the centromeric region towards the ends of the chromosome (Fig. 3). The same procedure was followed for estimating MQTLs density in 10 (Fig. 4). The genetic loci highly associated with grain Zn, Fe, and Mn (as the most targeted mineral elements in biofortification studies) in the A. thaliana GWAS study was retrieved and also compared with the MQTLs and the results of GWAS based on SNP markers 40 . The genomic positions of MQTLs and associated loci reported in the A. thaliana GWAS study 40 were compared based on the A. thaliana genome.

Functional candidate genes (CGs) in MQTLs interval.
To uncover functional CGs controlling the ionomic traits (especially for grain Zn and Fe content), the sequences of the flanking markers for each MQTL were used to identify the corresponding genomic interval on the A. thaliana genome. The sequences of flanking markers were retrieved from the TAIR10 genome and BLAST was used to identify the precise genomic location. For those flanking markers lacking a definite position on the A. thaliana genome, the closest marker on the consensus map was selected to identify the MQTL position on the genome. The corresponding gene annotation was retrieved from the gff file containing gene coordinates for TAIR10 in addition to the Uniprot ID corresponding to the A. thaliana reference genome (TAIR10).
Trait analysis within MQTL. To identify ionomic traits within MQTL regions, the MQTLs result was converted into binary scores (0 or 1) based on the absence/presence of an individual trait QTL within an MQTL region. Hence, the number of times a trait was present within an MQTL, the number of times more than one QTL for a trait was present within MQTL (implying confirmation of the QTL), and the number of times traits co-localized within an MQTL were calculated. A chi-squared test with a degree of freedom (DF) = 1 was per-

Results
Genetic consensus map and meta-QTL analysis. The original linkage map consisted of five linkage groups with an average of 118 DNA markers and covered 485.8 cM of the A. thaliana genome ( Table 2). The constructed consensus genetic map contained 725 markers that covered 514.01 cM of the A. thaliana genome ( Table 2). All the 483 initial QTLs were successfully projected onto the consensus genetic map and used for meta-QTL analysis. The meta-QTL analysis confined these QTLs into 33 MQTLs with QTLs originating from at least two populations for all the ionomic traits. The distribution of MQTLs for each ionomic trait on each chromosome is presented in Table 3  www.nature.com/scientificreports/  www.nature.com/scientificreports/ studies followed by MQTL-7/Chr1, MQTL-1/Chr3, and MQTL-3/Chr5 with 32, 28 and 28 initial QTLs derived from nine, seven, and nine independent population studies, respectively ( Table 3). The results showed that the "CH.192L-Col", "HH.357L", "AD.191L-Col", "SNP135", "SNP105", "DF.252L", "m291" and "BH.127L" markers were linked to these stable MQTLs. The current meta-QTL analysis reduced the average confidence interval (CI) up to eight folds with an average of 1.30 cM in MQTLs in comparison to the mean CI of the original QTLs (10.88 cM). Among the detected MQTLs, the CI of 12 MQTLs (36.4%) was reduced to < 1 cM (Table 3).

Gene density and distribution of MQTLs and QTLs on chromosomes. The distribution patterns of
QTLs on the consensus genetic map and MQTLs on the genetic and physical maps were investigated and compared with gene density (Figs. 3, 4). The overview of the distribution of the assessed QTLs for different ionomic traits exhibited that the non-telomeric regions possessed majority of QTLs and MQTLs in their intervals. Chromosomes 2, and 4 harbored relatively all of the QTLs on their non-telomeric regions (Figs. 2, 3b). A large number of QTLs were projected over the sub-telomeric region of chromosomes 1, 3, and 5 (Figs. 2, 3b). Similarly, the same pattern was observed for MQTLs distribution over the genetic consensus map. Chromosomes 2 and 4 possessed the majority of MQTLs in their non-telomeric intervals and chromosomes 1, 3, and 5 harbored several MQTLs in their sub-telomeric regions (Figs. 2, 3c). In contrast, the result of MQTLs distribution over A. thaliana physical map indicated a relatively high number of MQTLs in sub-telomeric regions (Fig. 3d).  www.nature.com/scientificreports/ The lowest number of QTLs and MQTLs were detected at the centromeric intervals for the ionomic traits (Figs. 2, 3). The MQTL-1/Chr2, MQTL-5/Chr3, and MQTL-5/Chr5 were located near the centromeric region of chromosomes 2, 3, and 5 of the physical and genetic consensus maps, respectively (Figs. 2, 3c,d). Furthermore, there was a significant correlation between the number of QTLs and MQTLs over the A. thaliana chromosomes (r = 0.97, P = 0.007) (Fig. 3a). The gene density over A. thaliana chromosomes showed a similar pattern with QTLs and MQTLs distributions (Fig. 4). Table 3. All the annotated genes located in each MQTL interval and the functions of the potential candidate are presented in Supplementary Table S1. The highest number of candidate genes in the identified meta-QTLs was observed in MQTL-2/Chr4 (774 genes) followed by MQTL-4/Chr2 and MQTL-2/Chr5 with 693 and 692 gene number, respectively ( Table 3). The lowest number of candidate genes belongs to MQTL-5/Chr3 with nine genes.

Discussion
The construction of a new consensus genetic map is a prerequisite for the projection of QTLs in the QTL database and running meta-QTL analysis. The results of the current MQTL analysis suggested that the MQTL-3/ Chr2, MQTL-7/Chr1, MQTL-1/Chr3, and MQTL-3/Chr5 with the highest number of initial QTLs identified in independent populations were the most viable, stable, and robust QTLs under different experimental conditions. Identification of meta-QTLs significantly increases the power and precision of our ability to map important traits that will provide greater precision for future fine mapping and marker development 48 . In the current MQTL analysis, the average confidence interval (CI) for QTLs of ionomic traits reduced up to eight folds compared to the mean CI of the original QTLs. The ability of the MQTL analysis methods to reduce the QTL confidence interval by taking advantage of pooling QTLs helps to increase the resolution of selected candidate genes 22 and consequently increases the precision of detecting functional candidate genes. Two complementary approaches for genetic mapping, linkage mapping (QTL mapping) and association mapping (LD mapping) have led to the  www.nature.com/scientificreports/ successful dissection of complex traits in many plant species by detecting QTLs and identifying marker-trait associations 49 . Often, the genetic and environmental background, number of traits, phenotypic variation, marker density, and map regulation are the most limiting factors in linkage mapping approaches. However, the meta-QTL analysis helps to identify the most stable QTLs regardless of limitations in traditional QTL mapping studies 49 .
The LD mapping uncovers genomic regions underlying quantitative traits with higher accuracy compared to the linkage mapping that is particularly efficient in species with low linkage disequilibrium (LD) 50,51 . Besides, the population structure of the studied panel in LD mapping can lead to false-positive association discovery 52,53 .
Despite the limitations of QTL-and LD-mapping methods, the statistical power to detect QTLs through the use of these approaches is particularly affected by the sample size 54 . To overcome these limitations, meta-analysis has been used to combine results from multiple studies and increase the power of genetic mapping 49 . The QTL meta-analysis method is a horizontal way to integrate information for the same trait from different experiments and populations. Thus, a meta-analysis allows for comparative genomics and is a valuable tool to complement information obtained from vertical integration 55 .
Our meta-QTL analysis over the A. thaliana chromosomes revealed that chromosomes 4 and 2 possessed a low number of QTLs and MQTLs, which could be attributed to fewer markers on these chromosomes. This result was in line with the results of the Serin et al. 56 study. Our results revealed a similar pattern for the distribution of QTLs and MQTLs over the A. thaliana chromosomes based on the genetic map used. Besides, the projection of QTLs and MQTLs across A. thaliana genome was in accordance with the constructed gene density map. It has been shown that QTLs result from the genetic segregation of sequence polymorphisms at functional elements such as regulatory sequences upstream of genes and/or coding sequences 57,58 . Therefore, it is expected that QTL density is related to the gene density, polymorphism rate at functional sites in genic regions and the frequency of recombination in the genome 32 . These outcomes showed for the first time that majority of QTLs and MQTLs controlling ionomic traits are located in the non-telomeric regions of the A. thaliana genetic map. However, higher QTL density at the sub-telomeric regions was apparent when physical distances were considered for meta-QTL analysis (Figs. 2, 3) 32 . Particularly, the distribution of MQTLs was higher in the gene dense regions located in non-telomeric and sub-telomeric regions in the A. thaliana genome. This result was in accordance with the result of the Martinez et al. 32 study showing QTLs appeared to more densely map on non-telomeric regions of genetic maps and sub-telomeric regions of A. thaliana physical map.
Identifying genes for ionomic traits control is a fundamental step for better understanding of ionome transportation, accumulation and mineral element changes in plants. Results of our meta-QTL analysis for ionomic traits showed that several important genes identified in the MQTLs were associated with micronutrient homeostasis. A better understanding of responsible genes involved in micronutrient homeostasis in plants is among the most important proceedings for breeding programs to combat "micronutrient malnutrition" 6 . A useful gene can be used to target multiple crops through genomic engineering approaches. The bHLH-encoding genes family including bHLH34 (8283003-8285081 bp) was identified in the MQTL-4/Chr3 of our meta-QTL study. Several members of bHLH-encoding genes are associated with the Fe homeostasis control 59 . The bHLH34 belonging to the clade IVc involved in A. thaliana Fe deficiency response 59 . The MQTL-2/Chr2 contains the NRAMP3 (9856149-9858778 bp) gene that is associated with iron deficiency in both roots and aerial parts of A. thaliana and localize to the vacuolar membrane, indicating its contribution to intracellular metal homeostasis [60][61][62] . The MQTL-3/Chr2 harbored another important gene, FIT (12004658-12006276 bp), which plays a predominant role in the regulation of iron mobilization [63][64][65][66] . The role of FIT in Fe acquisition had been documented 67 . The ZIP4 (19612020-19615383 bp) and IRT3 (22445299-22447299 bp) genes that located in the MQTL-7/Chr5 and MQTL-7/Chr1 intervals of our study are involved in Zn, Fe, Mn, and Cd transport 68 . Further, ZIP proteins that express in shoot tissues contribute to Zn uptake from the soil translocation of Zn throughout the plant 69,70 . In A. thaliana, the ZIP4 is highly expressed under zinc deficiency conditions 68,71 . The expression of AtIRT3 can respond to Zn 72 , which implies a possible role for IRT3 in Zn transport and homeostasis 73 . Overexpressing of AtIRT3 in A. thaliana increases accumulation of Zn in the shoot and Fe in the root of transgenic lines suggesting the role of IRT3 as a Zn and Fe-uptake transporter 73 . Detection of the MTP11 (16471560-16474060 bp) gene in the MQTL-5/Chr2 intervals in the current meta-analysis provides an evidence for Mn 2+ -specific transport activity of AtMTP11 that implicate the pre-vacuolar compartments in both Mn 2+ tolerance and Mn 2+ homeostasis The results of our meta-QTL analysis suggested that QTLs for K, Zn, P, Mn, Mg, and Ca were the most frequent distributed QTLs in the detected MQTLs regions that could be due to high heritability and the role of these ions in physiological pathways [80][81][82][83] . Our MQTL-analysis revealed co-localization between QTLs for Zn, Mn, and Fe ions. The co-localization of QTLs may reflect the correlations observed between the ionomic traits. Based on A. thaliana QTL analysis, different ionomic traits demonstrated significant correlation which most robust was the correlations between Zn, Mn, and Fe traits 16,42 . The result of co-localization frequency of ionomic QTLs on detected MQTLs regions in this study revealed the possible signs of MQTLs interaction, which impacts the association of ionomic traits at the genomic level. Co-localized QTLs or pleiotropic genes are of great significance for breeders and geneticists, for assessing the effects of selective genetic improvement programs, or for understanding pathways of genetic activity 84 . For biofortification, especially with Zn and Fe, large-effect and consistent QTLs are of particular importance and merit further examination 85 . Furthermore, the co-localization of QTLs of different traits in the same chromosomal regions suggests the existence of physiological and/or genetic relationships between traits 86 . Thus, the co-localized QTLs are important for breeding programs, QTL pyramiding, and simultaneous improvement of several ionomic traits.
Association analysis helps to unravel the genetic control of quantitative complex traits and validation of their corresponding genes in independent backgrounds 87,88 . In our study, four MQTLs (MQTL-3/Chr3, MQTL-4/ Chr3, MQTL-3/Chr5, and MQTL-3/Chr4) were collinear with nine highly associated SNPs responsible for Zn, Fe, and Mn quantitative traits. Further, the ITPK3 (5163220-5167192 bp) gene, which involved in micronutrients chelating and accumulation, was among the identified functional candidate genes detected through the current MQTL analyses and the GWAS study 47,89 . The ability of the MQTL analysis method to validate identified loci/ QTLs using GWAS have been previously documented [27][28][29] . Results of our meta-QTL analysis represent the coherence and consistency of MQTL analyses and GWAS for the identification of genomic regions corresponding to the studied ionomic traits.

Conclusions
To our knowledge, this is the first MQTL analysis report that identified several major genomic regions associated with ionomic traits in A. thaliana. This analysis defines a genome-wide landscape on the most stable genomic regions (MQTL-3/Chr2, MQTL-7/Chr1, MQTL-1/Chr3, and MQTL-3/Chr5) along with reliable genetic markers (CH.192L-Col, HH.357L, AD.191L-Col, SNP135, SNP105, DF.252L, m291, and BH.127L) that provide a robust tool for breeding ionomic traits through marker-assisted selection. Measuring trace elements in large sample size population is costly for breeders and identifying linked DNA markers accelerates pre-breeding and variety screening for biofortification programs and improving quality traits in plants. Results of our study suggested that the genomic positions of QTLs and MQTLs were mainly located in non-telomeric regions based on the newly constructed genetic consensus map, whilst the physical map revealed the projection of higher MQTLs in sub-telomeric regions. This was in line with the results of gene density in sub-telomeric and sub-centromeric regions of the Arabidopsis chromosomes. The results of the analysis for co-localization frequency revealed co-map QTLs controlling ionomic traits suggesting the existence of physiological and/or genetic relationships between these traits which lead to possible improvement of multiple ionome traits simultaneously. Further, the current study demonstrates that meta-QTL analysis method is cheaper than, and as informative as GWAS and then is powerful enough to provide greater resolution for future fine mapping without GWAS. Overall, the identified candidate genes at the detected MQTLs regions will provide a better understanding of ionomic variation in the A. thaliana genome. Besides their function can be generalized to other plant species which provides the raw material for breeders in various breeding programs such as biofortification.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.