# Source paper archive - Source title: Resolving Gene-Metabolism Links through Multi-Metabolite mQTL Commonality Analysis - Author(s): Rosental, L., Zach, Y. Y., Lapidot-Cohen, T., Weil, S., Tabachnik, C., Sternlicht, E., Niger, T., Borovikova-Sheinker, A., Panda, S., Rogachev, I., Shahaf, N., Ouyang, B., Aharoni, A., Dong, S., Sade, N., Brotman, Y. - Source landing URL: https://www.biorxiv.org/content/10.64898/2026.11.03.763898v1 - PDF URL: https://www.biorxiv.org/content/10.64898/2026.11.03.763898v1.full.pdf - PDF SHA-256: ea5fc5632802330eec68da1b37be824192dc8305cf77cdc70a558f839c96e436 - Converter: pymupdf4llm 0.0.27 - Limits: Automated PDF text extraction; reading order, equations, tables and formatting may be incomplete or incorrect. Images are omitted, and no OCR was performed. Markup and links are neutralized for safe downloading. Consult the original PDF for the authoritative version. --- [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1 **Resolving Gene-Metabolism Links through Multi-Metabolite mQTL Commonality** 2 **Analysis** 3 4 L eah Rosental [1] *, Yair Yehoshua Zach [2], Taly Lapidot ‐ Cohen [1], Simy Weil [1], Claire 5 Tabachnik [1], Ezra Sternlicht [1,2], Tali Niger [1], Anna Borovikova-Sheinker [2], Sayantan 6 Panda [3], Ilana Rogachev [3], Nir Shahaf [3], Bo Ouyang [4] Asaph Aharoni [3] S huchao Dong [5], Nir 7 Sade [2], Yariv Brotman [1,2] 8 9 1 Department of Life Sciences, Ben Gurion University of the Negev, Beersheva, Israel; 2 10 School of Plant Sciences and Food Security, Institute for Cereal Crops Research, Tel 11 Aviv University, Tel Aviv, 69978 Israel; [3] Department of Plant and Environmental 12 Sciences, Weizmann Institute of Science, Rehovot, Israel; [4] National Key Laboratory for 13 Germplasm Innovation and Utilization of Horticultural Crops, Huazhong Agricultural 14 University, Wuhan 430070, China; [5] Institute of Vegetable Crop, Jiangsu Key Laboratory 15 for Horticultural Crops Genetic Improvement, Jiangsu Academy of Agricultural Sciences, 16 Nanjing 210014, China 17 18 *Corresponding author: Yariv Brotman. Department of Life Sciences, Ben Gurion 19 University of the Negev, Beersheva, Israel and School of Plant Sciences and Food 20 Security, Institute for Cereal Crops Research, Tel Aviv University, Tel Aviv, Israel. 21 Email: **b** **[rotman@tauex.tau.ac.il](mailto:brotman@tauex.tau.ac.il)** 22 *Corresponding author: Leah Rosental. Department of Life Sciences, Ben Gurion 23 University of the Negev, Beersheva, Israel. Email **leahrose@post.bgu.ac.il** 24 1 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 25 Abstract 26 High-throughput metabolomics detects so many compounds that most studies use 27 only a fraction of them, focusing on the few metabolites with an annotated structure or 28 substructure. Metabolite genome-wide association studies of specialized metabolism 29 detect many metabolite quantitative trait loci (mQTLs), but most associated metabolite 30 features remain unannotated, limiting gene discovery. In leaves of 362 tomato 31 accessions under control and water-deficit conditions, only 13 of the 826 metabolite 32 features associated with 77 reproducible mQTLs were annotated with good confidence. 33 To study the mQTLs with the largest impact, we developed commonality analysis, an 34 approach that does not depend on well-studied pathways. 35 Commonality analysis uses biologically defined groups, such as the features 36 associated with one mQTL, to infer relatedness between metabolites. By identifying the 37 fragments, neutral losses or mass transformations that are common in a group, it links 38 the group to a substructure, pathway or reaction, and can extend annotation within the 39 group. Most high-impact mQTLs had a significantly distinct fragment commonality, 40 whose outstanding fragments suggested the pathway involved and its candidate gene. 41 We demonstrate the approach in three representative mQTLs: (i) common fragments 42 pointed to a sinapoyl-hexose pathway and a serine carboxypeptidase-like candidate 43 acyltransferase; (ii) common fragments marked the features associated with a purine 44 permease as nitrogen-containing, which isotope labeling then supported; (iii) a common 45 2.02 Da transformation pointed to a double-bond reduction and a putative reductase. 46 Overexpression of each gene reproduced the commonality patterns, supporting the 47 candidate genes and showing that the approach also applies to experimentally defined 48 metabolite groups. 2 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 49 **Introduction** 50 Plants produce an extraordinary diversity of metabolites that support growth, 51 development, defense, and environmental adaptation. Tomato ( _Solanum lycopersicum_ ) 52 provides an excellent system for studying genetic variation and crop metabolism. 53 Extensive genomic resources and germplasm collections [1–4] have facilitated the genetic 54 dissection of traits related to domestication, fruit quality, flavor, and metabolite 55 composition [5–9] . However, most tomato metabolomic studies have focused on fruit 56 associated metabolites, whereas the genetic regulation of the leaf metabolome remains 57 largely unexplored [10], despite the central role of leaves in photosynthesis, productivity, and 58 environmental response [11,12] . 59 Using the natural variation of metabolite abundances as phenotypic traits in metabolite 60 genome-wide association studies (mGWAS) enables the identification of metabolite 61 quantitative trait loci (mQTLs) in the genome and facilitates the discovery of metabolic 62 genes, pathways, and regulatory networks [13–15] . Environmental perturbation expands the 63 repertoire of metabolites produced and consequently the power of mGWAS to reveal 64 mQTLs that are weak or absent under other conditions, as well as those stable in both [16–] 65 20 . However, moving from statistical association to biochemical interpretation usually 66 requires information on the identity of the metabolites involved. 67 Despite advances in metabolite detection by liquid chromatography-mass 68 spectrometry (LC-MS), most specialized metabolites remain unannotated [21,22] . Current 69 annotation approaches commonly match MS/MS (or MS [2] ) spectra to open-source 70 reference libraries such as MassBank and GNPS [23,24] or vendor-specific libraries such as 71 mzCloud and METLIN. While effective, direct library matching excludes the identification 72 of metabolites missing from reference libraries. To expand the annotation of metabolites, 73 multiple approaches integrate the metabolite information obtained by LC-MS/MS in 74 various ways. For example, Feature-Based Molecular Networking (FBMN) links 75 metabolites according to MS [2] spectral similarity, enabling propagation of substructure 76 annotations within related chemical families [24,25] . The SIRIUS software uses MS [2] spectra 77 to predict chemical fingerprints of molecular substructures and molecule structures [26] . The 3 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 78 MS2LDA algorithm identifies recurring fragment and neutral-loss motifs that facilitate 79 substructure annotation of unknown compounds [27] . Although the number of approaches 80 and _in silico_ tools for expanding metabolite annotation is increasing, their accuracy and 81 suitability for annotating _de novo_ metabolites and chemical families are still limited. For 82 example, a systematic evaluation of fragment-ion annotations in three major _in silico_ MS [2] 83 libraries (HMDB, METLIN, and mzCloud) revealed substantial inaccuracies, highlighting 84 the enduring limitations of current annotation approaches [28] . 85 Beyond spectral similarity, biochemical relationships between metabolites can provide 86 an additional source of information. Candidate substrate-product pair (CSPP) networks 87 identify metabolites separated by characteristic mass differences corresponding to 88 enzymatic transformations [29], which can be used in combination with molecular 89 networking to identify recurrent biochemical modifications across metabolomes [30] . Tools 90 such as MetNet [31], NetID [32] and MetDNA [33] use biologically relevant mass differences and 91 metabolic connectivity to improve metabolite interpretation. 92 In the context of mGWAS, several methods have begun to integrate biochemical 93 information directly into genetic analyses. For example, chemical-tag-based semi 94 annotated metabolomics demonstrated that recurring modifications can support 95 candidate-gene discovery and pathway reconstruction without requiring full metabolite 96 identification [34] . Shortest reactional distance (SRD) analysis, which is based on KEGG 97 reactions, incorporated metabolic pathway information to prioritize gene-metabolite 98 relationships beyond conventional significance thresholds [35] . Together, these studies 99 illustrate the growing effort to move from statistical associations toward mechanistic 100 interpretation of mGWAS datasets. Nevertheless, most studies still rely on annotated 101 structures to guide the selection of genes for functional characterization [36–38] . Considering 102 that only ~1% of plant specialized metabolites have been annotated, a recent perspective 103 argued that the next major advance in plant metabolomics will require overcoming the 104 strong bias toward the limited number of annotated metabolite families [39] . Further, these 105 methods do not systematically exploit the co-association of multiple metabolite features 106 as an information source in their analysis. 4 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 107 Here, we present a large-scale analysis of leaf metabolism in a genetically diverse 108 tomato population grown under control and prolonged water-deficit conditions. We 109 provide a resource comprising transcriptomes, hundreds of chromatograms, data 110 dependent tandem mass spectra (ddMS [2] ), curated metabolite annotations and mQTL 111 associations. 112 To confront the scale and complexity of the data and to study metabolic gene 113 regulation in an unbiased manner, we developed commonality analysis, which extracts 114 shared mass descriptors (fragment ions, neutral losses and transformations) from the 115 features associated with each mQTL. We show that most high-impact mQTLs, those with 116 enough features to statistically test the mass descriptors, have a significantly distinct 117 commonality, which helps identify candidate pathways, reactions, and genes underlying 118 metabolic variation. We present three examples of how this strategy improves metabolite 119 annotation and enables identification and functional characterization of genes underlying 120 mQTLs. 121 122 5 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 123 **Results** 124 To investigate the genetic regulation of tomato leaf metabolism, a diverse tomato 125 collection was grown in a greenhouse experiment. The "TS" collection includes 362 126 accessions with representatives from _S._ _pimpinellifolium,_ _S._ _lycopersicum_ var. 127 cerasiforme and diverse _S. lycopersicum_ tomato varieties including many heirloom, 128 transitional and modern cultivars [2,6] . Because water stress substantially alters leaf 129 metabolism [40–43], we evaluated the population under two conditions: typical agricultural 130 production conditions (control) and reduced water availability (water deficit; WD). 131 Comprehensive evaluations of gene expression and metabolism were integrated with 132 high-density genomic polymorphism data [2] to enable thorough investigation (Fig. 1A). 133 Gene expression and metabolic variation under two conditions 134 We first set out to evaluate the severity of the imposed stress and the changes it 135 caused. The applied WD caused a significant and substantial drop in leaf water content 136 (Fig. 1B), but allowed plant survival, albeit with strongly reduced vegetative productivity 137 (Fig. 1C) by the end of the growing season. Transcriptomic analysis quantified relative 138 expression of 11,191 genes in 539 samples (Supp. Table 1). As expected, a global shift 139 in the gene expression pattern between control and WD plants was found (Fig. 1D), 140 explaining 48% of the variation in the first principal component. Most genes differed 141 significantly but mildly (low population-wide fold change) between conditions. Still, some 142 genes changed significantly and substantially: 251 (2.2%) genes decreased (fold 143 change < 0.5) and 199 (1.7%) increased (fold change > 2) under WD (Fig. 1F). 144 Population-wide variance in gene expression was greater among WD samples, as can 145 be seen in the wider spread of samples along the first and second principal components 146 and in the increase in expression variation of individual genes within the population (Fig. 147 1E). Gene Ontology (GO) analysis revealed that genes increased in WD were enriched 148 in secondary metabolite biosynthetic process (GO:0044550) genes, whereas general 149 metabolic process (GO:0008152) and other central metabolic process related GO terms 150 were depleted (Supp. Figure 1; Supp. Table 2). 151 Leaf metabolomic profiling was conducted using an untargeted LC-MS method 152 which detects polar and semi-polar metabolites. Over 100,000 metabolite features 6 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 153 (peaks with a unique mass-to-charge ratio (m/z) and retention time (RT) after isotope 154 clustering, hereafter features) were detected in positive and negative ionization modes 155 combined. After applying quality control filters, 78,138 were retained for the following 156 analyses. 157 Growth under WD greatly altered the metabolic composition of the leaves (Fig. 1G). 158 Accordingly, the vast majority of features were significantly increased or decreased in 159 WD. Although the significant changes were often mild, over 3.5% of features changed 160 by more than 2-fold between the conditions (Fig. 1H). In addition, we detected 15,419 161 (19%) features exclusively in control conditions and 19,267 (24%) in WD, highlighting 162 the benefit of including multiple conditions for a comprehensive analysis (Fig. 1I). 163 Leaf metabolite annotations 164 One major bottleneck in metabolomics is assigning structural information to detected 165 features [44,45] . We therefore combined multiple annotation methods to maximize 166 annotation coverage (Supp. Table 3). 167 Thirty diverse samples from our GWAS experiment were chosen for comprehensive 168 feature ddMS [2] fragmentation, leading to 13,781 ddMS [2] spectra obtained in positive 169 ionization mode and 12,298 in negative mode. We constructed a feature-based 170 molecular network (FBMN) from the ddMS [2] spectra and matched them to the GNPS [46] 171 library (Supp. Figure 2). The network included 548 well-matched library annotations, but 172 many matched multiple features at different RTs. Next, we searched several public and 173 in-house libraries using Compound Discoverer (Thermo Fisher Scientific). High 174 confidence matches to the in-house run library (based on MSMLS, Mass Spectrometry 175 Metabolite Library of Standards, Sigma-Aldrich) identified 81 features as 64 unique 176 metabolites. External library searches matched 458 features to metabolite MS [2] spectra. 177 Additional comparisons to published LC-MS based annotations gave putative 178 annotations to 161 features, but at lower confidence because they were based on 179 fragments co-locating in MS [1] rather than on ddMS [2] matches. The same set of 30 180 samples was analyzed with WeizMass [47] using its optimized LC-MS method [9], yielding 181 high-confidence metabolite annotations for 122 features. Each method contributed 182 some new annotations, but as expected many features had redundant identifications 183 from multiple sources. Overall, we curated metabolite annotations for 1,442 features 7 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 184 (Supp. Table 3), with varying annotation confidence levels [48], which were later used to 185 interpret the mGWAS results. Figure 1: **Experiment overview and stress-induced variation** A. Overview of experiment setup, datasets produced and the main analysis and results to which they led. B. Leaf water content, as the ratio of leaf weight lost after full oven drying to the fresh leaf weight, in leaves collected three weeks after treatment onset from control (light blue) and water deficit (WD; orange) treated plants. N control =155, N WD =54. C. The fresh vegetative weight of whole plants measured after last fruit harvest. N control =656, N WD =589. D. Global effect of WD on leaf gene expression. PCA based on relative gene expression of leaves collected from control (blue) and water deficit (dark orange) treated plants. E. Distribution of variation (standard deviation) in the expression of each gene across the population in either growth condition. F. Gene expression in response to WD. Volcano plot displaying the significance of a t-test comparing transcript levels between water deficit and control conditions to the fold change (after log10 transformation) of the same comparison. Transcripts with significant increase (p-value <α=0.05, Bonferroni correction) and fold change exceeding double expression in water deficit (dark orange) or control (blue) are colored. G. Global effect of WD on metabolic profile. PCA based on metabolic profiles of leaves collected from control (light blue) and water deficit (orange) treated plants with metabolite features detected in both positive and negative ionization modes. H. Feature changes in response to WD. Volcano plot displaying the significance of a paired t-test comparing feature levels of each accession between WD and control conditions to the median fold change (after log2 transformation). Features with significant increase (p-value <α=0.05, Bonferroni correction) and fold change exceeding two-fold expression in water deficit (orange) or control (light blue) are colored. I. The number of features detected exclusively in control (blue circle), water deficit (orange circle) or in both (overlap). 186 8 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 187 Overview of mQTLs detected in control and in water deficit conditions 188 The abundance of each feature was used as a phenotype for GWAS [2], with two 189 independent repetitions for each condition (control and WD). Adjacent genomic 190 associations of the same feature within a specific window of 100 kb were collapsed and 191 defined as a single mQTL. Features associated with overlapping genomic regions were 192 assigned to the same mQTL and together form its mQTL feature group. 193 In total, 719 unique mQTLs were found (Fig. 2A; Supp. Table 4). In control 194 conditions 410 mQTLs were found, with 192 repeatedly confirmed in the second control 195 replicate. For WD 401 mQTLs were found, and 186 were confirmed by repeated 196 association in WD. Many mQTLs were detected only once. These tended to have very 197 few associated features, reflecting the risk of random sporadic associations in such 198 high-throughput analysis. We therefore considered only mQTLs for which the feature 199 association was reproduced in the second repetition. We found that 36 mQTLs were 200 constitutive with consistent feature associations in both repetitions of the two conditions, 201 while 13 or 16 mQTLs had features repeatedly associated only in WD or in control, 202 respectively. Interestingly, 12 mQTLs had features repeatedly associated in WD and 203 other features repeatedly associated in control. This could be caused by mistaken co 204 association because of closely linked polymorphisms or by activity on features present 205 exclusively in either condition. Together, we detected 77 reproducible mQTLs. Again, 206 the condition-specific mQTLs show the benefit of investigating different conditions for 207 obtaining a more complete picture of genetic control of leaf metabolism. 208 Across all mQTLs, 826 features had repeated genomic associations, often with 209 multiple features co-associating with the same mQTL. Of all features with any significant 210 genomic associations, we were able to obtain ddMS [2] spectra for 1,445 (55%) in 211 positive, and 750 (48%) in negative ionization mode. However, only 13 features in 212 mQTL feature groups were annotated with good confidence. 213 Our approach was to be as inclusive as possible, to make maximal use of current 214 high-throughput detection methods. In this way, the results are not biased toward well 215 studied pathways and may reveal insights into previously uncharacterized compound 216 families. We see value in publishing detailed searchable genomic associations of 217 unannotated features, so that the scientific community can use them as knowledge of 9 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 218 compound structures expands. Indeed, we were able to obtain and analyze large 219 amounts of data (Fig. 1A). However, this multitude of data required new approaches to 220 highlight and prioritize relevant information for study. We therefore set out to develop 221 methods for functional characterization of feature groups based on their genomic 222 associations, alongside functional gene characterization based on associated feature 223 groups. Figure 2: **mQTL overview and scenarios** A. Each detected mQTL is represented as a column placed according to physical location on the tomato genome (x axis). The width of the column ranges from the first to the last significant SNP in the mQTL. The height of the column represents the number of metabolite features in the mQTL feature group, scale noted on the left y axis. The dot accompanying each mQTL column represents the significance of the highest SNP-feature association in the mQTL, scale shown on the right y axis. Results of mapping in control conditions are in the top panel in blue, and from water deficit in mirror image in the bottom panel in orange. The mQTL and genes studied in the following sections are marked by arrows. B. Graphical representation of the scenarios discussed in the text for mQTLs with large mQTL feature groups. Structures and colors represent arbitrary molecular substructures. Inset plots are simplified commonality spectra expected from various scenarios. In scenario IV, the pleiotropic effects are not expected to manifest in a typical commonality spectrum, so it is not drawn. 10 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 224 225 Commonality analysis: concept and method development 226 QTLs are usually caused by polymorphisms linked to one causal polymorphism, 227 which primarily influences one causal gene [5,13,15] . Although not always the case, we 228 assumed one causal gene for each mQTL in our working model. When multiple features 229 co-associate with one mQTL, this could arise from various scenarios: I) changes to an 230 enzyme in a pathway which influences upstream and downstream metabolites. II) 231 changes to an enzyme within an interconnected metabolic network. III) changes to an 232 enzyme with multiple substrates. IV) changes to a regulator which indirectly affects 233 metabolism, such as a transcription factor. These various scenarios are expected to 234 produce distinct patterns among the features of the mQTL feature group (Fig. 2B). In 235 detail: 236 Scenario I) The mQTL-influencing gene is an enzyme which catalyzes a reaction 237 within a defined, linear or weakly branched specialized metabolism pathway. In such a 238 case we expect to see altered abundance of all the compounds in the pathway, which 239 will then be co-associated with the mQTL. For example, a loss-of-function or reduced 240 activity allele would be associated with an accumulation of precursors and reduction in 241 downstream products. Consecutive reactions in specialized metabolism pathways 242 usually add/remove/alter various parts of the compounds, but the backbone structure is 243 usually preserved and would therefore manifest in the same fragment appearing in the 244 ddMS2 spectra of the features in the mQTL feature group (Fig. 2B I). Identifying these 245 highly shared fragments would define the central substructures and help focus 246 annotation efforts on the basal pathway influenced by the mQTL. 247 Scenario II) The gene is an enzyme involved in primary or interconnected and 248 redundant metabolic pathways. As in scenario I, altered activity in a specific reaction will 249 ripple and cause coordinated changes in upstream and downstream metabolites. 250 However, redundancy may buffer the allele influence and cause multiple low effect 251 mQTL associations, as often seen in mGWAS of primary metabolites [37,49] . Still, if the 252 polymorphism has a very strong effect, we may see multiple feature-mQTL 253 associations. Due to interconnectivity, different groups of metabolites may be affected. 254 In this scenario shared mass descriptors may be found for features with tight 11 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 255 associations to the polymorphism, but we may expect several common descriptors 256 shared by a lower proportion of the features in the mQTL feature group (Fig. 2B II). 257 Scenario III) The gene is an enzyme which catalyzes a reaction on similar 258 compounds involved in several pathways. In this case we would see changes in 259 compounds of various families. Even promiscuous enzymes usually specialize in a 260 specific type of reaction that is often associated with a specific substructure the enzyme 261 can recognize and bind [50,51] . In this case the recognized substructure may result in a 262 fragment shared by many of the features in the mQTL feature group. 263 Moreover, the specific reaction would add/subtract a substructure with a specific 264 molecular mass which would manifest in a specific m/z difference between CSPP [29,31] . 265 For example, acetylation or deacetylation (C 2 H 2 O/H) would add or subtract m/z=42.01, 266 respectively. Multiple CSPP features would be co-associated with the mQTL, and the 267 reaction-causing m/z difference would be frequent among all calculated paired m/z 268 differences, named transformations hereafter (Fig. 2B III). A transformation is a 269 calculated mass difference and does not necessarily demonstrate a biochemical 270 conversion. 271 Scenario IV) If the affected gene is a regulator (for example, a transcription factor) it 272 may affect multiple metabolic pathways, usually ones involved together in a coordinated 273 developmental program or response to changing conditions. They may, however, 274 involve very different reactions and compound structures (Fig. 2B IV). Such regulators 275 may be identified by their central role in gene-gene or gene-metabolite correlation 276 networks but are beyond the scope of this paper. 277 We hypothesize that mQTL feature groups will have a significantly distinct 278 commonality of descriptors, either fragments, neutral losses (NL) or transformations, 279 compared to random groups of features. Further, these distinct commonly shared 280 descriptors will assist in identifying the mQTL causal gene, reaction and pathway. We 281 set out to test this hypothesis by grouping ddMS [2] spectra of the features in each mQTL 282 feature group and identifying the frequency of descriptors in the group. We visualized 283 the commonality, the set of descriptor frequencies in a group, as commonality spectra: 284 plots that display, for each m/z value (x axis), the frequency of that descriptor, i.e. the 285 proportion of features in the group that carry it (y axis). Because of the different nature 12 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 286 of features and fragments detected in positive or negative ionization mode, the analysis 287 and representation were done separately for each ionization mode. Next, we developed 288 a statistical framework to define the significance of mQTL-based commonality 289 compared to ungrouped features (Supp. Note 1: development of the statistical 290 framework for commonality analysis). 291 13 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 292 Descriptor frequencies in mQTL feature groups differ from those in all leaf features 293 _Fragments_ 294 Prior to studying the commonality of mQTL feature groups, we identified the most 295 common fragments in all collected ddMS [2] spectra (Supp. Table 5). In positive mode the 296 most common fragments 85.03, 92.20, 127.04, 97.03 were found in 22%, 21%, 18%, 297 17% of all features with ddMS [2] spectra, respectively (Fig. 3A). In negative mode the 298 most common fragments 59.01, 101.02, 89.02, 71.01, 113.02 are hexose-derived 299 fragments (Fig. 3B). They were found in over 25% of all features with ddMS [2] spectra, 300 showing the prevalence of hexosylated (usually glycosylated [52] ) metabolites. In addition, 301 the 92.19 fragment of unknown molecular composition is highly common and found in 302 over 25% of spectra. When comparing these overall fragment frequencies to all features 303 influenced by natural variation in our experiment ( _i.e._ associated with any mQTL), the 304 pattern was very similar (Supp. Figure 3), indicating the diversity of metabolites and 305 pathways influenced by natural variation. The presence of the most common fragments 306 in such a large proportion of features may cause them to overshadow the mQTL 307 specific fragments. To focus on mQTL-specific common fragments, these overall 308 common fragments are shown in gray in the following plots and are not counted as most 309 common or outstanding for mQTLs. The frequencies of the top 1,000 fragments are 310 listed in Supp. Table 5. 311 We treated the frequencies of fragments in all spectra as the expected values to test 312 whether mQTL feature groups had a significantly distinct commonality (Supp. Note 1). 313 Confident analysis required groups of at least 20 ddMS [2] spectra. Almost all high-impact 314 mQTLs, those above this size, had a significantly distinct fragment commonality, 315 compared to the α<0.05 significance threshold determined by random permutations 316 (Fig. 3CD), which would match scenarios I or II. Moreover, most high-impact mQTLs 317 (24/27 in positive mode and 19/23 in negative) also had outstanding fragments, whose 318 frequency exceeded the mean top frequency of size-matched random groups (Fig. 319 3EF), thereby supporting our hypothesis. 14 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 320 Figure 3: **Fragment commonality of mQTL feature groups** A and B. The most common fragments found in all ddMS² spectra obtained from tomato leaf metabolite features detected in positive (A) or in negative (B) ionization mode. The m/z value of each fragment is indicated above the corresponding bar, whose height indicates the frequency of ddMS² spectra which contain the fragment out of all detected ddMS² spectra. Due to the excess of fragments some numbers overlap, the ranked list of fragments and their frequency can be found in Supp. Table 5. Hexose (gray hexagon) derived fragments are marked with gray underline. C and D. The significance level of mQTL fragment commonality indicating how different the commonality of an mQTL feature group is compared to the commonality of all ddMS² spectra. Shown for features detected in positive (C) or in negative (D) ionization mode. The pink line denotes the significance threshold derived from the 5 [th] percentile of random ddMS² spectra sampling with corresponding number of ddMS² spectra. Dashed gray lines note the suggested minimum number of spectra for a confident analysis. E and F. The frequency of the most common fragment of each mQTL feature group (excluding the overall most common fragments from panels A-B), for features detected in positive (E) ionization mode or in negative (F). The light blue line denotes the outstanding threshold, the mean top frequency of random ddMS² spectra sampling with corresponding number of ddMS² spectra. Dashed gray lines note the suggested minimum number of spectra for a confident analysis. Colored arrows indicate the mQTLs studied in the following sections, in accordance with Fig. 2. 15 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 321 _NL_ 322 The fragments detected in ddMS [2] spectra are only the charged fragments of the 323 molecule. However, the uncharged fragments can be estimated by calculating the NLs 324 within each spectrum [34] . We calculated NLs for each spectrum as the m/z differences 325 between all fragments including the precursor m/z. The resulting set of NLs for each 326 spectrum is its NL pseudospectrum. These NLs can also be treated as descriptors of 327 features, and the most common NL could supplement the fragmentation commonality in 328 highlighting metabolites' uncharged backbone substructures. NLs have been 329 successfully used to annotate metabolite substructures and identify metabolic genes 330 associated with them [34] . Similarly to fragments, we identified the frequency of NLs in all 331 collected ddMS [2] spectra and tested the difference in frequency of NLs in mQTL feature 332 groups (Supp. Figure 4; Supp. Table 6). 333 Because NLs are calculated and not directly measured, NL pseudospectra contain 334 more noise derived from calculated differences which do not actually represent parts of 335 a molecule. For this reason, the statistical power of the commonality analysis is 336 reduced, and requires at least 50 ddMS [2] spectra (Supp. Note 1) for confident analysis. 337 It is tempting to consider only the annotated NLs (Supp. Table 7), which are likely true 338 molecular substructures. However, one strength of commonality analysis is that it 339 highlights biologically relevant NLs over random artifacts of the calculation. This can be 340 seen in high-impact mQTLs, which had a significantly distinct NL commonality and 341 outstanding NLs compared to random feature groups (8/16; Supp. Figure 4). For 342 example, an NL of unknown structure with m/z of 19.01 is highlighted as an outstanding 343 NL in nine different mQTLs (Supp. Table 4). 344 _Transformations_ 345 In scenario III, where the mQTL enzyme is involved in multiple pathways, fragments 346 or NL commonality may not be able to highlight a specific backbone. However, the 347 recurring m/z difference, or transformation, between pairs of features in an mQTL 348 feature group may highlight the mQTL-related reaction. Transformations, too, are a 349 calculated descriptor which can have noise due to the "blind" nature of the calculation. 350 Still, the most common transformations found between all pairs of features were known 16 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 351 metabolite modifications (Fig. 4AB; Supp. Table 8). For features detected in positive 352 mode, the most common transformations among all measured features are 353 de/hydration, de/hydrogenation and methoxylation (Fig. 4A). In negative mode the most 354 common transformations are methoxylation, de/hydration and acetylation (Fig. 4B). In 355 both modes, these are followed by the calculated transformation of 12.00. However, this 356 does not represent a real biochemical reaction, but rather the many pairs of molecular 357 formulas that differ by one carbon atom. This demonstrates one of the necessary 358 cautions when interpreting the calculated differences. Other common transformations 359 are additions/subtractions of methanol or ethyl groups and sugar (hexose, pentose) 360 modification. The statistical method development (Supp. Note 1) showed that the 361 transformation commonality of mQTL feature groups may confidently be interpreted as 362 significantly distinct from the expected, when applied to groups of at least 50 features 363 (Fig. 4CD). The most common transformation also exceeded the outstanding threshold 364 in most high-impact mQTLs (20/28; Fig. 4EF). 365 These results support our hypothesis that mQTL feature groups have distinct 366 commonly shared descriptors. In the following sections we demonstrate the utility of this 367 approach for highlighting previously uncharacterized metabolic gene functions. 368 17 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 369 Figure 4: **Transformation commonality for identifying mQTLs' underlying reaction** A and B. The most common transformations (m/z differences) found for pairs of all metabolite features detected in positive (A) or in negative (B) ionization mode. The m/z value of each transformation is indicated above the corresponding bar, whose height indicates the number of occurrences out of all non-redundant feature pairs. Due to the excess of transformations some numbers overlap, the full ranked list of frequent transformations can be found in Supp. Table 8. When available, putative annotations are noted in purple, according to Supp. Table 7. C and D. The significance level of mQTL transformation commonality indicating how different feature pairs of an mQTL feature group are compared to pairs of all features, for features detected in positive (C) or in negative (D) ionization mode. The pink line denotes the significance threshold derived from the 5 [th] percentile of random sampling with corresponding number of feature m/z values. Dashed gray lines note the minimum number of m/z values suggested for a confident analysis. E and F. The frequency of the most common transformation of each mQTL feature group (excluding the overall most common transformations from panels A-B), for features detected in positive (E) or in negative (F) ionization 18 mode. The blue line denotes the outstanding threshold, the mean top frequency of the most common transformation in random m/z samplings, with a corresponding number of feature m/z values. Dashed gray lines note the minimum number of m/z values suggested for a confident analysis. Colored arrows indicate the mQTLs studied in the following sections, in accordance with Fig. 2. [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 370 Highly common fragments in an mQTL feature group point to a sinapoyl-hexose 371 acyltransferase 372 An mQTL on chromosome 12 had the most significantly distinct fragment 373 commonality and the most outstanding fragments among features detected in negative 374 mode, and also a significantly distinct commonality and high fragment frequencies in 375 positive mode (blue arrow in Fig. 3C-F). The mQTL was associated with 134 features, 376 mostly in control conditions and very few in WD (Supp. Table 9). We present an in-depth 377 study of this mQTL as an example of how mQTL commonality can reveal novel 378 functions of genes in metabolic pathways. 379 The mQTL spanned 63,464,146-64,404,847 bp on chromosome 12 (Fig. 5A). The 380 majority of features were most significantly associated with single nucleotide 381 polymorphism (SNP) S12.2059255 at 63,718,397 bp. 382 To better understand the relations between the multiple mapped features and how 383 they may be affected by the mQTL alleles, all features were grouped by pairwise 384 correlation hierarchical cluster analysis (HCA) (Fig. 5B). The first and largest group 385 (designated group A) was mostly negatively correlated to group B. However, many 386 features of group B (group B1) which had mixed correlation to group A also had an 387 mQTL peak with SNP S12.2053749 (position 63,540,503 bp) and may not be directly 388 related to the mapping of other features. Group C had poor correlations with the other 389 groups, and inconsistent mapping patterns which may also be a result of closely located 390 associations mistakenly grouped in the same mQTL. The opposite relations of the HCA 391 groups to each other and to SNP S12.2059255 alleles (Fig. 5C) may suggest they are 392 grouped into the precursors and the downstream products of the mQTL-associated 393 reaction. 394 To focus annotation efforts on the substructures most central to the pathway affected 395 by the mQTL, we present commonality spectra of the mQTL feature group for features 396 detected in positive (Fig. 5D) or negative ionization mode (Fig. 5E). The outstanding 397 fragments in negative mode were 205.05 and 223.06, present in 68% of mapped 398 features, and multiple hexose fragments (Fig. 5E, in gray). Several of the most common 399 fragments in positive and negative could be paired by the typical mass difference of 2 19 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 400 between [M+H]+1 and [M-H]-1 ions. Consistently, fragment 207.07 was the most 401 common in positive mode, present in 38% of features. 402 403 Based on experimental [47,53] and predicted fragmentation [54], the fragment pairs 404 205.05/207.07 and 223.06/225.05 are typical of a sinapoyl substructure. The expected 405 225.05 had low abundance in our data. The fragments 367.10 and 369.11, also among 406 the most common fragments in negative and positive modes, respectively, match 407 sinapoyl-hexose (Fig. 5E, inset). The presence of hexose (162.05) in many of the 408 features in the mQTL feature group is also observable in the NL commonality in both 409 positive and negative modes (Supp. Figure 5AB), as is the common presence of the 410 sinapoyl (206.06, 224.07) and sinapoyl-hexose (368.11) in negative spectra NLs (Supp. 411 Figure 5B). Figure 5: **Study of chromosome 12 mQTL based on common fragments** A. A. Manhattan plot of the significant mQTL region in control conditions. Each color represents the mGWAS significance of one of the metabolite features associated with the mQTL. The lower panel shows the linkage disequilibrium analysis of all SNPs in the extended mQTL region. The red triangle indicates the linkage group of the mQTL significant SNPs. Orange rectangles denote the genomic location of serine carboxypeptidase-like (SCPL) family genes. B. Feature-feature correlations for all features of the mQTL feature group in control conditions. Features are ordered and squares are drawn according to the similarity calculated by hierarchical clustering analysis (HCA). Color scale represents the R [2] coefficient for feature pairs. HCA-based group names were defined according to correlations in control conditions and are noted on the side bar. C. The effect of the most significantly associated SNP, S12.2059255, alleles on relative feature abundance of a representative feature from each HCA group. D and E. Commonality spectra of the chromosome 12 mQTL feature group highlighting the ddMS [2] fragments detected in positive (D) or in negative (E) ionization mode. The fragments most common among all features are presented in gray. In E, fragments in gray are hexose-derived. Inset in E: molecular structure of sinapoyl-glucose. F. Correlations between expression of genes in the mQTL region (in rows) and features (in columns) significantly associated with the mQTL in control conditions. Color scale represents the R [2] coefficient, as in panel B. Yellow triangles point to the selected candidate gene, Solyc12g088250. G. Genome wide association (Manhattan plot) of Solyc12g088250 relative expression in control conditions. H. Phylogenetic clustering of genes from the SCPL family based on protein sequence alignment. The amino acids of the conserved family motif are marked by a red square. Genes expressed in the tomato leaf samples are highlighted in green. I. Pairwise gene expression correlation of the leaf-expressed genes of the SCPL family. Color scale represents the R [2] coefficient, as in panel B. J. Heatmap of the changes in features of the mQTL feature group in plants overexpressing Solyc12g088250 (golden top bar) compared to wild-type (green top bar). Features are ordered by HCA groups from panel B. (figure in next page 20 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 412 21 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 413 Of the features associated with the mQTL, only one was directly annotated: it 414 matched 1,2-di-O-sinapoyl-beta-D-glucose in the GNPS library (Supp. Table 9), and four 415 others were linked to sinapoyl-hexose in the FBMN (Supp. Table 10). However, we were 416 able to identify the typical fragment in 58 of the features in the group (of the 97 with 417 available ddMS2 spectra) and could therefore extend the substructure annotation to 418 43% of the mQTL feature group. 419 The mQTL region included 28 genes, but none was annotated to catalyze an 420 expected metabolic reaction. Correlations between features of the mQTL feature group 421 and relative expression of genes within the mQTL region pointed to a gene, 422 Solyc12g088250, with many significant correlations (p-value < 10-7, |r|>0.4; Fig. 5F). 423 Solyc12g088250 is located at 63,710,096 bp, right at the peak of the mQTL and is 424 annotated (ITAG2.4) as a serine carboxypeptidase, one of four tandem genes with the 425 same annotation in the mQTL region. Analysis of gene expression levels by GWAS 426 revealed a cis expression QTL (eQTL) co-locating with the mQTL (Fig. 5G), suggesting 427 that a polymorphism affecting Solyc12g088250 expression underlies the mQTL. 428 We searched in Plant Metabolics Network (http://www.plantcyc.org) databases and 429 the literature [55–57] for enzymes involved in sinapoyl-hexose metabolism and found 430 multiple mentions of enzymes of the serine carboxypeptidase-like (SCPL) family. We 431 therefore performed a detailed homolog phylogenetic analysis (Supp. Figure 5C) which 432 revealed that Solyc12g088250 is in fact part of the SCPL subfamily, a group of 433 metabolic genes with a conserved motif of GDSYS (S is the essential catalytic serine) [56], 434 with no peptidase activity. The redefined gene family annotation and common 435 substructure annotation supported the involvement of Solyc12g088250 in sinapoyl 436 hexose transferase activity. We therefore set out to validate the gene's activity by 437 overexpressing (OE) Solyc12g088250 in tomato plants. 438 439 _Validation of gene-metabolite association by Solyc12g088250 overexpression_ 440 The Solyc12g088250 OE plants of five independent transformation lines (Supp. 441 Figure 6A) were visually comparable to the wild-type (WT) plants (accession Ailsa 442 Craig), but metabolic analysis revealed clear differences (Supp. Figure 5D). Of the 134 443 features associated with the mQTL, 59 features were also found in the OE experiment. 22 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 444 Most of them (43/59) were significantly altered in OE samples compared to WT (Fig. 445 5J). 446 Features of the large HCA group A were significantly increased by OE, in agreement 447 with the positive correlation between gene expression and their feature abundance seen 448 in the natural population. The features in group B2 were strongly associated with the 449 mQTL lead SNP, and the four which were found in the OE experiment displayed 450 decreased abundance in OE samples, as expected from their negative correlation to 451 group A and to Solyc12g088250 gene expression. 452 Additional features which were not associated with the mQTL were affected by the 453 OE. Partial least squares discriminant analysis (PLS-DA) was used to identify which 454 features contributed most to the difference between OE and WT samples. We selected 455 the top 1% of contributing features (168 in positive and 85 in negative ionization modes) 456 (Supp. Figure 5E). These were also significantly different (t-test) between OE and WT. 457 Commonality analysis of these features revealed hexose fragments and sinapoyl 458 hexose fragments (205.05, 223.05 and 367.1) among the most common in negative 459 mode and sinapoyl (207.07) in positive mode (Supp. Figure 5F). However, in addition to 460 the expected fragments, reproducing the central mQTL substructures, additional 461 fragments were common among the features most distinguished between OE and WT. 462 Possibly, strong OE reveals a more extended effect of the gene on metabolites whose 463 abundance changes in the natural variation are mitigated by interconnected metabolic 464 networks. Alternatively, OE of a gene may at times produce unnatural promiscuous 465 activity, which extends beyond the naturally occurring reactions [58,59] . Overall, the OE 466 experiment confirmed that the SCPL gene Solyc12g088250 affects the features of the 467 mQTL feature group within the sinapoyl-hexose pathway. 468 469 Shared and divergent activity of two SCPL genes 470 Twelve genes in the tomato genome fell within the SCPL subfamily clade of the 471 phylogenetic tree (Fig. 5H). Four of them were expressed in leaf samples in this 472 experiment (Fig. 5H, green highlight). Solyc12g019570 encodes a truncated protein 473 without the GDSYS motif and is probably a pseudogene. Solyc10g049210 is an 23 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 474 ortholog of Sopen10g020280 from tomato relative _S. pennellii_, previously characterized 475 as an isobutyrylglucose:isobutyrylglucose acyltransferase (GAT) [60] . 476 Interestingly, the expression of Solyc04g077650 was positively correlated with that 477 of Solyc12g088250 (Fig. 5I). The sequence similarity and coordinated expression raised 478 the question of whether these genes were redundant or served distinct metabolic 479 reactions. We hypothesized that commonality analysis could be utilized to highlight 480 important descriptors for features grouped by methods other than co-association. Here, 481 we used it in a reverse-genetics approach to identify overlap and specificity between 482 related SCPL enzymes. 483 To this end, we produced plants overexpressing Solyc04g077650 (Supp. Figure 6B). 484 Again, PLS-DA was used to identify which features contributed most to the difference 485 between OE and WT samples (Fig. 6A). The fragment commonality of this OE-affected 486 feature group was significantly (p-value < 10 [-16] ) distinct from the expected. The most 487 common fragments were similar to the ones found in commonality of the 488 Solyc12g088250 OE, albeit different in order of frequency. To further refine the way 489 commonality analysis helps us understand feature-gene relations, we divided the group 490 of OE-affected features into subgroups: features increased by OE, expected to be 491 downstream products of the enzyme reaction, and those decreased by OE, expected to 492 be upstream precursors. Side-by-side comparison of the commonality of predicted 493 precursors and products of the two genes is presented in Fig. 6. 494 For both genes many precursors and products contain hexose fragments. The 495 typical sinapoyl-hexose fragments are found in the products of both genes but in the 496 precursors of Solyc12g088250 only, suggesting that this enzyme may add sinapoyl and 497 hexose to molecules already containing these groups (Supp. Figure 5G). In addition, 498 apart from hexose fragments, the fragments found in more than 10% of precursors are 499 different between the genes, suggesting they perform a similar function of sinapoyl 500 acyltransferase with different substrate specificity. 501 502 Overall, the most common fragments point to sinapoyl-hexose as a core 503 substructure in the pathway altered by the mQTL. Correlation and co-locating eQTL of 504 Solyc12g088250 gene expression suggest that genomic variation alters expression, in 24 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 505 turn influencing enzyme activity levels and metabolite abundance in the sinapoyl 506 hexose pathway. Solyc12g088250 was misannotated by automated annotation. Without 507 the rapid identification of the sinapoyl-hexose substructure by commonality analysis, 508 this connection would have been harder to make. We further demonstrate the utility of 509 commonality analysis to quickly highlight metabolic pathways affected by experimental 510 genetic modifications, which may be used in forward and reverse-genetics approaches 511 to identify metabolic gene functions. Figure 6: **Fragment commonality of predicted precursors and products of two serine** **carboxypeptidase-like genes** A. A. Volcano plot displaying the significance of a t-test comparing metabolite feature abundance between plants overexpressing (OE) Solyc04g077650 and wild-type (WT; Ailsa Craig accession) against the fold change (after log2 transformation) of the same comparison. Features which contributed most (top 1%) to the separation of OE and WT in PLS-DA are colored golden (increased in OE) or green (decreased in OE). B-E. Commonality spectra of fragments of features decreased by OE of Solyc04g077650 (B, D) or Solyc12g088250 (C, E) in positive (B, C) or negative (D, E) ionization modes. F-H. Commonality spectra of fragments of features increased by OE, expected to be the products of the reaction as denoted by the arrows, of Solyc04g077650 (G) or Solyc12g088250 (F, H) in positive (F) or negative (G, H) ionization modes. 512 513 514 25 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 515 WD-specific regulation of a putative transporter influences accumulation of multiple 516 nitrogen-containing metabolites 517 The next mQTL had one of the most significantly distinct fragment commonalities in 518 positive ionization mode, but its interpretation was more challenging since the main 519 fragment structure is unknown. The mQTL spanned 40,443,154-40,914,962 bp on 520 chromosome 2. It was significantly associated with 49 features (Fig. 7AB). In WD, 35 521 features were associated, five of which were also significantly associated in control, but 522 with lower significance (Supp. Table 11; Supp. Figure 7). In contrast, 14 features were 523 significantly associated only in control (Fig. 7A) but showed the highest association with 524 a SNP 60 kb away from the WD-associated features. 525 The 49 features of the mQTL feature group were categorized into three groups by 526 correlation-based HCA (Fig. 7C). Groups A and B were negatively correlated to each 527 other and showed opposite effects of the most significantly associated SNP 528 S02.1211713 (Fig. 7D). 529 Of the 95 genes in the mQTL region, 34 were expressed in the leaf samples. One 530 gene, Solyc02g071100, was highly correlated (p-value < 10 [-7], |r|>0.4) with most features 531 that were associated with the mQTL in WD conditions (Fig. 7E). The expression of 532 Solyc02g071100 across all samples was elevated in WD samples compared to the 533 control (Fig. 7F). The alleles of the lead SNP of the mQTL had a significant effect on the 534 gene expression in WD, but not in control samples (Fig. 7F). This association was 535 further supported by an eQTL for Solyc02g071100 expression co-locating with the 536 mQTL exclusively in WD conditions (Fig. 7G). 537 Solyc02g071100 is annotated as a purine permease gene with a GO term of purine 538 nucleobase transmembrane transporter activity. It was found to be highly similar to nine 539 other genes in the tomato genome annotated as purine permeases, three of them within 540 the same mQTL region (Fig. 7B). Neither adenine nor guanine (Supp. Table 3) was 541 associated with the mQTL. In fact, all associated features had much higher molecular 542 weight, pointing to more complex compounds. 543 The majority of associated features were detected in positive ionization mode, so 544 commonality analysis was only possible for these features. The most common 545 fragments in the mQTL feature group were 222.11, 72.08, 163.04, 223.12 and 129.14 26 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 546 (Fig. 7H), but the substructure they are derived from is not known. Among the common 547 fragments of the mQTL feature group many had even nominal m/z values, suggesting 548 nitrogen-containing compounds (nitrogen rule for protonated ions) [61] . The only annotated 549 feature in the group was feruloylputrescine, and two others had some similarity to this 550 compound. Together with the putative annotation of Solyc02g071100 as a purine 551 permease and the association with feruloylputrescine, this led us to hypothesize that the 552 mQTL is involved in the metabolism of nitrogen-containing compounds. 553 _Annotation support by stable isotope labeling_ 554 To test this hypothesis and improve annotation, we conducted an isotope labeling 555 experiment using [15] N-labeled plants and [13] C-labeled calli. The feruloylputrescine 556 annotation, with the C 14 H 20 N 2 O 3 predicted molecular formula, was confirmed by a shift 557 of the precursor ion m/z of 265.1547 by 1.9941 in [15] N-labeled samples and by 14.0469 558 in [13] C labeled samples (Fig. 7I). We continued investigating all features of the mQTL 559 feature group which were found in the labeled samples in a similar fashion (Supp. 560 Figure 8). Because specialized metabolite production depends on environmental 561 conditions, it is expected to differ in tissue culture. Of the 41 features of the group in 562 positive mode, 20 were found in the [15] N-labeled samples, and fewer in the [13] C-labeled 563 calli. All detected features were demonstrated to contain between one and four nitrogen 564 atoms in their molecular formula (Supp. Figure 8). 565 Of the features in the group, 22 contained the 222.11 fragment. Of these, nine were 566 found in the [15] N-labeled samples and were confirmed to contain 2-4 nitrogen atoms. We 567 conclude that the mQTL influences nitrogen metabolism, and next validated the 568 involvement of the candidate gene. 569 _Validation of nitrogen metabolism function by Solyc02g071100 overexpression_ 570 We next set out to validate the role of Solyc02g071100 as the gene underlying the 571 mQTL by OE (Supp. Figure 6D). Since the effect of the mQTL was more pronounced in 572 WD conditions, the OE lines and control WT were subjected to reduced water 573 availability prior to metabolic profiling. As expected, the developmental stage, sampling 574 time and water availability influenced the general metabolic profiles (Supp. Figure 7C). 575 Remarkably, the overexpression of Solyc02g071100 also had a substantial effect on the 27 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 576 general leaf metabolic profile (Supp. Figure 7C). Of the 41 features of the mQTL feature 577 group detected in positive mode, 25 were found with high confidence in the OE 578 experiment samples. All but three displayed a strong reduction in OE samples (Fig. 7J), 579 as expected from the negative correlation between gene expression and feature 580 abundance in the GWAS population. The effect of OE was high regardless of the 581 watering regime, probably due to constitutive high expression by the 35S promoter 582 overriding any condition-dependent gene regulation. 583 Following the strong effect of Solyc02g071100 OE on the metabolic profile, we 584 identified the top 1% of features most influencing the separation between WT and OE 585 using PLS-DA (Supp. Figure 7F). Again, the most common nominal m/z values in this 586 group of features were even, suggesting nitrogen-containing compounds (Supp. Figure 587 7G). Fragment 222.11 was highly common in OE-decreased features, as in features 588 negatively correlated with gene expression in the natural population, but absent from 589 features increased by OE. 590 Interestingly, many features which increased in OE had very high m/z values (800 591 1486), while features with lower abundance in OE had mostly m/z values under 600, 592 indicating that the consequence of this gene's activity may involve synthesis of more 593 complex compounds rather than degradation. Also, the high m/z values of the most 594 common transformations between features increased versus decreased by OE point to 595 the addition of high-molecular-weight groups to the compounds (Supp. Figure 7H). In 596 support, the only three features in HCA group A (positively correlated with gene 597 expression) have the three highest m/z values among the features of the mQTL feature 598 group (Supp. Table 11). 599 Taken together, we propose Solyc02g071100 as a WD-induced transporter of 600 nitrogen-containing compounds and hypothesize that it transports the compounds into a 601 compartment where they are metabolized to form high-molecular-weight compounds. 602 603 604 28 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 605 Figure 7: **Study of chromosome 2 mQTL and candidate gene** A-B. Manhattan plot of the significant mQTL region in control (A) and water deficit (B) conditions. Each color represents the GWAS results of one of the metabolite features associated with the mQTL. The lower panel in B shows the linkage disequilibrium of all SNPs in the extended mQTL region. The red triangle indicates the linkage group which includes the mQTL significant SNPs. Orange rectangles denote the genomic location of purine permease family genes. C. Feature-feature correlations for all features of the mQTL feature group in water deficit conditions. Features are ordered and grouping squares are drawn according to the similarity calculated by hierarchical clustering analysis (HCA). Color scale represents the R [2] coefficient for feature pairs. HCA-based group names are noted on the side bar. D. The effect of the most significantly associated SNP, S02.1211713, alleles on relative feature abundance of representative features from the HCA groups in water deficit. E. Correlations between expression of genes in mQTL region (in rows) and features (in columns) with significant association to the mQTL in water deficit conditions. Color scale represents the R [2] coefficient, when significant. The yellow triangle points to the selected candidate gene, Solyc02g071100. F. The effect of the lead SNP S02.1211713 and of the environment on Solyc02g071100 gene expression. G. Genome wide association (Manhattan plot) of Solyc02g071100 relative expression in water deficit conditions. H. Fragment commonality spectrum of the mQTL feature group, for features detected in positive ionization mode. I. Spectra at retention time 4.44 min derived from [15] N and [13] C isotope labeling experiments. The comparison of spectra of this representative feature shows the increase in m/z corresponding to the number of labeled atoms, validating the predicted (C 14 H 20 N 2 O 3 ) molecular formula of feruloylputrescine. J. Heatmap of the changes in features of the mQTL feature group in plants overexpressing Solyc02g071100 (golden top bars) compared to wild-type (green top bars), under fully watered or water deficit conditions at the indicated time points. Feature side colors indicate the HCA groups from panel C. (figure on next page) 606 29 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 607 608 30 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 609 Outstanding transformation in an mQTL feature group points to a double bond 610 reductase 611 We next tested the ability of transformation commonality to reveal the reaction 612 underlying an mQTL. One of the mQTLs with the highest frequency of an outstanding 613 transformation in positive and negative modes was found on chromosome 9. It was 614 consistently associated with a high number of features, 271 and 126 in control and WD 615 conditions, respectively (Supp. Table 12). 616 The mQTL spanned 70,800,000-71,600,000 bp (Fig. 8AB). Detailed examination 617 revealed a strong association with a single SNP in the mQTL located at 70,943,627 bp, 618 with very low linkage (LD<0.1) to nearby SNPs. The SNP was positioned on the last 619 base of the predicted stop codon of Solyc09g091700 (Fig. 8D), eliminating the predicted 620 stop codon (ITAG2.4) and adding four amino acids, Trp-Lys-Lys-Lys, to the end of the 621 protein. Expression levels differed substantially in the population. However, expression 622 was not associated with the lead SNP alleles or correlated to the abundance of features 623 in the mQTL feature group, suggesting the coding polymorphism, not expression 624 differences, causes this mQTL. The reference allele A (genome version SL2.50) was 625 ubiquitous in _S. pimpinellifolium_ and _S. lycopersicum_ var. cerasiforme, but allele G 626 appeared among heirloom, transitional and modern accessions at a rate of 0.36, 0.18 627 and 0.27, respectively. The position of the stop codon in the A allele also matches that of 628 other homologous proteins (Supp. Figure 9A) and is likely the more active version. 629 The correlation-based HCA of features associated with this mQTL showed two major 630 groups, A and B, negatively correlated to each other, and a third group C, with mixed 631 relations to the first two groups (Fig. 8C). Groups A and B had consistently opposite 632 association with the lead SNP alleles in both conditions (Fig. 8E). Because group A 633 features increase with the likely active allele A, they are probably the products of the 634 reaction, and group B features are likely precursors. 635 The transformation commonality revealed 2.02 as the most frequent transformation 636 among all feature pairs of the mQTL feature group (Supp. Figure 9BC). Although 2.02 is 637 one of the common transformations found among all features (Fig. 4AB), it stands out 638 with an exceptionally high frequency among features associated with this mQTL. To 31 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 639 further refine the transformations most relevant for the enzymatic reaction driven by the 640 mQTL, we calculated transformations between the likely precursors and products. We 641 calculated the commonality of m/z differences between HCA group A and group B 642 features. In this directional analysis, 2.02 stood out at an even higher frequency in both 643 positive and negative modes (Fig. 8FG), suggesting that the mQTL-altered reaction 644 adds 2.02 Da to the precursor compounds. The 2.02 Da m/z difference is consistent 645 with that of 2 hydrogen atoms, or H 2 . Once we identified 2.02 as the transformation 646 central to the mQTL, we also noticed the same mass difference of 2.02 between 647 common fragments of the mQTL feature group (Supp. Figure 9DE). 648 Based on spectral matches to databases we were able to annotate seven features of 649 the mQTL feature group (Supp. Table 12). An additional seven features had partial or 650 substructure matches, three of which included a carveol substructure. The carveol main 651 library fragments were among the fifty most common fragments of the mQTL feature 652 group: 135.11, 119.085 and 137.096 (Supp. Figure 9D), and at least one was present in 653 46 (56%) of the features that had a ddMS [2] spectrum in positive ionization mode. All 654 structures we annotated for features of the mQTL feature group include a cyclohexene 655 or cyclohexane ring (Supp. Figure 9J). 656 The gene affected by the lead SNP, Solyc09g091700, is annotated as an alcohol 657 dehydrogenase, zinc-containing gene, with two specific GO terms: 2-alkenal reductase 658 [NAD(P)] activity and alcohol dehydrogenase (ADH) [NAD] activity. Solyc09g091700 659 was putatively assigned to ten alcohol dehydrogenase (EC 1.1.1.1) reactions in the 660 SolCyc database [62], producing aldehydes. However, none of the annotated features in 661 the group included an aldehyde. The closest Arabidopsis homolog, AT5G16970 (top 662 BLASTP hit in TAIR, 49% identity, 69% positives), is annotated as a 2-alkenal reductase 663 (EC 1.3.1.102). A similar gene, _Nicotiana tabacum_ double bond reductase ( _NtDBR_ ), has 664 been shown to catalyze reduction of double bonds on a variety of substrates [63], with 665 similarity to the structures we found (Supp. Figure 9J). These reactions incorporate H 2 666 into the molecule, consistent with the common 2.02 transformation of the mQTL. 667 _Validation of Solyc09g091700 activity by overexpression_ 668 Next, we set out to validate this hypothesis by OE of Solyc09g091700, using the A 669 allele coding sequence (Supp. Figure 6C). We found very high similarity between 32 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 670 commonality spectra of the mQTL feature group and of the top 1% of features 671 contributing most to the separation between WT and OE samples (Supp. Figure 9D-I). 672 This allowed a quick confirmation that Solyc09g091700 affects the features associated 673 with the mQTL. Of the 355 features which were associated with the mQTL, 83 were 674 found again in the OE experiment. The majority of them (70/83) were significantly 675 different between OE plants and WT (accession Micro-Tom) (Fig. 8H). Features in HCA 676 group A significantly increased in OE, while those in group B decreased. This agrees 677 with A being the functional allele, with group A features as products and group B 678 features as precursors. Directional transformation analysis between features increased 679 by OE and those decreased highlighted 2.02 as the most frequent transformation, as in 680 the directional transformation analysis of the mQTL. 681 Taken together, we propose that Solyc09g091700 functions as a double bond 682 reductase, highlighted by the 2.02 transformation. It likely accepts multiple substrates, 683 which often include cyclohexene or cyclohexane substructures. The coding variation in 684 this gene alters enzyme activity and influences the abundance of hundreds of features 685 in both WD and control conditions. 686 Figure 8: **Study of chromosome 9 mQTL based on high transformation commonality** A and B. Manhattan plot of the significant mQTL region in control (A) and water deficit (B) conditions. Each color represents the mGWAS significance of one of the metabolite features associated with the mQTL. C. Feature-feature correlations for all features of the mQTL feature group in water deficit conditions. Features are ordered and squares are drawn according to the similarity calculated by hierarchical clustering analysis (HCA). Color scale represents the R [2] coefficient for feature pairs. HCA-based group names were defined according to correlations in control conditions and are noted on the side bar. D. Focus on the genomic location of the single highly associated SNP, within the context of the gene Solyc09g091700. E. The effect of the highly associated SNP, S09.2367600, alleles on relative feature abundance of representative features from each HCA group in control (light blue) and water deficit (orange). F and G. Commonality of transformations between features in HCA groups A and B, which are negatively correlated to each other (panel C), for features detected in positive (F) or in negative (G) ionization mode. Inset: molecular structures of annotated metabolites with the typical double-bond of the mQTL feature group in positive (F) or in negative (G) ionization mode. H. Heatmap of the changes in features of the mQTL feature group in plants overexpressing Solyc09g091700 (golden top bar) compared to wild-type (green top bar). Features are ordered by HCA groups from panel C. (figure next page) 33 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 687 688 34 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 689 **Discussion** 690 Application and limitations of commonality analysis 691 In this study, we set out to reveal novel metabolic regulation in tomato leaves using 692 natural variation and environmental stress. The combination of genomic, transcriptomic 693 and metabolomic variation in natural populations has been successfully used in recent 694 studies to identify metabolic gene functions [8] . However, studies have often been limited 695 to annotated or partially annotated compounds. Some previous studies have identified 696 and quantified thousands of features but proceeded to map and study only hundreds of 697 them [36,37,64] . Currently, new discoveries are less limited by data acquisition and more by 698 processing and prioritizing deeper investigation. In this era of big data, the commonality 699 analysis presented here can facilitate gene-metabolite association and functional 700 annotation of previously understudied metabolites. 701 It has been shown that substructure annotation or semi-annotation is useful for 702 discovery of novel gene-metabolite associations [34] . Zhu _et al._ studied genomic hotspots 703 where multiple metabolites sharing the same identified substructure, a fragment or NL, 704 co-associate. The FBMN approach groups features by fragment similarity to identify the 705 underlying structures and group related metabolites. However, when grouping based 706 solely on similarity, not all pathway metabolites which share the same fragments will be 707 tied together in the FBMN, as seen by the dispersed placing of sinapoyl containing 708 features linked to the Solyc12g088250 mQTL. The commonality analysis combines the 709 advantages of both approaches, grouping by genomic co-association and using the 710 outstanding fragments to define the pathway backbone structure. Grouping by genomic 711 or other biological association may provide a more direct route to biological 712 interpretation. 713 We show that most high-impact mQTLs have a significantly distinct fragment 714 commonality and outstanding fragments suggested to belong to the backbone of 715 affected pathways. Still, commonality analysis is a prioritization and interpretation tool, 716 not a replacement for metabolite identification or functional validation. The simplicity of 717 the proposed calculations, however, allows for quick evaluation of fragment, NL and 718 transformation commonality for each mQTL or other studied group. Nevertheless, the 35 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 719 method is not universally effective. Depending on the genetic polymorphism and the 720 structure of the metabolic pathway, some mQTLs may have no outstanding descriptor. 721 The major limitation of the method is that it requires a large mQTL feature group with 722 ddMS [2] spectra for its features. Another is that the statistical framework does not 723 currently correct for the testing of multiple mQTLs (or other groups) in the experiment. 724 The use of NL commonality comes as a complementary method to using fragments 725 for semi-annotation. However, the presence of noise from calculated differences which 726 are not substructures can hamper the analysis. It therefore requires an even higher 727 number of ddMS [2] spectra for a meaningful analysis. If available, the commonality 728 method has the advantage of identifying highly recurring, and therefore possibly 729 important, NLs, even if they are not annotated. For example, nine mQTLs have an 730 outstanding NL of 19.01 in 37-61% of the features in the group (Supp. Table 4), with an 731 unknown formula, but perhaps worth investigating. 732 The transformation commonality analysis builds upon the CSPP approach [29,31], 733 where the m/z difference between two features points to the biochemical reaction 734 turning one to the other. However, instead of suggesting a reaction for each pair of 735 features, it draws upon the dominant transformation for the whole group, or between 736 groups, to suggest the reaction influencing them all. Transformation analysis benefits 737 from using precursor ion m/z values and therefore does not require ddMS [2] spectra. 738 However, it shares the disadvantage of using calculated descriptors which do not all 739 represent true metabolic differences. 740 The poorly correlated subgroups of features observed at several loci illustrate both 741 the risk of automated grouping by co-association and the ability of dominant common 742 descriptors to remain detectable despite some unrelated features. 743 The highlighted descriptors can then be evaluated in combination with other 744 available information, such as gene expression correlations, gene annotations and 745 similarity to annotated compounds to identify targets for investigation. 746 747 Common sinapoyl fragments reveal an SCPL-associated metabolic pathway 748 A clear example can be found with the mQTL on chromosome 12. The two 749 outstanding fragments quickly pointed to the sinapoyl-hexose pathway. The gene 36 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 750 expression correlation implicated Solyc12g088250, which we identified as an SCPL 751 family gene based on the common fragments. Sinapoyl-glucose and other cinnamic 752 acid derivatives were found in many studies to be substrates or products of SCPL 753 enzymes [50,60,65–79] . Classification of Solyc12g088250 as an SCPL gene was then 754 confirmed by the conserved motif. 755 Despite the original discovery in _S. pennellii_, the vast potential in metabolite 756 production, involvement of SCPL genes and their products in biotic stress response [80,81] 757 and some bioinformatic studies [82,83], the SCPL genes have not been experimentally 758 characterized in tomato or the closely related _Solanaceae_ . Characterized reactions of 759 SCPL genes mostly involve small metabolites. However, in this study we found several 760 large soluble metabolites related to sinapoyl derivatives that, to the best of our 761 knowledge, were previously unknown. 762 The SCPL gene family is fast expanding [83,84] and is often found in tandem in genomic 763 clusters. Tissue-specific expression, as seen previously and in our results [57], contributes 764 to diversification. Furthermore, despite close sequence similarity, genes were found to 765 have diverse substrates [57,72,74], with only a few amino acid changes sufficient to alter 766 substrate binding abilities [50] . We suggest fragment commonality as a tool to highlight 767 enzyme redundancy and specificity in substrates and products. As we show, features 768 increased by OE of Solyc12g088250 and of Solyc04g077650 share some of the most 769 common fragments, but other fragments are different. Together with tissue-specific 770 expression this can be used to study redundancy and divergence of metabolic gene 771 families. Notably, comparing commonality is possible also when direct feature 772 comparison is not straightforward, for example between experiments run with different 773 UPLC methods. 774 775 Water deficit exposes a purine-permease-associated network of nitrogen-containing 776 metabolites 777 On chromosome 2, we studied an mQTL with a clear WD-dependent regulation. 778 Both feature and gene expression (eQTL) associations co-located to the same mQTL 779 and were more significant in WD conditions (Fig. 7), indicating that the same 780 polymorphism likely regulates both gene expression and feature levels. The remarkably 37 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 781 consistent correlations between feature abundance and Solyc02g071100 gene 782 expression were also seen predominantly in WD, exemplifying the importance of 783 experiments exposing plants to diverse conditions to reveal more of the world of 784 specialized metabolites and their regulation. 785 Solyc02g071100 is annotated as a purine permease, a member of the purine uptake 786 permeases (PUP) gene family. PUP proteins are membrane-localized proton 787 symporters with multiple transmembrane domains [85] . First identified as adenine 788 transporters, PUPs also transport pyridoxine (vitamin B6) [86], kinetin, alkaloids including 789 caffeine [87,88] and nicotine (thereby influencing alkaloid metabolism [89,90] ), and 790 benzylisoquinoline alkaloid precursors and opiates [91] . Another important substrate of this 791 gene family is cytokinin [92], as found in rice [93] and Arabidopsis [94] . Recently reviewed by 792 Demurtas et al. [95], the PUP family has been shown to transport different purine ring 793 containing compounds, with distinct substrate preferences for each gene. All known 794 substrates contain nitrogen, which was also highlighted in our commonality analysis and 795 validated by the [15] N labeling experiment. 796 One feature associated with the mQTL was annotated as feruloylputrescine, and two 797 others were similar to it. First found and named in citrus [96], feruloylputrescine is a 798 hydroxycinnamic acid amide. It shows an organ-specific distribution in several plants [97] 799 including tomato [98,99], where it is abundant in developing fruit [100–102] . It is present in 800 leaves and accumulates upon exposure to _Tuta absoluta_ [99] and _Pseudomonas_ 801 _syringae_ [98] . A direct link between feruloylputrescine and abiotic stress has rarely been 802 reported in the literature [103,104], and we found no previous reports of water-stress-related 803 regulation. 804 Features of the larger subgroup of the mQTL feature group (HCA group B), including 805 feruloylputrescine, were negatively correlated to gene expression and decreased in 806 leaves with Solyc02g071100 OE, and are therefore expected to be the precursors in the 807 metabolism associated with PUP activity. They are characterized by common fragment 808 222.11, which is not found in the products increased by Solyc02g071100 OE. These 809 products are mostly high-molecular-weight compounds, larger than any compounds 810 described in the literature to be transported by PUP, suggesting that they are formed in 811 the subcellular compartment into which Solyc02g071100 transports the precursors. 38 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 812 Compartmentalization is known to be very important for cellular regulation of 813 metabolite production and accumulation [105,106] . PUP family genes have gene-specific 814 subcellular localization [107], in either the plasma membrane, tonoplast or endoplasmic 815 reticulum (ER). Like other purine permeases, Solyc02g071100 is predicted to be a 816 transporter with multiple transmembrane domains (Supp. Figure 7D), but its membrane 817 localization is not known. Subcellular localization prediction suggests a membrane 818 localization, either plasma membrane, ER, tonoplast or Golgi (Supp. Figure 7E). 819 To conclude, although the full mechanism of Solyc02g071100 transport and 820 metabolic influence could not be elucidated by our analysis, commonality analysis and 821 15 N labeling point to a role for this gene in the transport of multiple nitrogen-containing 822 compounds larger than previously described. 823 824 Transformation commonality links _Solyc09g091700_ to double-bond reduction 825 On chromosome 9, we found a strong mQTL that was consistent across conditions 826 (Fig. 8). The mQTL was associated with a single SNP with very low LD to others, which 827 is unusual for this population with dense SNPs and LD decay of up to 800 kb [2] . This is 828 probably because the SNP arose relatively recently: it first appears in heirloom varieties, 829 at most hundreds of years ago [2] . 830 Solyc09g091700 is homologous to alkenal double bond reductase (DBR), or 2 831 alkenal reductase (AER) [108,109] family genes. These NADPH-dependent enzymes have 832 been studied in various plant species and share two general motifs AXXGXXG and 833 GXXGXXG [109], found also in Solyc09g091700. Enzymes of this family catalyze the 834 reduction of the α,β-unsaturated bond of reactive carbonyls [108,110], nitro-fatty acids (NO2 835 FAs) [111] and the reduction of cinnamaldehydes to 836 hydroxydihydrocinnamaldehydes [63,109,112] . The outstanding 2.02 transformation between 837 group B and group A features points to a reduction reaction, matching the proposed 838 function of the gene as a DBR that reduces a C=C double bond by adding H 2 . 839 Consistently, the common fragments of HCA group A features were 2.02 m/z units 840 higher than the common fragments found for B group features (Supp. Figure 9DE). 841 A tomato AER, encoded by Solyc05g005480, reduces the double bond adjacent to 842 the carbonyl group of Furaneol-like flavor compounds [113], which have a five-atom ring. In 39 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 843 contrast, Solyc09g091700 appears to have a different substrate specificity: its 844 associated features include a cyclohexene or cyclohexane ring. Within this gene family, 845 the enzymes whose substrates are most similar to these metabolites, including a 846 cyclohexene ring, were found in _Nicotiana tabacum_ [63] and peppermint [114] . The _Nt_ DBR 847 converts perillaldehyde to dihydroperillaldehyde, similar to the perillyl alcohol and 848 carveol substructures found in the mQTL-associated compounds [63] . Carveol, which is a 849 substructure found in some associated metabolites, is a common volatile [115] . It has an 850 inhibitory effect on food-borne harmful bacteria [116], a fungus [117], and harmful insects [118,119] . 851 β-ionone, one of the annotated associated features, is a volatile involved in plant-insect 852 interactions [120] and gives tomato fruit a pleasant scent [121,122] . Although scarcely 853 mentioned in the literature, the Solyc09g091700 gene expression was up-regulated in a 854 tolerant genotype upon exposure to _Tuta absoluta_ [123], and is co-expressed with 855 insecticidal acyl-sugar acyltransferases [124] . Understanding the role of these compounds 856 and genes in biotic stress responses may help improve plant interactions with insects 857 and pathogens. The proposed function of Solyc09g091700, with its broad range of 858 putative substrates, is a step in this direction. 859 Broader uses and outlook 860 Although commonality analysis was developed here around mQTL feature groups, 861 its underlying principle is broader: features can be grouped by any biologically 862 meaningful criterion and then examined for enriched structural or transformation 863 descriptors. Examples include features that respond significantly to environmental 864 perturbation, or features that change significantly in plants with a specific gene knockout 865 or OE, as shown for the four genes in this work. Features significantly and tightly 866 correlated under any set of variable conditions or present in a specific set of samples, 867 such as an exotic phylogenetic group, can also be subjected to commonality analysis to 868 highlight shared structures. With such grouping, stricter significance thresholds give 869 cleaner groups. With a less stringent group, sub-grouping by fragment similarity by 870 FBMN or another biologically relevant method may be useful. Another benefit of the 871 approach is that it can be applied in retrospect to existing data from existing projects, 872 provided enough MS [2] spectra are available. 40 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 873 Many mQTLs not examined here contain candidate genes and recurrent descriptors 874 that can be revisited as annotation resources and analytical methods improve. The 875 fragmentation and gene expression information still holds potential for further analysis. 876 The systematic analysis of mQTL commonality showed that a significantly distinct 877 commonality and outstanding descriptors, either fragments, NLs or transformations, are 878 found in most high-impact mQTLs. The three loci studied here demonstrate how 879 commonality analysis can help glean investigation-worthy descriptors and associations 880 from untargeted metabolite profiling. Testing different descriptor types provides 881 complementary applications: recurrent fragments exposed a shared pathway backbone, 882 recurrent mass differences suggested an enzymatic reaction, and condition-specific 883 fragment patterns guided isotope-based characterization of transporter-associated 884 metabolites. Used together with genetic, transcriptomic, and experimental evidence, 885 commonality analysis provides a flexible route from large, incompletely annotated 886 metabolomic datasets to biologically relevant features, independent of prior annotation. 887 888 889 **Materials and methods** 890 Plant material for genetic mapping - the TS collection 891 Seeds of the collection used for GWAS and the genomic variation information for the 892 accessions were kindly provided by Harry Klee (University of Florida). The collection 893 includes 15 _S._ _pimpinellifolium,_ 83 _S._ _lycopersicum_ var. cerasiforme and 300 diverse _S._ 894 _lycopersicum_ tomato varieties including many heirloom and modern cultivars, of which 895 we used 362 accessions. The population was previously genotyped and the population 896 structure has been studied and published [2] . This population was used for mapping of 897 metabolic traits related to fruit flavor [6] . 898 Greenhouse experiment setup and treatment regime 899 Seeds of all accessions of the TS collection were sown in seedling trays and placed 900 in a commercial nursery for germination and establishment. Four weeks after 901 germination, the seedlings were transferred to the experimental greenhouse. The 41 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 902 experiment was designed in complete randomized blocks, with one plant of each 903 accession in each block. Three complete randomized blocks of the population were 904 grown in control conditions, akin to commercial growing conditions. Three additional 905 blocks received water deficit (WD) treatment. After two weeks in the greenhouse, 906 irrigation for the WD blocks was gradually reduced until there was a measurable 907 difference in stress measures between WD and control plants (see section below). At 908 this point, WD plants received 12% of the water given to control plants, and this rate 909 was maintained for the remainder of the experiment. Exact watering amounts during the 910 full experiment are reported elsewhere [125] . 911 912 Leaf water content and vegetative production 913 Leaves from one block of both control and WD were evaluated for water content 914 three weeks after initiation of the WD treatment. The full fourth open leaf was picked 915 and immediately weighed at the greenhouse using a semi-analytical scale. Leaves were 916 then oven dried at 60 °C for 24 h then weighed with an analytical scale. The water 917 content was then calculated by: (fresh weight - dry weight)/fresh weight. 918 At the end of the growing season, after all fruits were harvested, the stem and leaf 919 biomass was weighed (fresh weight) to estimate vegetative production. 920 921 Leaf sampling for metabolomics and transcriptomics 922 Leaf sampling was conducted during a four-hour period in the middle of the day 923 (10:00 - 14:00), to reduce variation due to diurnal metabolic cycles. For metabolic 924 analysis, two independent complete blocks from each condition were sampled and 925 analyzed separately to evaluate reproducibility. On each sampling day, a pair of control 926 and WD blocks were sampled to ensure sampling at appropriate hours. 927 For each sampled plant three distal leaflets from the third, fourth and fifth leaves 928 were pooled together in a vial and immediately flash frozen in liquid nitrogen. Samples 929 were kept completely frozen and stored at -80 °C until extraction. 930 Samples were ground to a fine powder using a Mixer Mill-MM 200 (Retsch), and 100 931 mg of frozen powder was weighed for extraction. 42 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 932 Data processing 933 All data processing was done in R (version 4.1.3) [126] in the RStudio environment 934 (2022.07.1 build 554) using packages as described or custom code available as listed 935 under Data availability. 936 Transcriptomic profiling 937 RNA extractions 938 Of the same frozen leaf powder used for metabolic extraction, 50 mg frozen powder 939 were weighed for RNA extraction. The extraction followed the TRIzol protocol (Merck, 940 TRI Reagent®, Cat no. 93289) with some modifications [127] . Briefly, 1 ml chilled TRIzol 941 was added to each frozen sample, it was then vortexed and placed on ice. After 5 min 942 incubation at room temperature samples were centrifuged 8 minutes at 24000 g at 4 °C. 943 Then chloroform separation was performed twice: the clear upper phase was 944 transferred to a new tube, 0.25 ml chloroform was added, samples were vortexed, 945 incubated 15 minutes at room temperature, centrifuged 8 minutes at 24000 g at 4 °C, 946 and 0.4 ml of the aqueous phase was transferred to a new tube. To that 0.3 ml chilled 947 isopropanol was added and samples were gently mixed and placed at -20 °C overnight. 948 Samples were centrifuged 25 minutes at 24000 g at 4 °C to pellet RNA. The pellet was 949 washed twice with 75% ethanol. Then it was re-suspended in 100 µl double distilled 950 water and 10 µl sodium acetate (3M, pH=5.2) before 0.3 ml chilled absolute ethanol was 951 added and samples were placed at -20 °C overnight. Following, RNA was pelleted and 952 washed twice as described above. Finally, the clean RNA pellet was re-suspended in 30 953 µl double distilled water. RNA quality was evaluated by agarose gel separation stained 954 with SYBR™ Safe (Invitrogen™, cat no. S33102) and Nanodrop (DeNovix) reading of 955 concentration (260nm) and purity (260/280 and 260/230 >1.8). The concentration of all 956 samples was adjusted to 18 ng/µl for TranSeq library preparation. 957 958 RNA sequencing by TranSeq 959 The detailed method of library preparation, sequencing and results processing of the 960 TranSeq method was described by Tzfadia et al. [127] . Briefly, RNA was fragmented using 43 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 961 NEB fragmentation buffer (NEB cat. no. E6150S) at 94 °C for 4 minutes. After quick 962 cooling to -20 °C, 3' ends were pulled down by poly-T beads. Selected poly-A containing 963 fragments were barcoded and converted to cDNA by barcode-poly-T primers for reverse 964 transcription. After barcoding, cDNA samples were pooled before RNAse treatment. 965 Then cDNA was cleaned from reagents with SPRI beads (Beckman, AMPure XP) and 966 PCR adaptors were ligated to the cDNA. Following a second SPRI bead cleaning, cDNA 967 samples were amplified by PCR. Quality of resulting library was evaluated by Qubit 968 (Invitrogen™, Qubit™ dsDNA HS Cat no. Q32851) for concentration reading and 969 TapeStation (Agilent) for evaluating library size distribution. Six libraries, each pooled 970 from 94 RNA samples (including three quality controls in each library) were sequenced 971 by 100 bp single reads on a NovaSeq machine (Illumina). Reads were aligned to tomato 972 genome version SL2.50 ITAG2.4. 973 974 TranSeq - filtering, batch effect normalization and eQTL 975 Raw counts of each gene were normalized using vst in the R package DESeq [128] . 976 Genes absent from more than 90% of the samples were removed. Genes which had 977 median count in the quality control samples lower than 3.4 (empirically determined 978 threshold) were also removed. In order to reduce variation caused by batches of library 979 preparation and sequencing, each gene's reads were normalized to the median of 980 batches from the same treatment. The resulting normalized reads were evaluated by 981 PCA to assess batch effect and check normalization of quality control samples using 982 pcaMethods [129] . 983 Normalized gene counts were used for GO term enrichment, feature-gene 984 correlations and expression QTL (eQTL). eQTL were obtained by using gene 985 expression as phenotype for GWAS mapping with the same method described below for 986 mQTLs. 987 988 Gene Ontology (GO) term enrichment 989 For each gene expressed in the leaf samples, expression was compared between 990 WD and control samples across the full population. The fold change and significance 44 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 991 (Student's two-tailed t-test) of the change in WD compared to control were calculated in 992 R. The genes with p value lower than Bonferroni corrected α=0.05 and fold change 993 higher than 1.5 were considered up regulated. Conversely, genes with p value lower 994 than Bonferroni corrected α=0.05 and fold change lower than 0.5 were considered down 995 regulated. 996 The lists of genes either up or down regulated were subjected to GO term 997 enrichment analysis according to the instructions in the Gene Ontology website 998 [http://geneontology.org/docs/go-enrichment-analysis using Fisher's exact test with FDR](http://geneontology.org/docs/go-enrichment-analysis) 999 correction. The full list of expressed genes found in our experiment was used as a 1000 reference list. 1001 1002 Metabolic profiling 1003 Metabolomic extraction 1004 Extraction protocol was previously described [130] . In brief: 1 ml pre-cooled (−20 °C) 1005 extraction buffer (homogenous methanol/methyl-tert-butyl-ether [1:3] mixture) was 1006 added to frozen plant sample. After 10 min incubation at 4 °C and sonication for 10 min 1007 in an ultra-sonication bath, 500 μl of methanol/water [1:3] mixture was added. Samples 1008 were vigorously mixed by vortex, then centrifuged (5 min, 14,000 g), resulting in 1009 separation of the lipophilic and polar phases. Phases were aliquoted for analysis in 1010 three chromatography-mass spectrometry platforms. Three hundred microliters of the 1011 (upper) lipophilic phase were collected for analysis of lipophilic compounds. Of the polar 1012 phase, 300 µl were taken for analysis by LC-MS. All aliquoted samples were fully dried 1013 under a vacuum and stored at -80 °C until resuspension. Prior to LC-MS analysis dried 1014 samples were re-suspended with 100 μl methanol water (1:4) mixture. 1015 LC-MS profiling- polar metabolites 1016 LC-MS data were obtained using a Waters Acquity UPLC system (Waters) fitted with 1017 a C18 column, coupled to a Q-Exactive Orbitrap mass spectrometer (Thermo Fisher 1018 Scientific). Instrumental settings for quantification in full MS were previously described 1019 by Giavalisco _et al._ [131] . In detail, the UPLC separation was performed using an HSS T3 1020 C18 reversed-phase column (100 mm × 2.1 mm × 1.8 μm particles; Waters). The mobile 45 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1021 phases were 0.1% formic acid in H2O (BufferA, UPLC grade; BioSolve, 1022 http://www.biosolve-chemicals.com) and 0.1% formic acid in acetonitrile (BufferB, UPLC 1023 grade; BioSolve). A 2 μl re-suspended sample was loaded per injection. The gradient 1024 flow rate was 400 µl/min, and consisted of: 1 minute 99% A, 1-11 minute linear gradient 1025 from 99% A to 60% A, 11-13 minutes linear gradient from 60% A to 30% A, 13-15 min 1026 linear gradient from 30% to 1% A, hold 1%A until 17 min, 17-17.5 minutes linear 1027 gradient from 1% to 99% A, and re-equilibrate the column for 2.5 min (20-min total run 1028 time). 1029 The UPLC-separated sample was directly injected into a Q-Exactive mass 1030 spectrometer. Each sample was analyzed in both positive and negative polarities via 1031 separate, consecutive LC-MS runs. The spectra for quantification of all samples were 1032 recorded full-scan and covering a mass range from 100 to 1500 m/z. The resolution was 1033 set to 70000, with 3 scans per second, restricting the loading time to 100 ms. The 1034 capillary voltage was set to 3.5 kV with a sheath gas flow value of 60 and an auxiliary 1035 gas flow of 20 (values are in arbitrary units). The capillary temperature was set to 275 1036 °C, whereas the drying gas in the heated electrospray source was set to 350 °C. The 1037 skimmer voltage was set to 25 V, whereas the tube lens was set to a value of 130 V. 1038 The spectra were recorded from 1 min to 20 min of the UPLC gradients. 1039 Chromatogram processing, quantification, filtering, feature consolidation and 1040 normalization 1041 The Xcalibur Qual browser (Thermo Fisher Scientific) was used to visualize and 1042 explore chromatograms. 1043 REFINER MS 11 (GeneData) software was used for chromatogram processing and 1044 quantification. Processing of MS data included peak alignment, removal of isotopic 1045 peaks, and chemical noise. Peak identifiers (Cluster ID), retention time (RT) and m/z 1046 were obtained along with peak area in each of the samples. The peak area was chosen 1047 to quantify the relative abundance of corresponding features in the samples. To reduce 1048 the redundancy, all isotopic peaks from the same compound were clustered, and each 1049 cluster, represented by its monoisotopic m/z, is referred to as a metabolite feature 1050 (hereafter feature). Several features can derive from one compound. 46 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1051 Features with median area intensity lower than 100,000, or which were present in 1052 the blank negative controls with median intensity higher than 1,000,000, were removed. 1053 Features found in less than 30% of the samples in either condition were also removed. 1054 The resulting list of features was consolidated in order to remove redundant 1055 information caused by adducts or possible multiple detected fragments (in-source 1056 fragmentation) of the same compound, using a custom R script. The algorithm is based 1057 on similar RT and high correlation (r>0.9) between features across samples. 1058 The resulting feature abundance values were log2 transformed and normalized to 1059 reduce machine-caused variance using each sample's median intensity and median 1060 machine run batch intensity, using R. The process of leaf sampling, extraction, LC-MS 1061 analysis, data processing, filtering and normalization was done separately for two 1062 repetitions of the population grown in control and WD conditions. 1063 Annotation by 30 representative samples 1064 For detailed annotation of all detected features, 30 representative samples were 1065 subjected to multiple methods of annotation. The selected samples were genetically 1066 varied from both conditions. Extracted aliquots from the samples were sent for analysis 1067 in the lab of Prof. Aharoni for high confidence annotation using WeizMass [47] and lab in 1068 house libraries. The resulting annotations were compared to ddMS [2] fragmentation 1069 spectra obtained from the same samples by methods described in the next section. 1070 LC-MS for ddMS [2] and list-based ddMS [2] acquisition 1071 For high-level metabolite annotation, it was crucial to obtain feature-specific 1072 fragmentation spectra. For high coverage ddMS [2] fragmentation of all features, the 30 1073 representative samples were analyzed by LC-MS, using the exact same 1074 chromatographic method previously used but with ddMS [2] Top10 method in the MS. In 1075 detail: MS acquisition was during 0.4-16.55 min of the full chromatography run, either 1076 polarity was acquired separately. For full MS: resolution - 70000, AGC target - 1e6, max 1077 IT - 200 ms, scan range - 100-1500 m/z. For ddMS [2] : resolution - 17500, AGC target 1078 1e5, max IT - 50 ms, loop count - 10, isolation window - 3.0 m/z, CE - 25, 35 (stepped), 1079 spectrum data type - centroid. For dd settings: minimum AGC - 2.5e3, apex trigger - 1 1080 to 4 s, charge exclusion - unassigned, dynamic exclusion - 5.0 s. 47 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1081 Despite the high coverage, not all features mapped by GWAS had ddMS [2] spectra 1082 using the Top10 method. Therefore, features associated with detected mQTLs were 1083 reanalyzed in the LC-MS, this time with a ddMS [2] inclusion list method. A list of features, 1084 and their RT and m/z identifiers were included for specific fragmentation. In detail: MS 1085 acquisition was during 0-17.58 min of the full chromatography run. For full MS, either 1086 polarity was acquired separately: resolution - 70000, AGC target - 3e6, max IT - 200 ms, 1087 scan range - 100-1500 m/z. For parallel reaction monitoring (PRM): resolution - 17500, 1088 AGC target - 5e5, max IT - 100 ms, isolation window - 3.0 m/z, NCE - 25, 35 (stepped). 1089 Compound Discoverer 1090 To maximize annotation methods, the 30 samples were analyzed with Compound 1091 Discoverer 3.3 (Thermo Fisher Scientific). 1092 Processing parameters: All spectra between 0.45-16.5 min were acquired. 1093 Chromatogram RT was aligned according to the quality control pool sample. For detect 1094 compounds: mass tolerance 5 ppm, min peak intensity - 100000, S/N threshold - 1.5, no 1095 baseline removal, Zig-Zag index threshold - 0.2, jaggedness threshold - 0.4, modality 1096 threshold - 0.9, considered ions- [2M+FA-H]-1; [2M+H]+1; [2M+K]+1; [2M+Na]+1; 1097 [2M+NH4]+1; [2M-H]-1; [2M-H+HAc]-1; [M+FA-H]-1; [M+H]+1; [M+H+H2O]+1; [M+K]+1; 1098 [M+Na]+1; [M+NH4]+1; [M-H]-1; [M-H+H2O]-1; [M-H+HAc]-1; [M-H-H2O]-1, base ions 1099 [M+H]+1; [M-H]-1 . For group compounds: mass tolerance 5 ppm, RT tolerance 0.1 min, 1100 align peaks - false, preferred ions: [M+FA-H]-1; [M+H]+1; [M+Na]+1; [M+NH4]+1; [M-H] 1101 1; [M-H+HAc]-1, area integration - most common ion, peak rating contribution was: area 1102 -3, CV -10, FWHMN -5, Zig-Zag index - 5, jaggedness-5, modality threshold - 5, peak 1103 rating threshold - 4, number of files -1. For fill gaps: mass tolerance 5 ppm, S/N 1104 threshold - 1.5. For QC correction the SERRF algorithm was used: min QC coverage 1105 50%, max QC area RSD - 30%, max corrected QC area RSD 50%, max # files between 1106 QC files - 15, # batches -1, interpolate gap filled QC areas - false. 1107 Detected features were compared to several libraries for annotation: full MS m/z 1108 were compared to ChemSpider [132], an in-house-run central metabolites chemical 1109 standards library (MSMLS, Mass Spectrometry Metabolite Library of Standards, Sigma 1110 Aldrich), and literature-derived lists [20,133–136] as mass lists, ddMS [2] fragmentation spectra 1111 were compared to mzCloud online database [137] and multiple acquired databases: 48 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1112 Bamba laboratory 598 polar metabolites, Fiehn Vaniya natural product library, 1113 MassBank of North America (MoNA) (export 2022), and ReSpect [138] . Priority in 1114 annotation was given to in-house-run standards, ddMS [2] based matching, followed by 1115 precursor-ion-based mass-list matching (Supp. Table 3). 1116 GNPS and feature-based molecular networking 1117 Global similarity of all detected features was used for extending annotation from 1118 annotated features to similar ones. Chromatograms from the 30 representative samples 1119 were used, analyzed with the regular full MS method and the Top10 ddMS [2] method, in 1120 either positive or negative mode. First, chromatogram files were converted from raw to 1121 mzML format, using ProteoWizard file converter (v. 3.0.21112) [139] . Next, chromatogram 1122 and spectra, for each polarity separately, were aligned and processed using MS-DIAL 1123 (v. 4.60) [140], with the following specifications: Ionization type - Soft ionization, Separation 1124 type - Chromatography, MS method type - Conventional LC/MS, Data type (MS1) 1125 Centroid, Data type (MS/MS) - Centroid, Ion mode - Positive/Negative, Target omics 1126 Metabolomics. For Data collection: Mass accuracy tolerance for MS1- 0.01 Da, for MS2 1127 0.025 Da, RT range - 0.5-16.5, MS1 mass range - 100-2000 Da, MS/MS mass range 0 1128 2000 Da, maximum charge - 2. For peak detection: minimum peak height 1000 1129 amplitude, mass slice width - 0.1 Da, smoothing method - Linear weighted moving 1130 average, smoothing level - 3 scans, minimum peak width - 5 scans. For MS2 1131 Deconvolution: sigma window value - 0.5, keep isotopic ions until 0.5 Da. For 1132 identification: RT tolerance - 0.5 min, accurate mass tolerance for MS1 - 0.01 Da, for 1133 MS2 - 0.05 Da, and the rest as default. Adducts considered in positive polarity: [M+H]+, 1134 [M+NH4]+, [M+Na]+, [M+CH3OH+H]+ and [M+K]+. Adducts considered in negative 1135 polarity: [M-H]-, [M-H2O-H]-, [M+Na-2H]-, [M+FA-H]-, [M+Hac-H]-, [2M-H]-, [2M+FA-H] 1136 and [2M+Hac-H]-. For alignment: reference file was set as the quality control pool 1137 sample, RT tolerance - 0.05 min, MS1 tolerance - 0.015 Da, and the rest as default. 1138 Alignment results were exported for "GNPS export" in the "mgf" format. 1139 In the Global Natural Products Social Molecular Networking (GNPS) [46] module for 1140 FBMN (v. 28.2) a new project was started for each polarity. The MS-DIAL alignment 1141 table, the spectra file in mgf format and all chromatogram files in mzML format were 1142 uploaded using an ftp client, WinSCP (v. 5.17.10). FBMN was run with default settings, 49 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1143 except Quantification Table Source was set to MS-DIAL, Precursor Ion Mass Tolerance 1144 0.75 Da, Fragment Ion Mass Tolerance - 0.5 Da. 1145 Results were visualized within the GNPS browser and exported in "Cytoscape data" 1146 format for further analysis. Cytoscape (v3.9.1) [141] was used for further network 1147 exploration. 1148 CFM-ID in silico fragmentation prediction 1149 Predicted fragmentation patterns of formula proposed by Xcalibur (Thermo Fisher 1150 Scientific) were generated using CFM-ID 3.0 online software [142] . 1151 Calculation of expected mass from molecular formula 1152 While comparing predicted compositions and library annotations, molecular formulas 1153 were used to calculate predicted average molecular weight and monoisotopic mass by 1154 R package "OrgMassSpecR" [143] . 1155 GWAS mapping with GAPIT 1156 The normalized relative abundance of each feature in each condition was used as 1157 the phenotype for mQTL mapping by GWAS. Population single nucleotide 1158 polymorphism (SNP) variations were kindly provided by Harry Klee (University of 1159 Florida). The GWAS SNP locations and all genome information are based on SL2.50 [1] in 1160 SGN browser (https://solgenomics.net/jbrowse_solgenomics/?data=data/json/SL2.50), 1161 and gene predictions are based on ITAG2.4, unless otherwise stated. 1162 GWAS analysis was performed using the Genome Association and Prediction 1163 Integrated Tool (GAPIT) algorithm implemented in R [144,145] . Generally, the default 1164 settings were used, specifically the Mixed Linear Model (MLM) method was used and 1165 VanRaden kinship [144,146] and three principal components were used to control for 1166 population structure. 1167 1168 mQTL co-location grouping and comparisons 1169 After gathering the results from mapping of all features in all conditions, co-locating 1170 mQTLs of all features were identified using a custom script. mQTLs were regarded as 1171 one co-locating mQTL if significant SNPs (-log(p)>8.3, Bonferroni correction for α=0.05 50 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1172 for 1.8 million SNPs) were at the most 100 kb apart. The set of features associated with 1173 one mQTL is referred to as its mQTL feature group. A feature can belong to more than 1174 one mQTL feature group (mQTLs on different chromosomes). Comparison of mQTL 1175 feature groups and significance level at each condition was performed. 1176 To evaluate the reproducibility of the mQTL mapping, all features which had 1177 significant mapping in either condition of the first control-WD block pair were listed. The 1178 corresponding features in the second block pair's dataset were identified by matching 1179 RT and m/z identifiers. The relative abundance of the matched features in the second 1180 pair was mapped with GAPIT, as described above, and resulting mQTLs were 1181 compared for co-mapping in either condition as described above (but with significance 1182 threshold of -log(p)>7). 1183 1184 Commonality analysis 1185 The mgf files made by MS-DIAL (v. 4.60) [140] from all experiments were combined in 1186 R, using the "MsBackendMgf" package [147] . The ddMS2 spectra were processed with the 1187 "CluMSID" package [148] and custom scripts. Redundant spectra from features with close 1188 RT within 0.05 min and m/z within 0.005 Da, and spectra similarity cosine score above 1189 0.6 were removed. The many fragments with negligible counts are a heavy burden on 1190 downstream calculations, therefore, fragments with counts lower than 2% of the 1191 highest-count fragment were removed. All spectra with less than 2 fragments were 1192 removed, since they could not be used in downstream calculations. 1193 The commonality, the set of descriptor frequencies in a group, was calculated using 1194 custom scripts. In brief, all ddMS [2] spectra in the defined group were assembled. 1195 Fragments were binned with an m/z tolerance of 0.005 Da, and the number of detected 1196 fragments in all group spectra was calculated. The frequency of a fragment is the 1197 proportion of the features with ddMS2 spectra in the group whose spectrum contains it. 1198 The frequency of the top 50 fragments in the ddMS [2] group was plotted as a 1199 commonality spectrum. The frequency of each fragment in the commonality of all ddMS [2] 1200 spectra was defined as the expected frequency for further statistical analysis. The G 1201 test [149] was used to statistically compare commonality of tested groups to the expected 51 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1202 frequency. We defined the statistical threshold of significance based on 2,000 randomly 1203 sampled ddMS [2] spectra groups. We found that the G-test power was greatly affected by 1204 the size of the group (Supp. Note 1), therefore, subsequent tests were compared to the 1205 threshold based on similar group size to determine significance. A group whose 1206 commonality differed significantly from the expected commonality by this test is 1207 described as having a significantly distinct commonality. Separately, a fragment whose 1208 frequency in a group exceeded the mean top frequency of size-matched random groups 1209 is described as an outstanding descriptor. This is a descriptive threshold, not a 1210 significance test. High-impact mQTLs were defined as mQTLs with enough features to 1211 statistically test a descriptor type in one ionization mode: at least 20 ddMS [2] spectra for 1212 fragments, 50 ddMS [2] spectra for NLs, and 50 features for transformations (Supp. Note 1213 1). 1214 The neutral losses (NLs) for each ddMS [2] spectrum were calculated as the 1215 differences between all fragments, including the full m/z of the precursor ion, rounded to 1216 0.01 Da. All commonality and statistical analysis were then conducted in the same way 1217 as described above for fragments. NL annotation was based on a list curated from the 1218 literature (Supp. Table 7). 1219 Transformation calculations used the m/z of the precursor ions. A transformation is 1220 the m/z difference between the precursor ions of two features in a group. It is a mass 1221 difference and does not by itself demonstrate a biochemical conversion. For general 1222 transformation commonality the differences between all m/z were binned with tolerance 1223 of 0.01 Da and the frequency of each difference was calculated out of the number of 1224 non-redundant pairs and plotted. Transformations were compared to calculations from 1225 the MetNet R package [31] . For directed transformation commonality, only pairwise m/z 1226 differences between opposing groups were included. The expected transformation 1227 frequencies were calculated based on m/z differences of all detected precursor ions. 1228 The G-test [149] was used to statistically compare transformation commonality of tested 1229 groups to the expected frequency. We defined the statistical threshold of significance 1230 based on 2,000 random permutations of randomly sampled precursor ion m/z values. 1231 We found that, here too, the G-test power was greatly affected by the size of the group 52 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1232 (Supp. Note 1), therefore, subsequent tests were compared to the threshold based on 1233 similar group size to determine significance. Transformation annotation was based on a 1234 list curated from the literature (Supp. Table 7). 1235 1236 Feature and gene expression data analysis 1237 Principal component analysis (PCA) 1238 PCA was used to gain initial overview of the samples' metabolomic profile and 1239 evaluate normalization performance, using the R package pcaMethods [129] . 1240 Feature-feature correlations 1241 The features of each mQTL feature group were evaluated by correlations. 1242 Correlation matrices were made with the R package corrplot [150] . Features were grouped 1243 according to hierarchical clustering and major groups were assigned a letter from A-C. 1244 Correlation matrices were made for each condition separately. 1245 Feature-gene correlations 1246 For each mQTL, a correlation matrix of the relative abundance of all features in the 1247 mQTL feature group to the normalized expression of all genes in the mQTL region was 1248 produced. Correlations with p-value <0.05 were plotted. Correlation matrices were made 1249 for each condition separately. Correlations with p value <10 [-7] and |r|>0.4 were 1250 considered significant. 1251 Phenotype single allele interactions 1252 To look at the effect of the mQTL alleles on feature abundance, the genotype of the 1253 mQTL lead SNP in each accession was extracted, and the effect of the allele on feature 1254 abundance in the samples was tested by Student's t-test. 1255 Linkage disequilibrium (LD) 1256 To assess LD of SNPs in the mQTL region, all SNPs in the mQTL region and 1257 extending 5 kb before and after the mQTL were used. Calculations and plotting were 1258 performed with the R package LDheatmap [151] . 53 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1259 Isotope labeling 1260 15N isotope labeling 1261 To provide the seeds the material required for their growth to mature plants, 20 S. 1262 lycopersicum seeds were germinated on 1/2 Murashige and Skoog (MS) medium (2.15 1263 g/L) and 3% glucose for 30 d. To label the nitrogen within the plant all nitrogen 1264 components within the MS medium were replaced with labeled alternatives: 1650 mg/L 1265 Ammonium nitrate ( [15] NH 4 [15] NO 3, Sigma-Aldrich: 366528) and 1900 mg/L Potassium 1266 nitrate (K [15] NO 3, Merck: 57654-83-8) in the labeled media. The MS medium was then pH 1267 balanced to 5.8 ± 0.1 before adding the glucose and autoclaving for 20 min at 121 °C. 1268 The plants grew 30 days in 16-8 light conditions and 24 °C before mature leaves were 1269 harvested and snap-frozen using liquid nitrogen for UPLC-MS/MS metabolite analysis. 1270 13C isotope labeling 1271 For 13C isotope labeling _S. lycopersicum_ calli were grown in tissue culture from 1272 lesions in mature leaves in 0.8% agar MS medium supplemented with 0.1 mg/L indole 1273 3-acetic acid (IAA) for 30 days in 16-8 light conditions and 24 °C, according to 1274 Moscatiello _et al._ [152] . The resulting callus cells were then roughly homogenized using a 1275 pipette and grown in media consisting of 4.3 g/L MS medium, 1 μg/mL 1 1276 naphthaleneacetic acid (NAA), and 0.1 μg/mL of both 2,4-dichlorophenoxyacetic acid 1277 and kinetin. 1278 To label carbon, the medium was supplemented with 3 g/L of either C 6 H 12 O 6 glucose 1279 (control), or [13] C 6 H 12 O 6 (Sigma: 389374). These cell cultures were suspended in a 1280 shaking incubator at 25 °C for 32 d before harvesting. The cells were harvested 1281 according to a UTHealth Houston cell sample preparation protocol. To process the cells 1282 for UPLC-MS/MS the samples were centrifuged for 2 min at 1,500 rpm and the 1283 supernatant discarded. 1 mL PBS was then added to the pellets, which were then 1284 roughly ground using a pipette and centrifuged for 2 min at 2,000 rpm before removing 1285 PBS. Afterwards, samples were snap-frozen using liquid nitrogen and stored at -80 °C 1286 until UPLC-MS analysis. 54 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1287 Isotope labeled metabolomic detection and analysis 1288 The labeled samples were extracted and prepared and run on the UPLC-MS in the 1289 same way as described for leaf samples. Chromatograms were visualized in Xcalibur 1290 Qual browser (Thermo Fisher Scientific). For each feature first the expected peak was 1291 identified by matching m/z and RT in control samples. Next, the labeled samples were 1292 tested for presence of a peak with the expected m/z and absence or reduction of the 1293 original unlabeled m/z. Elemental compositions of unannotated features were estimated 1294 in Xcalibur. Expected m/z values for either [13] C or [15] N labeled compounds were 1295 calculated using https://bmrb.io/metabolomics/mol_mass.php. 1296 Putative gene annotations 1297 Putative gene annotations were initially taken from ITAG 2.4, and investigated 1298 further by BLAST (TAIR, NCBI) to find homologous genes or motifs with known 1299 functions. 1300 Geneious 1301 All sequence analysis, viewing, alignment, neighbor-joining trees and BLAST 1302 searches, were done in Geneious Prime® 2022.1.1 Software, unless otherwise stated. 1303 Validation by overexpression (OE) 1304 Cloning of target genes for OE 1305 Genes of interest were amplified from leaf cDNA using HS TAQ MIX RED (PCR 1306 BIOSYSTEMS, London, UK) or PrimeSTAR HS DNA Polymerase (TAKARA, Shiga, 1307 Japan) followed by digestion with KpnI and XbaI restriction enzymes (New England Ltd., 1308 Hitchin, UK) and ligation using T4 DNA Ligase (New England Ltd., Hitchin, UK) into the 1309 multiple cloning site of the _pBTEX_ plasmid driven by the CaMV 35S promoter. All primers 1310 used in this project are listed in Supp. Table 13. Plasmids were introduced into 1311 _Agrobacterium tumefaciens_ strain _GV3101_ by heat shock. Cotyledons isolated from _S._ 1312 _lycopersicum cv. M82_ (for Solyc02g071100 OE) were incubated with _A. tumefaciens_, 1313 harboring the _pBTEX::35S::gene_ vectors following the protocol described by 1314 McCormick [153], to generate transgenic OE lines. Transformation of Ailsa Craig (for 55 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1315 Solyc12g088250 and Solyc04g077650 OE) and Micro-Tom (Solyc09g091700 OE) was 1316 performed as described by Ouyang _et al_ . [154] . 1317 RNA extraction, cDNA synthesis and qPCR 1318 RNA was extracted from leaves using the Total RNA Mini Kit (Plant) from Geneaid 1319 (New Taipei, Taiwan). The success of the extraction and the integrity of the RNA was 1320 validated using NanoDrop (MaestroNano, New Taipei, Taiwan) and gel electrophoresis. 1321 cDNA strands were synthesized using a designated kit (qScript cDNA Synthesis Kit; 1322 BIOSYSTEMS, London, UK) and PCR (SimpliAmp™ Thermal Cycler; ThermoFisher 1323 Scientific, Waltham, MA, USA), according to the manufacturer's instructions. 1324 Quantitative PCR was performed using the PikoReal 96 Real-Time PCR System 1325 (Thermo Scientific, MA, USA) and qPCRBIO SyGreen Blue Mix Hi-ROX kit by PCR 1326 Biosystems (London, UK). The different sets of primers used for the amplification of the 1327 target genes are listed in (Supp. Table 13). Analysis of the relative gene expression was 1328 performed according to the comparative cycle threshold (2−ΔΔCT) method [155], calibrated 1329 using transcript values relative to the endogenous tomato actin gene. 1330 OE growing conditions, WD 1331 For metabolomic sampling, WT and OE lines were grown in the greenhouse under 1332 optimal irrigation. The fourth true leaf was sampled four weeks after transplanting 1333 seedlings to the greenhouse. For Solyc02g071100, OE plants were grown in optimal 1334 conditions (100% of previous day transpiration) for two weeks, then irrigation was 1335 dropped to 50% of daily transpiration for half of the plants for 10 days. Leaf samples 1336 were taken prior to applying the WD regime and 4 and 9 days after, from water-limited 1337 and fully irrigated (control) plants. 1338 1339 1340 56 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1341 **Author contributions/Acknowledgment** 1342 Author contributions 1343 L.R. designed and managed experiments, prepared samples, analyzed the data and 1344 wrote and revised the manuscript. Y.Y.Z., N.Sa. and A.B-S. prepared overexpression 1345 lines and conducted their experiments and sample preparation. T.L-C. and S.W. 1346 designed and managed experiments, prepared samples and metabolic data acquisition. 1347 C.T. prepared samples and data acquisition for isotope labeling experiments. E.S. 1348 prepared samples and developed chromatogram processing and metabolite annotation 1349 analysis pipelines and reviewed the manuscript. T.N. and S.P. prepared samples for 1350 transcriptomics. I.R., N.Sh., and A.A. acquired data and analyzed samples for 1351 metabolite annotation. B.O. and S.D. generated transgenic lines. Y.B. initiated the 1352 project, acquired funding, supervised all experiments and wrote and reviewed the 1353 manuscript. All authors read and approved the manuscript. 1354 1355 Acknowledgements 1356 We would like to acknowledge the immense support we received from the Ben 1357 Gurion University's Ilse Katz Institute for Nanoscale Science & Technology Mass 1358 Spectrometry and chromatography unit, Elena Doubijansk, Dr. Mark Karpasas, Dr. 1359 Albert Batushansky and Dr. Efrat Forti. This work would not be possible without you. We 1360 also thank Dr. Albert Batushansky for critical review of the manuscript. Our thanks are 1361 extended to Shai Shaked and the greenhouse workers for the high-quality professional 1362 cultivation of the tomato GWAS experiment. 1363 We would also like to thank Prof. Dr. Dirk Walther for discussing the statistical 1364 framework for the commonality analysis. 1365 1366 Funding 1367 This research was supported by grant number 16-37-0002 from the Chief Scientist 1368 of the Ministry of Agriculture and Rural Development, Israel, awarded to the Brotman 57 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1369 Lab for the project titled "Identifying the genetic components of the metabolic response 1370 to water deficit in tomato." 1371 Leah Rosental was a fellow of the Ariane de Rothschild Women's Doctoral Program 1372 and recipient of the scholarship for Ph.D. students (Ben Gurion University of the Negev) 1373 and Kreitman Fellowship for short-term post-doctoral research. She is currently 1374 supported by the Marie Skłodowska-Curie Actions Postdoctoral Fellowship (Project 1375 101153815, CPRrescue). 1376 Bo Ouyang is supported by National Natural Science Foundation of China 1377 (32472743) and the National Key Research and Development Program of China 1378 (2022YFE0100900). 1379 1380 Conflicts of interest 1381 The authors declare no conflict of interest 1382 AI disclosure 1383 AI-assisted language editing was used in preparation of the manuscript. No 1384 generative AI was used to write new content. 1385 **References** 1386 1. Consortium, T. T. G. The tomato genome sequence provides insights into fleshy 1387 fruit evolution. _Nature_ **485**, 635–41 (2012). 1388 2. Lin, T. _et al._ Genomic analyses provide insights into the history of tomato 1389 breeding. _Nat. Genet._ **46**, 1220–1226 (2014). 1390 3. Gao, L. _et al._ The tomato pan-genome uncovers new genes and a rare allele 1391 regulating fruit flavor. _Nat. Genet._ **51**, 1044–1051 (2019). 1392 4. Zhou, Y. _et al._ Graph pangenome captures missing heritability and empowers 1393 tomato breeding. _Nature_ **606**, (2022). 1394 5. Sauvage, C. _et al._ Genome-wide association in tomato reveals 44 candidate loci 1395 for fruit metabolic traits. _Plant Physiol._ **165**, 1120–1132 (2014). 58 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1396 6. Tieman, D. _et al._ A chemical genetic roadmap to improved tomato flavor. _Science_ 1397 _(80-. )._ **355**, 391–394 (2017). 1398 7. Li, X., Tieman, D., Liu, Z., Chen, K. & Klee, H. J. Identification of a lipase gene 1399 with a role in tomato fruit short-chain fatty acid-derived flavor volatiles by genome 1400 wide association. _Plant J._ **104**, 631–644 (2020). 1401 8. Zhu, G. _et al._ Rewiring of the Fruit Metabolome in Tomato Breeding. _Cell_ **172**, 1402 249-261.e12 (2018). 1403 9. Szymański, J., Bocobza, S., Panda, S., Sonawane, P. & Cárdenas, P. D. Analysis 1404 of Wild Tomato Introgression Lines Elucidates the Genetic Basis of Transcriptome 1405 and Metabolome Variation Underlying Fruit Traits and Pathogen Response. _Nat._ 1406 _Genet._ **52**, (2020). 1407 10. Nunes-Nesi, A. _et al._ Identification and characterization of metabolite quantitative 1408 trait loci in tomato leaves and comparison with those reported for fruits and seeds. 1409 _Metabolomics_ **15**, 1–13 (2019). 1410 11. Manuela M. Chaves, Maroco, J. P. & Pereira, J. S. Understanding plant 1411 responses to drought — from genes to the whole plant. _Funct. Plant Biol._ **30**, 1412 239–624 (2003). 1413 12. Claeys, H. & Inzé, D. The agony of choice: How plants balance growth and 1414 survival under water-limiting conditions. _Plant Physiol._ **162**, 1768–1779 (2013). 1415 13. Illig, T. _et al._ A genome-wide perspective of genetic variation in human 1416 metabolism. _Nat. Genet._ **42**, 137–141 (2010). 1417 14. Fang, C. & Luo, J. Metabolic GWAS-based dissection of genetic bases underlying 1418 the diversity of plant metabolism. _Plant J._ **97**, 91–100 (2019). 1419 15. Chen, J., Xue, M., Liu, H., Fernie, A. R. & Chen, W. Exploring the genic resources 1420 underlying metabolites through mGWAS and mQTL in wheat: From large-scale 1421 gene identification and pathway elucidation to crop improvement. _Plant Commun._ 1422 **2**, 100216 (2021). 59 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1423 16. Visscher, P. M. _et al._ 10 Years of GWAS Discovery: Biology, Function, and 1424 Translation. _Am. J. Hum. Genet._ **101**, 5–22 (2017). 1425 17. De Walsche, A. _et al._ metaGE: Investigating genotype x environment interactions 1426 through GWAS meta-analysis. _PLoS Genet._ **21**, 1–26 (2025). 1427 18. Wu, S. _et al._ Mapping the Arabidopsis Metabolic Landscape by Untargeted 1428 Metabolomics at Different Environmental Conditions. _Mol. Plant_ **11**, 118–134 1429 (2018). 1430 19. Luzarowska, U. _et al._ Hello darkness, my old friend: 3-KETOACYL-COENZYME A 1431 SYNTHASE4 is a branch point in the regulation of triacylglycerol synthesis in 1432 Arabidopsis thaliana. _Plant Cell_ **35**, 1984–2005 (2023). 1433 20. Zhang, F. _et al._ Genomic basis underlying the metabolome-mediated drought 1434 adaptation of maize. _Genome Biol._ **22**, 1–26 (2021). 1435 21. Yonekura-Sakakibara, K., Higashi, Y. & Nakabayashi, R. The Origin and 1436 Evolution of Plant Flavonoid Metabolism. _Front. Plant Sci._ **10**, 1–16 (2019). 1437 22. Tsugawa, H., Rai, A., Saito, K. & Nakabayashi, R. Metabolomics and 1438 complementary techniques to investigate the plant phytochemical cosmos. _Nat._ 1439 _Prod. Rep._ **38**, 1729–1759 (2021). 1440 23. Horai, H. _et al._ MassBank: A public repository for sharing mass spectral data for 1441 life sciences. _J. Mass Spectrom._ **45**, 703–714 (2010). 1442 24. Wang, M. _et al._ Sharing and community curation of mass spectrometry data with 1443 Global Natural Products Social Molecular Networking. _Nat. Biotechnol._ **34**, 828– 1444 837 (2016). 1445 25. Nothias, L.-F. _et al._ Feature-based molecular networking in the GNPS analysis 1446 environment. _Nat. Methods_ **17**, 905–908 (2020). 1447 26. Dührkop, K. _et al._ SIRIUS 4: a rapid tool for turning tandem mass spectra into 1448 metabolite structure information. _Nat. Methods_ **16**, 299–302 (2019). 1449 27. van der Hooft, J. J. J., Wandy, J., Barrett, M. P., Burgess, K. E. V & Rogers, S. 60 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1450 Topic modeling for untargeted substructure exploration in metabolomics. _Proc._ 1451 _Natl. Acad. Sci._ **113**, 1–6 (2016). 1452 28. van Tetering, L. _et al._ A spectroscopic test suggests that fragment ion structure 1453 annotations in MS/MS libraries are frequently incorrect. _Commun. Chem._ **7**, 1–11 1454 (2024). 1455 29. Morreel, K. _et al._ Systematic structural characterization of metabolites in 1456 Arabidopsis via candidate substrate-product pair networks. _Plant Cell_ **26**, 929– 1457 945 (2014). 1458 30. Hartmann, A. C. _et al._ Meta-mass shift chemical profiling of metabolomes from 1459 coral reefs. _Proc. Natl. Acad. Sci._ **114**, 11685–11690 (2017). 1460 31. Naake, T. & Fernie, A. R. MetNet: Metabolite Network Prediction from High 1461 Resolution Mass Spectrometry Data in R Aiding Metabolite Annotation. _Anal._ 1462 _Chem._ **91**, 1768–1772 (2019). 1463 32. Chen, L. _et al._ Metabolite discovery through global annotation of untargeted 1464 metabolomics data. _Nat. Methods_ **18**, (2021). 1465 33. Shen, X. _et al._ Metabolic reaction network-based recursive metabolite annotation 1466 for untargeted metabolomics. _Nat. Commun._ **10**, (2019). 1467 34. Zhu, A. _et al._ Chemical-tag-based semi-annotated metabolomics facilitates gene 1468 identification and specialized metabolic pathway elucidation in wheat. _Plant Cell_ 1469 **36**, 540–558 (2024). 1470 35. Baron, C., Cherkaoui, S., Therrien-, S., Des, C. & Julie, G. Gene-metabolite 1471 annotation with shortest reactional distance enhances metabolite genome- wide 1472 association studies results. _iScience_ **26**, 1–16 (2023). 1473 36. Wen, W. _et al._ Combining quantitative genetics approaches with regulatory 1474 network analysis to dissect the complex metabolism of the maize kernel. _Plant_ 1475 _Physiol._ **170**, 136–146 (2016). 1476 37. Wei, W. _et al._ QTL analysis of important agronomic traits and metabolites in foxtail 61 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1477 millet (Setaria italica) by RIL population and widely targeted metabolome. _Front._ 1478 _Plant Sci._ **13**, 1–19 (2023). 1479 38. Brouckaert, M. _et al._ QT – GWAS : A novel method for unveiling biosynthetic loci 1480 affecting qualitative metabolic traits. _Mol. Plant_ **16**, 1212–1227 (2023). 1481 39. Aharoni, A., Goodacre, R. & Fernie, A. R. Plant and microbial sciences as key 1482 drivers in the development of metabolomics research. _Proc. Natl. Acad. Sci. U. S._ 1483 _A._ **120**, 1–11 (2023). 1484 40. Bowne, J. B. _et al._ Drought responses of leaf tissues from wheat cultivars of 1485 differing drought tolerance at the metabolite level. _Mol. Plant_ **5**, 418–29 (2012). 1486 41. Benhassaine-Kesri, G., Aid, F., Demandre, C., Kader, J. C. & Mazliak, P. Drought 1487 stress affects chloroplast lipid metabolism in rape (Brassica napus) leaves. 1488 _Physiol. Plant._ **115**, 221–227 (2002). 1489 42. Gigon, A., Matos, A. R., Laffray, D., Zuily-Fodil, Y. & Pham-Thi, A. T. Effect of 1490 drought stress on lipid metabolism in the leaves of Arabidopsis thaliana (Ecotype 1491 Columbia). _Ann. Bot._ **94**, 345–351 (2004). 1492 43. Egea, I. _et al._ The drought-tolerant Solanum pennellii regulates leaf water loss 1493 and induces genes involved in amino acid and ethylene/jasmonate metabolism 1494 under dehydration. _Sci. Rep._ **8**, 1–14 (2018). 1495 44. de Souza, L. P., Naake, T., Tohge, T. & Fernie, A. R. From chromatogram to 1496 analyte to metabolite. How to pick horses for courses from the massive web 1497 resources for mass spectral plant metabolomics. _Gigascience_ **6**, 1–20 (2017). 1498 45. Alseekh, S. & Fernie, A. R. Expanding our coverage: Strategies to detect a 1499 greater range of metabolites. _Curr. Opin. Plant Biol._ **73**, 102335 (2023). 1500 46. Aron, A. T. _et al._ Reproducible molecular networking of untargeted mass 1501 spectrometry data using GNPS. _Nat. Protoc._ **15**, 1954–1991 (2020). 1502 47. Shahaf, N. _et al._ The WEIZMASS spectral library for high-confidence metabolite 1503 identification. _Nat. Commun._ **7**, (2016). 62 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1504 48. Alseekh, S. _et al._ Mass spectrometry-based metabolomics: a guide for annotation, 1505 quantification and best reporting practices. _Nat. Methods_ **18**, 747–756 (2021). 1506 49. Alseekh, S. _et al._ Identification and mode of inheritance of quantitative trait loci for 1507 secondary metabolite abundance in tomato. _Plant Cell_ **27**, 485–512 (2015). 1508 50. Fraser, C. M. _et al._ Related arabidopsis serine carboxypeptidase-like 1509 sinapoylglucose acyltransferases display distinct but overlapping substrate 1510 specificities. _Plant Physiol._ **144**, 1986–1999 (2007). 1511 51. Khersonsky, O. & Tawfik, D. S. Enzyme promiscuity: A mechanistic and 1512 evolutionary perspective. _Annu. Rev. Biochem._ **79**, 471–505 (2010). 1513 52. Wang, S., Alseekh, S., Fernie, A. R. & Luo, J. The Structure and Function of 1514 Major Plant Metabolite Modifications. _Mol. Plant_ **12**, 899–919 (2019). 1515 53. Kim, S. _et al._ PubChem Substance and Compound databases. _Nucleic Acids Res._ 1516 **4**, D1202-13 (2016). 1517 54. Wang, F. _et al._ CFM-ID 4.0: More Accurate ESI-MS/MS Spectral Prediction and 1518 Compound Identification. _Anal. Chem._ **93**, 11692–11700 (2021). 1519 55. Mugford, S. T. & Osbourn, A. Evolution of serine carboxypeptidase-like 1520 acyltransferases in the Monocots. _Plant Signal. Behav._ **5**, 193–195 (2010). 1521 56. Stehle, F., Brandt, W., Stubbs, M. T., Milkowski, C. & Strack, D. 1522 Sinapoyltransferases in the light of molecular evolution. _Phytochemistry_ **70**, 1652– 1523 1662 (2009). 1524 57. Fraser, C. M., Rider, L. W. & Chapple, C. An expression and bioinformatics 1525 analysis of the Arabidopsis serine carboxypeptidase-like gene family. _Plant_ 1526 _Physiol._ **138**, 1136–1148 (2005). 1527 58. Aharoni, A. _et al._ Terpenoid Metabolism in Wild-Type and Transgenic Arabidopsis 1528 Plants. _Plant Cell_ **15**, 2866–2884 (2003). 1529 59. Gachon, C., Baltz, R. & Saindrenan, P. Over-expression of a scopoletin 1530 glucosyltransferase in Nicotiana tabacum leads to precocious lesion formation 63 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1531 during the hypersensitive response to tobacco mosaic virus but does not affect 1532 virus resistance. _Plant Mol. Biol._ **54**, 137–146 (2004). 1533 60. Li, A. X. & Steffens, J. C. An acyltransferase catalyzing the formation of 1534 diacylglucose is a serine carboxypeptidase-like protein. _Proc. Natl. Acad. Sci. U._ 1535 _S. A._ **97**, 6902–6907 (2000). 1536 61. Kind, T. & Fiehn, O. Seven Golden Rules for heuristic filtering of molecular 1537 formulas obtained by accurate mass spectrometry. **20**, 1–20 (2007). 1538 62. Foerster, H. _et al._ SolCyc: A database hub at the Sol Genomics Network (SGN) 1539 for the manual curation of metabolic networks in Solanum and Nicotiana specific 1540 databases. _Database_ **2018**, 1–13 (2018). 1541 63. Mansell, D. J. _et al._ Biocatalytic asymmetric alkene reduction: Crystal structure 1542 and characterization of a double bond reductase from Nicotiana tabacum. _ACS_ 1543 _Catal._ **3**, 370–379 (2013). 1544 64. Yin, B. _et al._ Dynamic metabolite QTL analyses provide novel biochemical 1545 insights into kernel development and nutritional quality improvement in common 1546 wheat. _Plant Commun._ **5**, 100792 (2024). 1547 65. Xiang, Y. _et al._ Functional divergence of tandem duplicate SCPL acyltransferase 1548 genes reveals convergent evolution of chicoric acid biosynthesis. _Sci. Adv._ **11**, 1549 eady8282 (2025). 1550 66. Ciarkowska, A., Ostrowski, M. & Jakubowska, A. A serine carboxypeptidase-like 1551 acyltransferase catalyzes synthesis of indole-3-acetic (IAA) ester conjugate in rice 1552 (Oryza sativa). _Plant Physiol. Biochem._ **125**, 126–135 (2018). 1553 67. Mugford, S. T. _et al._ A serine carboxypeptidase-like acyltransferase is required for 1554 synthesis of antimicrobial compounds and disease resistance in oats. _Plant Cell_ 1555 **21**, 2473–2484 (2009). 1556 68. Ahmad, M. Z. _et al._ Genome-Wide Analysis of Serine Carboxypeptidase-Like 1557 Acyltransferase Gene Family for Evolution and Characterization of Enzymes 1558 Involved in the Biosynthesis of Galloylated Catechins in the Tea Plant (Camellia 64 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1559 sinensis). _Front. Plant Sci._ **11**, 1–21 (2020). 1560 69. Chen, J., Li, W. & Yan-xia Jia. The Serine Carboxypeptidase-Like Gene SCPL41 1561 Negatively Regulates Membrane Lipid Metabolism in Arabidopsis thaliana. _Plants_ 1562 **9**, 696 (2020). 1563 70. Zhao, Y. _et al._ Flavan-3-ol Galloylation-Related Functional Gene Cluster and the 1564 Functional Diversification of SCPL Paralogs in Camellia sp. _J. Agric. Food Chem._ 1565 **71**, 488–498 (2023). 1566 71. Chen, X. _et al._ Deactivating mutations in the catalytic site of a companion serine 1567 carboxypeptidase-like acyltransferase enhance catechin galloylation in Camellia 1568 plants. _Hortic. Res._ **12**, (2025). 1569 72. Lehfeldt, C. _et al._ Cloning of the SNG1 gene of arabidopsis reveals a role for a 1570 serine carboxypeptidase-like protein as an acyltransferase in secondary 1571 metabolism. _Plant Cell_ **12**, 1295–1306 (2000). 1572 73. Shirley, A. M., McMichael, C. M. & Chapple, C. The sng2 mutant of arabidopsis is 1573 defective in the gene encoding the serine carboxypeptidase-like protein 1574 sinapoylglucose:choline sinapoyltransferase. _Plant J._ **28**, 83–94 (2001). 1575 74. Tohge, T. _et al._ Characterization of a recently evolved flavonol 1576 phenylacyltransferase gene provides signatures of natural light selection in 1577 Brassicaceae. _Nat. Commun._ **7**, (2016). 1578 75. Lee, S. _et al._ Benzoylation and sinapoylation of glucosinolate R-groups in 1579 Arabidopsis. _Plant J._ **72**, 411–422 (2012). 1580 76. Stehle, F., Brandt, W., Schmidt, J., Milkowski, C. & Strack, D. Activities of 1581 Arabidopsis sinapoylglucose:malate sinapoyltransferase shed light on functional 1582 diversification of serine carboxypeptidase-like acyltransferases. _Phytochemistry_ 1583 **69**, 1826–1831 (2008). 1584 77. Fu, R. _et al._ Versatility in acyltransferase activity completes chicoric acid 1585 biosynthesis in purple coneflower. _Nat. Commun._ **12**, 1562 (2021). 65 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1586 78. Qiu, F. _et al._ Functional genomics analysis reveals two novel genes required for 1587 littorine biosynthesis. _New Phytol._ **225**, 1906–1914 (2020). 1588 79. Huang, Y. _et al._ SCPL acyltransferases catalyze the metabolism of chlorogenic 1589 acid during purple coneflower seed germination. _New Phytol._ **243**, 229–239 1590 (2024). 1591 80. Wang, Y. _et al._ Genome-wide analysis of serine carboxypeptidase-like protein 1592 (SCPL) family and functional validation of Gh_SCPL42 unchromosome conferring 1593 cotton Verticillium der Verticillium wilt stress in Gossypium hirsutum. _BMC Plant_ 1594 _Biol._ **22**, 1–13 (2022). 1595 81. He, L., Liu, Q. & Han, S. Genome-Wide Analysis of Serine Carboxypeptidase-like 1596 Genes in Soybean and Their Roles in Stress Resistance. _Int. J. Mol. Sci._ **25**, 1597 (2024). 1598 82. Jamil, A. _et al._ A genome-wide approach to the comprehensive analysis of the 1599 serine carboxypeptidase-like protein family in Solanum lysopersicum L. _Genet._ 1600 _Resour. Crop Evol._ **72**, 10185–10211 (2025). 1601 83. Naake, T., D’Auria, J. C., Fernie, A. R. & Scossa, F. Phylogenomic and synteny 1602 analysis of BAHD and SCP/SCPL gene families reveal their evolutionary histories 1603 in plant specialized metabolism. _Philos. Trans. R. Soc. B Biol. Sci._ **379**, (2024). 1604 84. Bontpart, T., Cheynier, V., Ageorges, A. & Terrier, N. BAHD or SCPL 1605 acyltransferase? What a dilemma for acylation in the world of plant phenolic 1606 compounds. _New Phytol._ **208**, 695–707 (2015). 1607 85. Girke, C., Daumann, M., Niopek-Witz, S., Möhlmann, T. & Möhlmann, O. 1608 Nucleobase and nucleoside transport and integration into plant metabolism. _Front._ 1609 _Plant Sci._ **5**, 443 (2014). 1610 86. Kato, K., Shitan, N., Shoji, T. & Hashimoto, T. Tobacco NUP1 transports both 1611 tobacco alkaloids and vitamin B6. _Phytochemistry_ **113**, 33–40 (2015). 1612 87. Kakegawa, H. _et al._ Uptake of adenine by purine permeases of Coffea canephora. 1613 _Biosci. Biotechnol. Biochem._ **83**, 1300–1305 (2019). 66 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1614 88. Gillissen, B. _et al._ A new family of high-affinity transporters for adenine, cytosine, 1615 and purine derivatives in Arabidopsis. _Plant Cell_ **12**, 291–300 (2000). 1616 89. Kato, K., Shoji, T. & Hashimoto, A. Tobacco nicotine uptake permease regulates 1617 the expression of a key transcription factor gene in the nicotine biosynthesis 1618 pathway. _Plant Physiol._ **166**, 2195–2204 (2014). 1619 90. Hildreth, S. B. _et al._ Tobacco nicotine uptake permease (NUP1) affects alkaloid 1620 metabolism. _Proc. Natl. Acad. Sci._ **108**, 18179–18184 (2011). 1621 91. Dastmalchi, M. _et al._ Purine Permease-Type Benzylisoquinoline Alkaloid 1622 Transporters in Opium Poppy. _Plant Physiol._ **181**, 916–933 (2019). 1623 92. Hu, Y. & Shani, E. Tansley insight Cytokinin activity – transport and homeostasis 1624 at the whole plant, cell, and subcellular levels. _New Phytol._ **239**, 1603–1608 1625 (2023). 1626 93. Qi, Z. & Xiong, L. Characterization of a Purine Permease Family Gene OsPUP7 1627 Involved in Growth and Development Control in Rice. _J. Integr. Plant Biol._ **55**, 1628 1119–1135 (2013). 1629 94. Zürcher, E., Liu, J., Donato, M. di, Geisler, M. & Müller, B. Plant development 1630 regulated by cytokinin sinks. _Science (80-. )._ **353**, 1027–1030 (2016). 1631 95. Demurtas, O. C., Nicolia, A. & Diretto, G. Terpenoid Transport in Plants: How Far 1632 from the Final Picture? _Plants_ **12**, 634 (2023). 1633 96. Wheaton, T. A. & Stewart, I. Feruloylputrescine : Isolation and Identification from 1634 Citrus Leaves and Fruit. _Nature_ **206**, 620–621 (1965). 1635 97. Martin-Tanguy, J. The occurrence and possible function of hydroxycinnamoyl acid 1636 amides in plants. _Plant Growth Regul._ **3**, 381–399 (1985). 1637 98. Roumani, M., Besseau, S., Hehn, A. & Larbat, R. Functional characterization of a 1638 small gene family coding for putrescine hydroxycinnamoyltransferases, involved 1639 in phenolamide accumulation, in tomato. _Phytochemistry_ **229**, (2025). 1640 99. Roumani, M. _et al._ Transcriptomics and Metabolomics Analyses Reveal High 67 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1641 Induction of the Phenolamide Pathway in Tomato Plants Attacked by the 1642 Leafminer Tuta absoluta. _Metabolites_ **12**, (2022). 1643 100. Mintz-Oron, S. _et al._ Gene Expression and Metabolism in Tomato Fruit Surface 1644 Tissues. _Plant Physiol._ **147**, 823–851 (2008). 1645 101. Adato, A. _et al._ Fruit-Surface Flavonoid Accumulation in Tomato Is Controlled by a 1646 SlMYB12-Regulated Transcriptional Network. _PLoS Genet._ **5**, e1000777 (2009). 1647 102. Breitel, D. A. _et al._ AUXIN RESPONSE FACTOR 2 Intersects Hormonal Signals in 1648 the Regulation of Tomato Fruit Ripening. _PLoS Genet._ **12**, e1005903 (2016). 1649 103. Wang, W. _et al._ Temperature-relied defense: CsMYBs-CsHT9 modules activate 1650 N-feruloylputrescine biosynthesis against anthracnose in tea (Camellia sinensis). 1651 _Plant J._ **125**, 1–16 (2026). 1652 104. Nguyen, V. N. T. _et al._ OsABCG37 mediates root development and the stress 1653 response by modulating feruloylputrescine levels. _Plant Physiol._ **199**, 1–14 1654 (2025). 1655 105. Hasanuzzaman, M. _et al._ Regulation of ROS metabolism in plants under 1656 environmental stress: A review of recent experimental evidence. _Int. J. Mol. Sci._ 1657 **21**, 1–44 (2020). 1658 106. Sonawane, P. D., Jozwiak, A., Panda, S. & Aharoni, A. ‘Hijacking’ core 1659 metabolism: a new panache for the evolution of steroidal glycoalkaloids structural 1660 diversity. _Curr. Opin. Plant Biol._ **55**, 118–128 (2020). 1661 107. Nogia, P. & Pati, P. K. Plant Secondary Metabolite Transporters: Diversity, 1662 Functionality, and Their Modulation. _Front. Plant Sci._ **12**, 1–20 (2021). 1663 108. Mano, J. _et al._ Protection against photooxidative injury of tobacco leaves by 2 1664 alkenal reductase. Detoxication of lipid peroxide-derived reactive carbonyls. _Plant_ 1665 _Physiol._ **139**, 1773–1783 (2005). 1666 109. Youn, B. _et al._ Mechanistic and structural studies of apoform, binary, and ternary 1667 complexes of the Arabidopsis alkenal double bond reductase At5g16970. _J. Biol._ 68 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1668 _Chem._ **281**, 40076–40088 (2006). 1669 110. Zhu, H., Ding, Y., Sun, X. & Wu, Y. The α,β-Double Bond Reductases (DBRs): An 1670 Overview. _J. Evol. Biochem. Physiol._ **61**, 820–834 (2025). 1671 111. Mata-Pérez, C. _et al._ Endogenous Biosynthesis of S-Nitrosoglutathione From 1672 Nitro-Fatty Acids in Plants. _Front. Plant Sci._ **11**, (2020). 1673 112. Xiong, Z. _et al._ Functional characterization of a Colchicum autumnale L. double 1674 bond reductase (CaDBR1) in colchicine biosynthesis. _Planta_ **256**, 1–12 (2022). 1675 113. Klein, D., Fink, B., Arold, B., Eisenreich, W. & Schwab, W. Functional 1676 characterization of enone oxidoreductases from strawberry and tomato fruit. _J._ 1677 _Agric. Food Chem._ **55**, 6705–6711 (2007). 1678 114. Ringer, K. L., McConkey, M. E., Davis, E. M., Rushing, G. W. & Croteau, R. 1679 Monoterpene double-bond reductases of the (−)-menthol biosynthetic pathway: 1680 isolation and characterization of cDNAs encoding (−)-isopiperitenone reductase 1681 and (+)-pulegone reductase of peppermint. _Arch. Biochem. Biophys._ **418**, 80–92 1682 (2003). 1683 115. Mączka, W., Sołtysik, D., Wińska, K., Grabarczyk, M. & Szumny, A. Plant 1684 Mediated Biotransformations of S(+)- and R(–)-Carvones. _Appl. Sci._ **8**, 2605 1685 (2018). 1686 116. Yang, M.-T. _et al._ Authentication, phytochemical characterization and anti 1687 bacterial activity of two Artemisia species. _Food Chem._ **333**, 127458 (2020). 1688 117. ChongRen, S. _et al._ Effects of different plant inhibitors on the growth and toxin 1689 production of Aspergillus flavus. in _3rd International Conference on Food Quality_ 1690 _and Safety._ 2892–2897 (2017). 1691 118. Xiao, D. _et al._ Exogenous Application of a Plant Elicitor Induces Volatile Emission 1692 in Wheat and Enhances the Attraction of an Aphid Parasitoid Aphidius gifuensis. 1693 _Plants_ **11**, 3496 (2022). 1694 119. Cui, S. _et al._ Development of Fast E-nose System for Early-Stage Diagnosis of 69 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1695 Aphid-Stressed Tomato Plants. _Sensors_ **19**, 3480 (2019). 1696 120. Zhang, H. _et al._ Beta-ionone is a functional plant volatile that attracts the parasitic 1697 wasp, Microplitis pallidipes. _BioControl_ **67**, 1–13 (2022). 1698 121. Wang, Y., Xu, J. & Liu, A. Identification of the carotenoid cleavage dioxygenase 1699 genes and functional analysis reveal DoCCD1 is potentially involved in beta 1700 ionone formation in Dendrobium officinale. _Front. Plant Sci._ **13**, 967819 (2022). 1701 122. Simkin, A. J., Schwartz, S. H., Auldridge, M., Taylor, M. G. & Klee, H. J. The 1702 tomato carotenoid cleavage dioxygenase 1 genes contribute to the formation of 1703 the flavor volatiles β-ionone, pseudoionone, and geranylacetone. _Plant J._ **40**, 1704 882–892 (2004). 1705 123. Manzo, D. Integrated–omics approaches to explore tomato interaction with the 1706 leafminer Tuta absoluta. (UNIVERSITY OF NAPLES ‘FEDERICO II’’’, 2016). 1707 124. He, Q. _et al._ Functional diversification of WRINKLED3 integrates fatty acid 1708 metabolism with insecticidal acylsugar production in Solanaceae species. _Nat._ 1709 _Commun._ **17**, (2026). 1710 125. Lapidot-Cohen, T. Genetic regulation of agricultural traits and metabolites in 1711 tomato in control and water deficit environments. (Ben Gurion University of the 1712 Negev, 2021). 1713 126. Team, R. C. R: A language and environment for statistical computing. vol. 201 at 1714 https://www.r-project.org/ (2021). 1715 127. Tzfadia, O. _et al._ The ‘TranSeq’ 3′-end sequencing method for high-throughput 1716 transcriptomics and gene space refinement in plant genomes. _Plant J._ **96**, 223– 1717 232 (2018). 1718 128. Anders, S. & Huber, W. Differential expression analysis for sequence count data. 1719 _Nat. Preced._ (2010) doi:https://doi.org/10.1038/npre.2010.4282.1. 1720 129. Stacklies, W., Redestig, H., Scholz, M., Walther, D. & Selbig, J. pcaMethods-a 1721 bioconductor package providing PCA methods for incomplete data. _Bioinformatics_ 70 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1722 **23**, 1164–1167 (2007). 1723 130. Giavalisco, P., Köhl, K., Hummel, J., Seiwert, B. & Willmitzer, L. 13C isotope 1724 labeled metabolomes allowing for improved compound annotation and relative 1725 quantification in liquid chromatography-mass spectrometry-based metabolomic 1726 research. _Anal. Chem._ **81**, 6546–6551 (2009). 1727 131. Giavalisco, P. _et al._ Elemental formula annotation of polar and lipophilic 1728 metabolites using 13C, 15N and 34S isotope labelling, in combination with high 1729 resolution mass spectrometry. _Plant J._ **68**, 364–376 (2011). 1730 132. Pence, H. E. & Williams, A. ChemSpider: An Online Chemical Information 1731 Resource. _J. Chem. Educ._ **87**, 1123–1124 (2010). 1732 133. Schwahn, K., de Souza, L. P., Fernie, A. R. & Tohge, T. Metabolomics-assisted 1733 refinement of the pathways of steroidal glycoalkaloid biosynthesis in the tomato 1734 clade. _J. Integr. Plant Biol._ **56**, 864–875 (2014). 1735 134. Iijima, Y. _et al._ Metabolite annotations based on the integration of mass spectral 1736 information. _Plant J._ **54**, 949–962 (2008). 1737 135. Caprioli, G., Cahill, M., Logrippo, S. & James, K. Elucidation of the mass 1738 fragmentation pathways of tomatidine and β1-hydroxytomatine using orbitrap 1739 mass spectrometry. _Nat. Prod. Commun._ **10**, 575–576 (2015). 1740 136. Rodríguez-Pérez, C., Quirantes-Piné, R., Fernández-Gutiérrez, A. & Segura 1741 Carretero, A. Comparative characterization of phenolic and other polar 1742 compounds in Spanish melon cultivars by using high-performance liquid 1743 chromatography coupled to electrospray ionization quadrupole-time of flight mass 1744 spectrometry. _Food Res. Int._ **54**, 1519–1527 (2013). 1745 137. Wang, J. _et al._ A Platform to Identify Endogenous Metabolites Using a Novel High 1746 Performance Orbitrap and the mzCloud Library. vol. 25 2009 at (2009). 1747 138. Sawada, Y. _et al._ RIKEN tandem mass spectral database (ReSpect) for 1748 phytochemicals: A plant-specific MS/MS-based data resource and database. 1749 _Phytochemistry_ **82**, 38–45 (2012). 71 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1750 139. Kessner, D., Chambers, M., Burke, R., Agus, D. & Mallick, P. ProteoWizard: open 1751 source software for rapid proteomics tools development. _Bioinformatics_ **24**, 2534– 1752 2536 (2008). 1753 140. Tsugawa, H. _et al._ MS-DIAL: Data-independent MS/MS deconvolution for 1754 comprehensive metabolome analysis. _Nat. Methods_ **12**, 523–526 (2015). 1755 141. Shannon, P. _et al._ Cytoscape: A Software Environment for Integrated Models. 1756 _Genome Res._ **13**, 2498–2504 (2003). 1757 142. Allen, F., Pon, A., Wilson, M., Greiner, R. & Wishart, D. CFM-ID: a web server for 1758 annotation, spectrum prediction and metabolite identification from tandem mass 1759 spectra. _Nucleic Acids Res._ **42**, 94–99 (2014). 1760 143. Dodder, N. & Mullen, K. Package ‘OrgMassSpecR’. at 1761 http://orgmassspec.github.io/ (2017). 1762 144. Lipka, A. E. _et al._ GAPIT: Genome association and prediction integrated tool. 1763 _Bioinformatics_ **28**, 2397–2399 (2012). 1764 145. Tang, Y. _et al._ GAPIT Version 2: An Enhanced Integrated Tool for Genomic 1765 Association and Prediction. _Plant Genome_ **9**, plantgenome2015.11.0120 (2016). 1766 146. Zhang, Z. _et al._ Mixed linear model approach adapted for genome-wide 1767 association studies. _Nat. Genet._ **42**, 355–360 (2010). 1768 147. Rainer, J. _et al._ A Modular and Expandable Ecosystem for Metabolomics Data 1769 Annotation in R. _Metabolites_ **12**, (2022). 1770 148. Depke, T., Franke, R. & Brönstrup, M. CluMSID: An R package for similarity 1771 based clustering of tandem mass spectra to aid feature annotation in 1772 metabolomics. _Bioinformatics_ **35**, 3196–3198 (2019). 1773 149. McDonald, J. . G–test of goodness-of-fit. in _Handbook of Biological Statistics_ 53– 1774 58 (Sparky House Publishing, 2014). 1775 150. Wei, T. _et al._ Package ‘corrplot’. _Statistician_ **56**, e24 (2017). 1776 151. Shin, J.-H., Blay, S., McNeney, B. & Graham, J. LDheatmap: an R function for 72 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1777 graphical display of pairwise linkage disequilibria between single nucleotide 1778 polymorphisms. _J. Stat. Softw._ **16**, 1–9 (2006). 1779 152. Moscatiello, R., Baldan, B. & Navazio, L. Plant Cell Suspension Cultures BT 1780 Plant Mineral Nutrients: Methods and Protocols. in (ed. Maathuis, F. J. M.) 77–93 1781 (Humana Press, 2013). doi:10.1007/978-1-62703-152-3_5. 1782 153. McCormick, S. Transformation of tomato with Agrobacterium tumefaciens BT 1783 Plant Tissue Culture Manual: Supplement 7. in (ed. Lindsey, K.) 311–319 1784 (Springer Netherlands, 1997). doi:10.1007/978-94-009-0103-2_17. 1785 154. Ouyang, B. _et al._ Transformation of tomatoes with osmotin and chitinase genes 1786 and their resistance to Fusarium wilt. _J. Hortic. Sci. Biotechnol._ **80**, 517–522 1787 (2005). 1788 155. Livak, K. J. & Schmittgen, T. D. Analysis of Relative Gene Expression Data Using 1789 Real-Time Quantitative PCR and the 2−ΔΔCT Method. _Methods_ **25**, 402–408 1790 (2001). 1791 1792 1793 73 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1794 **Figure legends** 1795 1796 Figure 1. 1797 Figure 1: **Experiment overview and stress-induced variation** 1798 A. Overview of experiment setup, datasets produced and the main analysis and results 1799 to which they led. 1800 B. Leaf water content, as the ratio of leaf weight lost after full oven drying to the fresh 1801 leaf weight, in leaves collected three weeks after treatment onset from control (light 1802 blue) and water deficit (WD; orange) treated plants. N control =155, N WD =54. 1803 C. The fresh vegetative weight of whole plants measured after last fruit harvest. 1804 N control =656, N WD =589. 1805 D. Global effect of WD on leaf gene expression. PCA based on relative gene expression 1806 of leaves collected from control (blue) and water deficit (dark orange) treated plants. 1807 E. Distribution of variation (standard deviation) in the expression of each gene across 1808 the population in either growth condition. 1809 F. Gene expression in response to WD. Volcano plot displaying the significance of a t 1810 test comparing transcript levels between water deficit and control conditions to the fold 1811 change (after log10 transformation) of the same comparison. Transcripts with significant 1812 increase (p-value <α=0.05, Bonferroni correction) and fold change exceeding double 1813 expression in water deficit (dark orange) or control (blue) are colored. 1814 G. Global effect of WD on metabolic profile. PCA based on metabolic profiles of leaves 1815 collected from control (light blue) and water deficit (orange) treated plants with 1816 metabolite features detected in both positive and negative ionization modes. 1817 H. Feature changes in response to WD. Volcano plot displaying the significance of a 1818 paired t-test comparing feature levels of each accession between WD and control 1819 conditions to the median fold change (after log2 transformation). Features with 1820 significant increase (p-value <α=0.05, Bonferroni correction) and fold change exceeding 1821 two-fold expression in water deficit (orange) or control (light blue) are colored. 1822 I. The number of features detected exclusively in control (blue circle), water deficit 1823 (orange circle) or in both (overlap). 74 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1824 Figure 2. 1825 Figure 2: **mQTL overview and scenarios** 1826 A. Each detected mQTL is represented as a column placed according to physical 1827 location on the tomato genome (x axis). The width of the column ranges from the first to 1828 the last significant SNP in the mQTL. The height of the column represents the number 1829 of metabolite features in the mQTL feature group, scale noted on the left y axis. The dot 1830 accompanying each mQTL column represents the significance of the highest SNP 1831 feature association in the mQTL, scale shown on the right y axis. 1832 Results of mapping in control conditions are in the top panel in blue, and from water 1833 deficit in mirror image in the bottom panel in orange. The mQTL and genes studied in 1834 the following sections are marked by arrows. 1835 B. Graphical representation of the scenarios discussed in the text for mQTLs with 1836 large mQTL feature groups. 1837 Structures and colors represent arbitrary molecular substructures. Inset plots are 1838 simplified commonality spectra expected from various scenarios. In scenario IV, the 1839 pleiotropic effects are not expected to manifest in a typical commonality spectrum, so it 1840 is not drawn. 1841 Figure 3. 1842 Figure 3: **Fragment commonality of mQTL feature groups** 1843 A and B. The most common fragments found in all ddMS² spectra obtained from 1844 tomato leaf metabolite features detected in positive (A) or in negative (B) ionization 1845 mode. The m/z value of each fragment is indicated above the corresponding bar, whose 1846 height indicates the frequency of ddMS² spectra which contain the fragment out of all 1847 detected ddMS² spectra. Due to the excess of fragments some numbers overlap, the 1848 ranked list of fragments and their frequency can be found in Supp. Table 5. Hexose 1849 (gray hexagon) derived fragments are marked with gray underline. 1850 C and D. The significance level of mQTL fragment commonality indicating how 1851 different the commonality of an mQTL feature group is compared to the commonality of 1852 all ddMS² spectra. Shown for features detected in positive (C) or in negative (D) 1853 ionization mode. The pink line denotes the significance threshold derived from the 5 [th] 75 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1854 percentile of random ddMS² spectra sampling with corresponding number of ddMS² 1855 spectra. Dashed gray lines note the suggested minimum number of spectra for a 1856 confident analysis. 1857 E and F. The frequency of the most common fragment of each mQTL feature group 1858 (excluding the overall most common fragments from panels A-B), for features detected 1859 in positive (E) ionization mode or in negative (F). The light blue line denotes the 1860 outstanding threshold, the mean top frequency of random ddMS² spectra sampling with 1861 corresponding number of ddMS² spectra. Dashed gray lines note the suggested 1862 minimum number of spectra for a confident analysis. 1863 Colored arrows indicate the mQTLs studied in the following sections, in accordance 1864 with Fig. 2. 1865 1866 Figure 4. 1867 Figure 4: **Transformation commonality for identifying mQTLs' underlying reaction** 1868 A and B. The most common transformations (m/z differences) found for pairs of all 1869 metabolite features detected in positive (A) or in negative (B) ionization mode. The m/z 1870 value of each transformation is indicated above the corresponding bar, whose height 1871 indicates the number of occurrences out of all non-redundant feature pairs. Due to the 1872 excess of transformations some numbers overlap, the full ranked list of frequent 1873 transformations can be found in Supp. Table 8. When available, putative annotations 1874 are noted in purple, according to Supp. Table 7. 1875 C and D. The significance level of mQTL transformation commonality indicating how 1876 different feature pairs of an mQTL feature group are compared to pairs of all features, 1877 for features detected in positive (C) or in negative (D) ionization mode. 1878 The pink line denotes the significance threshold derived from the 5 [th] percentile of 1879 random sampling with corresponding number of feature m/z values. Dashed gray lines 1880 note the minimum number of m/z values suggested for a confident analysis. 1881 E and F. The frequency of the most common transformation of each mQTL feature 1882 group (excluding the overall most common transformations from panels A-B), for 1883 features detected in positive (E) or in negative (F) ionization mode. The blue line 1884 denotes the outstanding threshold, the mean top frequency of the most common 76 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1885 transformation in random m/z samplings, with a corresponding number of feature m/z 1886 values. Dashed gray lines note the minimum number of m/z values suggested for a 1887 confident analysis. 1888 Colored arrows indicate the mQTLs studied in the following sections, in accordance 1889 with Fig. 2. 1890 1891 Figure 5. 1892 Figure 5: **Study of chromosome 12 mQTL based on common fragments** 1893 A. A. Manhattan plot of the significant mQTL region in control conditions. Each color 1894 represents the mGWAS significance of one of the metabolite features associated with 1895 the mQTL. The lower panel shows the linkage disequilibrium analysis of all SNPs in the 1896 extended mQTL region. The red triangle indicates the linkage group of the mQTL 1897 significant SNPs. Orange rectangles denote the genomic location of serine 1898 carboxypeptidase-like (SCPL) family genes. 1899 B. Feature-feature correlations for all features of the mQTL feature group in control 1900 conditions. Features are ordered and squares are drawn according to the similarity 1901 calculated by hierarchical clustering analysis (HCA). Color scale represents the R [2] 1902 coefficient for feature pairs. HCA-based group names were defined according to 1903 correlations in control conditions and are noted on the side bar. 1904 C. The effect of the most significantly associated SNP, S12.2059255, alleles on 1905 relative feature abundance of a representative feature from each HCA group. 1906 D and E. Commonality spectra of the chromosome 12 mQTL feature group 1907 highlighting the ddMS [2] fragments detected in positive (D) or in negative (E) ionization 1908 mode. The fragments most common among all features are presented in gray. In E, 1909 fragments in gray are hexose-derived. Inset in E: molecular structure of sinapoyl 1910 glucose. 1911 F. Correlations between expression of genes in the mQTL region (in rows) and 1912 features (in columns) significantly associated with the mQTL in control conditions. Color 1913 scale represents the R [2] coefficient, as in panel B. Yellow triangles point to the selected 1914 candidate gene, Solyc12g088250. 77 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1915 G. Genome wide association (Manhattan plot) of Solyc12g088250 relative 1916 expression in control conditions. 1917 H. Phylogenetic clustering of genes from the SCPL family based on protein 1918 sequence alignment. The amino acids of the conserved family motif are marked by a 1919 red square. Genes expressed in the tomato leaf samples are highlighted in green. 1920 I. Pairwise gene expression correlation of the leaf-expressed genes of the SCPL 1921 family. Color scale represents the R [2] coefficient, as in panel B. 1922 J. Heatmap of the changes in features of the mQTL feature group in plants 1923 overexpressing Solyc12g088250 (golden top bar) compared to wild-type (green top 1924 bar). Features are ordered by HCA groups from panel B. 1925 1926 Figure 6 1927 Figure 6: **Fragment commonality of predicted precursors and products of two** 1928 **serine carboxypeptidase-like genes** 1929 A. A. Volcano plot displaying the significance of a t-test comparing metabolite 1930 feature abundance between plants overexpressing (OE) Solyc04g077650 and wild-type 1931 (WT; Ailsa Craig accession) against the fold change (after log2 transformation) of the 1932 same comparison. Features which contributed most (top 1%) to the separation of OE 1933 and WT in PLS-DA are colored golden (increased in OE) or green (decreased in OE). 1934 B-E. Commonality spectra of fragments of features decreased by OE of 1935 Solyc04g077650 (B, D) or Solyc12g088250 (C, E) in positive (B, C) or negative (D, E) 1936 ionization modes. 1937 F-H. Commonality spectra of fragments of features increased by OE, expected to be 1938 the products of the reaction as denoted by the arrows, of Solyc04g077650 (G) or 1939 Solyc12g088250 (F, H) in positive (F) or negative (G, H) ionization modes. 1940 1941 Figure 7. 1942 Figure 7: **Study of chromosome 2 mQTL and candidate gene** 1943 A-B. Manhattan plot of the significant mQTL region in control (A) and water deficit (B) 1944 conditions. Each color represents the GWAS results of one of the metabolite features 78 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1945 associated with the mQTL. The lower panel in B shows the linkage disequilibrium of all 1946 SNPs in the extended mQTL region. The red triangle indicates the linkage group which 1947 includes the mQTL significant SNPs. Orange rectangles denote the genomic location of 1948 purine permease family genes. 1949 C. Feature-feature correlations for all features of the mQTL feature group in water 1950 deficit conditions. Features are ordered and grouping squares are drawn according to 1951 the similarity calculated by hierarchical clustering analysis (HCA). Color scale 1952 represents the R [2] coefficient for feature pairs. HCA-based group names are noted on 1953 the side bar. 1954 D. The effect of the most significantly associated SNP, S02.1211713, alleles on 1955 relative feature abundance of representative features from the HCA groups in water 1956 deficit. 1957 E. Correlations between expression of genes in mQTL region (in rows) and features 1958 (in columns) with significant association to the mQTL in water deficit conditions. Color 1959 scale represents the R [2] coefficient, when significant. The yellow triangle points to the 1960 selected candidate gene, Solyc02g071100. 1961 F. The effect of the lead SNP S02.1211713 and of the environment on 1962 Solyc02g071100 gene expression. 1963 G. Genome wide association (Manhattan plot) of Solyc02g071100 relative expression 1964 in water deficit conditions. 1965 H. Fragment commonality spectrum of the mQTL feature group, for features detected 1966 in positive ionization mode. 1967 I. Spectra at retention time 4.44 min derived from [15] N and [13] C isotope labeling 1968 experiments. The comparison of spectra of this representative feature shows the 1969 increase in m/z corresponding to the number of labeled atoms, validating the predicted 1970 (C 14 H 20 N 2 O 3 ) molecular formula of feruloylputrescine. 1971 J. Heatmap of the changes in features of the mQTL feature group in plants 1972 overexpressing Solyc02g071100 (golden top bars) compared to wild-type (green top 1973 bars), under fully watered or water deficit conditions at the indicated time points. 1974 Feature side colors indicate the HCA groups from panel C. 1975 79 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.03.763898; this version posted November 3, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.03.763898) (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. 1976 Figure 8. 1977 Figure 8: **Study of chromosome 9 mQTL based on high transformation** 1978 **commonality** 1979 A and B. Manhattan plot of the significant mQTL region in control (A) and water deficit 1980 (B) conditions. Each color represents the mGWAS significance of one of the metabolite 1981 features associated with the mQTL. 1982 C. Feature-feature correlations for all features of the mQTL feature group in water 1983 deficit conditions. Features are ordered and squares are drawn according to the 1984 similarity calculated by hierarchical clustering analysis (HCA). Color scale represents 1985 the R [2] coefficient for feature pairs. HCA-based group names were defined according to 1986 correlations in control conditions and are noted on the side bar. 1987 D. Focus on the genomic location of the single highly associated SNP, within the 1988 context of the gene Solyc09g091700. 1989 1990 E. The effect of the highly associated SNP, S09.2367600, alleles on relative feature 1991 abundance of representative features from each HCA group in control (light blue) and 1992 water deficit (orange). 1993 F and G. Commonality of transformations between features in HCA groups A and B, 1994 which are negatively correlated to each other (panel C), for features detected in positive 1995 (F) or in negative (G) ionization mode. Inset: molecular structures of annotated 1996 metabolites with the typical double-bond of the mQTL feature group in positive (F) or in 1997 negative (G) ionization mode. 1998 H. Heatmap of the changes in features of the mQTL feature group in plants 1999 overexpressing Solyc09g091700 (golden top bar) compared to wild-type (green top 2000 bar). Features are ordered by HCA groups from panel C. 2001 80