Genome-wide association studies of seed metabolites identify loci controlling specialized metabolites in Arabidopsis thaliana

Plants synthesize specialized metabolites to facilitate environmental and ecological interactions. During evolution, plants diversified in their potential to synthesize these metabolites. Quantitative differences in metabolite levels of natural Arabidopsis thaliana accessions can be employed to unravel the genetic basis for metabolic traits using genome-wide association studies (GWAS). Here, we performed metabolic GWAS (mGWAS) on seeds of a panel of 315 A. thaliana natural accessions, including the reference genotypes C24 and Col-0, for polar and semi-polar seed metabolites using untargeted ultra-performance liquid chromatography-mass spectrometry. As a complementary approach, we performed quantitative trait locus (QTL) mapping of near-isogenic introgression lines between C24 and Col-0 for specific seed specialized metabolites. Besides common QTL between seeds and leaves, GWAS revealed seed-specific QTL for specialized metabolites indicating differences in the genetic architecture of seeds and leaves. In seeds, aliphatic methylsulfinylalkyl and methylthioalkyl glucosinolates associated with the GS-ALK and GS-OHP locus on chromosome 4 containing alkenyl hydroxyalkyl producing 2 (AOP2) and 3 (AOP3) and/or with the GS-ELONG locus on chromosome 5 containing methylthioalkyl malate synthase (MAM1) and MAM3. We detected two unknown sulfur-containing compounds that were also mapped to these loci. In GWAS, some of the annotated flavonoids (kaempferol 3-O-rhamnoside-7-O-rhamnoside, quercetin 3-O-rhamnoside-7-O-rhamnoside) were mapped to transparent testa 7 (AT5G07990), encoding a cytochrome P450 75B1 monooxygenase. Three additional mass signals corresponding to quercetin-containing flavonols were mapped to UGT78D2 (AT5G17050). The association of the loci and associating metabolic features were functionally verified in knockdown mutant lines. By performing GWAS and QTL mapping, we were able to leverage variation of natural populations and parental lines to study seed specialized metabolism. The GWAS data set generated here is a high-quality resource that can be interrogated in further studies.


Introduction
Two main phenotypic novelties have been critical during the transition from an aquatic to a terrestrial environment. The first of these innovations was the emergence of phenylpropanoid and lignin biosynthesis, allowing early terrestrial plants to acquire a relatively rigid body structure and colonize the land (Weng and Chapple, 2010). The second innovation consisted in the development of structures specialized for reproduction and dispersal, like pollen and seeds. These were essential for long-distance transport and successful colonization of the new environment by the offspring of primordial land plants (Linkies et al., 2010;Willis et al., 2014). Seeds, as a reproductive structure, also needed to be protected from adverse environmental conditions, including fungal attacks, insect feeding, or UV radiation. The chemical composition of seeds was thus selected not only to provide the essential nutrients during germination, but also to accumulate a number of specialized metabolites conferring protective properties against biotic and abiotic stresses (Debeaujon et al., 2000).
Arabidopsis thaliana is an ideal model to study the link between phenotypic and genomic variation, given the wealth of genomic resources available (Alonso-Blanco et al., 2016;Togninalli et al., 2018). The considerable genetic variation of Arabidopsis was employed to study local adaptation in collections of natural accessions (Seren et al., 2017). GWAS is a technique to leverage natural variation and was used in previous studies to detect adaptive traits (Atwell et al., 2010;Togninalli et al., 2018). GWAS assesses the effect of each genomic marker at the population level, represented by information on high-density SNPs, on a quantitatively-assessed phenotype with the likelihood of the association (Seren et al., 2017).
In A. thaliana, two major classes of specialized metabolites have been considered to confer protective properties to abiotic stress, namely flavonoids and glucosinolates. Flavonoids are arguably the best-characterized class of specialized metabolites that are universally distributed in the plant kingdom (Winkel-Shirley, 2001;Winkel-Shirley, 2002;Falcone Ferreyra et al., 2012;Tohge et al., 2017). By analyzing flavonoid-less mutants (ban, ap2, and transparent testa), Debeaujon et al. (2000) could show that a lack of flavonoids resulted in lower dormancy and structural aberrations in seeds (missing layers, modified epidermal layers). Tepfer et al. (2012) and Tepfer and Leach (2017) showed that flavonoid-less seed mutants exhibited lower survival rate when exposed to solar UV and cosmic radiation for 1.5 years. Generally, besides being involved in developmental and photoprotective processes, flavonoids convey antioxidative properties (Seyoum et al., 2006;Mierziak et al., 2014) and play a role in biotic stress (Treutter, 2005;Mierziak et al., 2014); as yet, however, no such information is available if the same is true in seeds. In A. thaliana seeds, a wide array of flavonoids can be found, mainly belonging to the subclasses of flavonols (mono-and diglycosylated quercetin, kaempferol, and isorhamnetin derivatives) and of flavan-3-ols (epicatechin monomers and procyanidin polymers, Routaboul et al., 2006).
Previous studies in our group revealed differences in seed glucosinolate levels of A. thaliana Col-0, C24 (unpublished data) in introgression lines (Törjék et al., 2008). Taken together with previous findings that showed differences in the accumulation of seed specialized metabolites in Arabidopsis ecotypes, we conducted an untargeted metabolic profiling analysis by UPLC coupled to high-resolution mass spectrometry (MS) on A. thaliana seed polar and semi-polar metabolites (covering several classes of specialized metabolites) to reveal quantitative differences of metabolites between accessions. To find putative novel gene candidates that control the accumulation of specialized metabolites, we conducted GWAS and, in a complementary approach, QTL mapping on the Arabidopsis IL population obtained from the cross between C24 and Col-0. We show here that (i) previously characterized metabolites (flavonoids and glucosinolates) associate with known loci, (ii) two unknown sulfur-containing metabolites map to glucosinolates-associated loci, and (iii) that the respective Arabidopsis SALK knockdown lines of the gene AT5G17050, previously selected from GWAS, showed quantitative changes in the levels of the associated quercetin-containing flavonol compounds.

Results and Discussion
Genome-wide association studies of untargeted seed metabolite analysis shows a large set of mass feature pairs associated with the same loci Genetic natural variation is an indispensable resource to find genes that are involved in the biosynthesis and regulation of plant specialized metabolites (Matsuda et al., 2015;Chen et al., 2016). Here, we determined the relative levels of polar and semipolar seed metabolites from about 300 A. thaliana ecotypes using UPLC-MS from two growing seasons (replicate 1 and replicate 2) and from one previously published set of leaf metabolites (Wu et al., 2018), and mapped the features to their associated genomic loci using the same GWAS approach we applied previously (Zhu et al., 2022). This approach encompasses mixed linear models to account for the amount of phenotypic covariance caused by the genetic relatedness, which should reduce confounding effects due to the population structure and kinship (Yu et al., 2006;Kang et al., 2008;Zhang et al., 2010;Vilhjálmsson and Nordborg, 2013). Due to computational constraints, we did not identify epistatic interactions, even though these will contribute to the observed phenotypes (Marchini et al., 2005;Cordell, 2009;Kam-Thong et al., 2011;Chen et al., 2014;Dong et al., 2015;Kerwin et al., 2015). Epistasis is the interaction of genetic variation at multiple loci that results in non-additive effects in the analyzed phenotypes (Soltis and Kliebenstein, 2015). In Arabidopsis, epistatic interactions typically involve the interaction of three or more loci (Wentzell et al., 2007;Rowe et al., 2008;Joseph et al., 2013a;Joseph et al., 2013b;Kerwin et al., 2015).
To compare metabolites across the different sets, we matched the alignments of mass features of seed replicate 1 and seed replicate 2/leaf based on their m/z deviation, retention time deviation, and the covariance between the two seed replicates resulting in 21007 features for the negative and 36194 features for the positive ionization mode. To further refine the accuracy, we imposed stricter matching rules, adjusting for retention time shift between replicate 1 and replicate 2 and a correlation of > 0.1, resulting in a total number of 9008 features for the negative and 12133 features for positive ionization mode (core set). 2882 In a next step, we created a table with the GWAS results of the two biological replicates of seeds and of the leaf samples. We reported for each joint mass feature the assigned QTL and LOD scores. From the core set, the two replicates from seed GWAS showed high positive that shared loci between the seed replicates showed a higher number (4884 for negative and 5688 for positive ionization mode) compared to that between seed and leaf GWAS (2242/1332 for negative and 2935/1942 for positive ionization mode). This may indicate that the reproducibility between the seed replicates of the core set is higher as when compared to the results from leaf GWAS. Alternatively, this may reflect a degree of variation in the genetic architecture. The comparison of the seeds and leaf datasets allowed the identification of tissue-specific QTL (Wu et al., 2018) and highlighted the different genetic architecture of these two tissues in controlling the accumulaton of some specific mass features. However, there were 1026 and 1247 loci controlling the mass features in the core sets that are shared between the two seed replicates and the leaf data set, indicating conserved loci controlling the levels of mass features across different tissues. The mass feature pairs that showed overlap between the two seed replicates, and for some of these also to the leaf data set, represents a highly valuable resource that we make fully available in the Supplemental Material. In the subsequent paragraphs we only focus on those associations related to glucosinolates and flavanols, as well as on some unknown mass features supposedly representing novel glucosinolates and flavonoids.
Alongside the shared loci, the majority of loci were not shared between the different mass feature sets (Figure 1 C and D). When analyzing the distribution of the proportion of mass features mapped to loci, the sets that did not show intersection (seed replicate 1, seed replicate 2, leaf) had a higher proportion of two or more of mass features mapped to the same loci compared to sets that show intersection (Supplementary Figure S3). This could be attributed to measurement errors, associations of non-causative markers with a given trait, driven by linkage with causative markers (Korte and Farlow, 2013), reflect environmental variance, and/or genotype-environment interaction effects. We assume that the significant associations from the intersection sets (Figure 1 C and D), being conserved between the different replicates, represent genuine QTL characterized by lower sources of error. Intersection sets of QTL for mass feature pairs with LOD > 5.3 for negative (C) and positive ionization mode (D). LOD: logarithm of odds; QTL: quantitative trait loci.
As a complementation, we performed QTL mapping using biparental NILs from Col-0 and C24 (Törjék et al., 2008). This population was useful in detecting saiginols (phenylacylated flavonols) in floral tissues (Tohge et al., 2016). In the GWAS population and NIL population of Arabidopsis seeds, saiginols were not detected.

Variation in glucosinolate levels in seeds is controlled by the GS-ELONG, GS-ALK, and GS-OHP loci
Based on previous studies, we annotated several metabolites in our data set, including amino acids, glucosinolates, and metabolites from the flavonoid and phenylpropanoid biosynthetic pathways (Supplementary Table S5 Figure   S4 B for 3-hydroxypropyl glucosinolate, GWAS population) indicating that the allelic variation at these target loci is responsible for the observed metabolite differences. Indeed, some of the SNPs for these genes involved in glucosinolate biosynthesis resulted in amino acid differences (Supplementary Table S9). The LD analysis for 3-hydroxypropyl glucosinolate revealed that the alleles on chromosome 4 are in high LD (standardized LD r 2 close to 1) spanning the genomic region containing AOP1, AOP2, and AOP3 (Figure 2 C, left panel). In A. thaliana, LD usually decays 50% within 5 kb (Gan et al., 2011;Korte and Farlow, 2013). Here, the loci containing AOP1, AOP2, and AOP3 showed wider LD. The situation on chromosome 5 marks a sharp decrease for 3-hydroxypropyl glucosinolate and peaks in the gene region of MAM1 (AT5G23010). Interestingly, MAM3 was in high LD (r 2 > 0.6) with the SNP showing the highest LOD in MAM1, but did not show as high r 2 values as neighboring genes within the locus on chromosome 4. This indicates that MAM1 is the main gene controlling 3-hydroxypropyl glucosinolate levels. Previously, these loci were also detected from GWAS of glucosinolate levels in leaves (Chan et al., 2011). The Arabidopsis gtr1 gtr2 double mutant, which lacks (or contains low amounts, depending here on the type of the mutations it carries) the nitrate/peptide transporters responsible for glucosinolate transport to seeds, did not accumulate glucosinolates in seeds and exhibited a tenfold over-accumulation in the source tissues leaves and silique walls (Nour-Eldin et al., 2012). Thus, it seems that the variation in glucosinolate levels is `inherited' from these source tissues.
Figure 2: Genome-wide association mapping for 3-hydroxypropyl glucosinolate (negative ionization mode). A: The Manhattan plot of 3-hydroxypropyl glucosinolate shows two peaks in each replicate on chromosomes 4 (highest LOD: 11.01) and 5 (16.35). These loci contain the genes AOP1, AOP2, and AOP3 (chromosome 4), MAM1 and MAM3 (chromosome 5) that are involved in glucosinolate biosynthesis. Only replicate 1 is shown here. B: Haplotype analysis of metabolite levels of 3-hydroxypropyl glucosinolate. The nucleotide sequence differences were statistically associated with the levels of 3hydroxypropyl glucosinolate (ANOVA q-value: 1.78e-20 for replicate 1). Only data for replicate 1 is shown in A and B. The data for replicate 2 is depicted in Supplementary Figure S4. C: LD analysis of the associated genomic regions on chromosome 4 and 5 for 3-hydroxypropyl glucosinolate. The locus on chromosome 4 shows LD for the genomic region containing the genes AOP1, AOP2, and AOP3, while the locus on chromosome 5 marks a sharper decrease in standardized LD (r 2 ) indicating that the MAM1 gene is mainly responsible for the natural diversity of 3-hydroxypropyl glucosinolate levels. AOP: alkenyl hydroxyalkyl producing; A.U.: arbitrary units; LD: linkage disequilibrium; LOD: logarithm of odds; MAM: methylthioalkylmalate synthase.
Unknown sulfur-containing metabolites map to GS-ELONG, GS-ALK, and/or GS-OHP loci in genome-wide association mapping  Table S4) also showed associations to these loci (Supplementary Figure S7).
The LD analysis revealed that for the unknown 626 the SNP with the highest LOD was located near or within the MAM1 gene. The standardized LD, r 2 , decreased sharply when moving away from the MAM1 gene (Supplementary Figure S6). To reveal the chemical composition of the two unknowns, we fed isotope-labeled 13 C and 34 S to the siliques and analyzed the metabolites by LC-QTOF-MS. The unknown 596 (m/z 596.1104 in negative ionization mode) and 626 (m/z 626.1032 in negative ionization mode) contain most probably 20 C atoms and 22 C atoms, respectively, based on isotope feeding experiments with 13 C. The MS analysis indicated for the feeding experiments with 34 S that the two unknown compounds contain two S atoms (Supplementary Figure S8). Interestingly, the QTL mapping between C24 and Col-0 introgression lines of the unknowns 596 and 626 identified an additional locus close, but not overlapping, to GS-ALK and GS-OHP (AT4G15733-AT4G24620, M502 in Supplementary Figure S5). From the GWAS analysis here, the candidate region for the unknown 596 is  Table S10) and SALK lines were individually cultivated in two replicates.
Seeds of the SALK lines were analyzed by UPLC-MS (Supplementary Figure S10). The resulting data set was analyzed in terms of presence and differential abundance of the associated mass features with respect to Col-0 and other SALK line seeds in negative and positive ionization mode. Only the line SALK_049338 (AT5G17050, encoding UDP-GLUCOSYL TRANSFERASE 78D2) showed differential abundance for several mass features compared to the control lines (Supplementary Table S11 and S12).
When growing the line SALK_049338 (AT5G17050), we observed shorter stature for all plants as compared to the wild type (Figure 3 A and Supplementary Figure S11), a phenotype also previously reported when mutating this gene (Yin et al., 2014 (Figure 3 B). Here, we focused on the analysis of the association involving candidate structural genes. In this paper we focused on candidate structural genes of the glucosinolate/flavonoid pathways, although we report in the Supplementary Material the full list of significant associations that may represent a resource to investigate the additional control these metabolites may have at the level of pathway regulation. Furthermore, the results from glucosinolates and some of the flavonoids indicated pleiotropic effects and collocating QTL for joint mass features. This analysis can be extended to a wider scale and to non-biosynthetic enzymes. Moreover, the core set exhibited differences in QTL between the mass features from the two seed replicates and the leaf replicate. Future studies will investigate the variation in the genetic architecture of traits controlling the levels of specialized metabolites across different tissues.

Conclusion
Here, we performed GWAS on metabolic mass features of two biologically independent replicates of seed from two growing seasons and one replicate of leaves obtained by untargeted UPLC-MS. As a complementary approach, we performed QTL mapping of NIL introgression lines between C24 and Col-0 for specific seed metabolites. By including GWAS of leaf metabolites, we detected 4884 and 5688 loci for mass feature pairs (negative and positive ionization mode) that were exclusively detected for seed GWAS, indicating differences in tissue-specific associations between seeds and leaves. On the other hand, 1026 and 1247 QTL for mass feature pairs (negative and positive ionization mode) were conserved across seed and leaf tissues in GWAS.
In seeds, aliphatic methylsulfinylalkyl and methylthioalkyl glucosinolates as well as two show changes in flavonol levels. We would like to draw the following conclusions regarding the genetic architecture of seed specialized metabolism: (i) seed specialized metabolism differs substantially from leaf metabolism as shown by the identification of QTL that differ between these tissues, but the two tissues also exhibit common genetics to some degree; (ii) AOP and MAM genes are key regulators for glucosinolate seed metabolite levels in seeds.
Aliphatic glucosinolates are presumably not synthesized in situ in seeds, but are transported from source tissues to seeds. The variation of aliphatic glucosinolates is 'inherited' from these source tissues; (iii) the alleles of UGT78D2 (AT5G17050) affect the levels of quercetincontaining flavonols in seeds. The natural GWAS population was shaped by processes of genetic adaptation and meiotic events during evolution. This results in greater phenotypic variance compared to the NIL population between Col-0 and C24 as exemplified by differences in flavonoid levels. However, the overlap suggests, as previously stated (Brog et al., 2019), that genome-wide association and QTL mapping are complementary techniques to study seed specialized metabolism.

Plant material
The for 10 min at room temperature (17,900 g), and the supernatant was transferred to a new tube.
The tubes were centrifuged for 10 min at room temperature (17,900 g). 135 μL of the supernatant were transferred to a new tube, dried by speed-vac for 2-3 h, filled with argon, and stored at -80°C. On the day of analyses, the samples were resuspended in 100 μL 80% MeOH and transferred to sample vials.
Metabolites from leaves were extracted from 50 mg leaf material (cooled in liquid N2) using 500 μL of the same extraction buffer as above. The same extraction protocol was followed as  To test for differences in SALK lines, the log2-normalized raw intensities were tested against Col-0 or the respective complement of SALK lines using limma (v3.50.3). To this end, linear models were fitted for each metabolic feature using lmFit and moderated t-statistics were computed by empirical Bayes moderation of the standard errors towards a global value using eBayes (trend = TRUE). p-values were adjusted using FDR via the Benjamini-Hochberg method. Since there was no corresponding second replicate available in positive ionization mode, the corresponding features in negative ionization mode were determined using correlation analysis and retention time window thresholding. If multiple features in negative ionization mode matched to the feature in positive ionization mode, the feature with highest correlation to the feature of positive mode (replicate 1) was selected.

Notes
The authors declare no competing financial interest.