Metabolite-driven mechanisms reveal chemical ecology of Lehmann Lovegrass (Eragrostis lehmanniana) invasion in North American semi-arid ecosystems

Introduction

Invasive plant species pose a significant global threat to biodiversity, ecosystem services, and human health by outcompeting native plants for resources, altering habitat structures, and releasing allergenic pollen and toxins¹. Currently, invasion ecology has focused on identifying distinct functional traits of invasive plants explaining their success in the introduced range compared to native plants^2,3. Invasive plants often exhibit superior nutrient acquisition, accelerated growth, enhanced reproductive and dispersal capabilities, and higher phenotypic plasticity compared to co-occurring native plants^4,5,6,7,8. However, functional traits are aggregates of multiple physiological processes, and the success of a certain invasive plant can be related to multiple interacting traits⁹. This complexity in trait interactions can mask the underlying physiological mechanisms driving plant invasion. Although trait-based frameworks are valuable for predicting invaders and assessing habitat resistance¹⁰, current trait-based approaches provide limited guidance for management decisions and poorly predict how control measures will affect co-occurring native plants¹¹.

The well-documented invasive grass Lehmann lovegrass (Eragrostis lehmanniana, hereafter ERLE) exemplifies the limitations of current functional trait approaches in explaining plant invasion dynamics. Introduced in 1932, ERLE has transformed 145,000 ha of semi-desert grasslands in the American Southwest, dominating >90% of grass biomass in invaded plots¹². Traditionally, ERLE’s success has been attributed to its superior functional traits, including efficient nutrient acquisition, grazing resistance, drought adaptation, and fire tolerance^13,14. However, comparative studies with Arizona cottontop (Digitaria californica, DICA) highlight inconsistencies. Despite ERLE’s seemingly advantageous characteristics—such as a nearly half the shoot C:N ratio suggesting superior nutrient competition, DICA demonstrates greater biomass increase under elevated nitrogen conditions¹⁵ and more rapid shoot regrowth after defoliation than ERLE¹⁶. Field studies also suggest that grazing intensity and extreme drought do not confer unique advantages to ERLE over DICA^17,18. While ERLE is known for quick fire recovery¹⁹, no evidence suggests fire negatively impacts DICA²⁰. These contradictions highlight the inadequacy of current functional trait approaches in fully explaining ERLE’s dominant invasion pattern, highlighting the need for mechanistic insights into plant invasion dynamics.

Recent metabolomic studies have revealed key mechanisms driving plant invasions, including allelopathy against native plants, pathogen and herbivore defenses, interactions with soil microbes, alterations in soil chemistry, responses to environmental stress, and metabolic plasticity^21,22,23, highlighting the potential of metabolomics to bridge the gap between invasive plant physiology and their competitive success. However, individual analytical methods often fall short of capturing the complexity of plant metabolomes. Proton nuclear magnetic resonance (¹H NMR), while providing structural, quantitative, and reproducible data on small metabolites such as root exudates²⁴, is limited in capturing the overall plant metabolome due to challenges in detecting low-concentration metabolites and annotating large-mass metabolites. Fourier transform ion cyclotron resonance (FT-ICR-MS) offers superior coverage (200–1200 Da) and resolution for separating similar m/z ratios²⁵ but remains constrained by its ability to provide only molecular formula assignments. Liquid chromatography tandem mass spectrometry (LC MS/MS), widely adopted in untargeted metabolomics for its structural elucidation capabilities within the 50–1000 Da range²⁶, struggles with feature annotation accuracy and quantitative reliability due to ion suppression issues²⁷. These technical limitations of individual metabolomic methods highlight the need for an integrated analytical approach to comprehensively identify and cross-validate plant metabolomes.

To address the current knowledge gaps in mechanistic insights of invasive plant success and limitations in individual metabolomic approaches, we developed a multimodal metabolomics strategy by integrating four analytical techniques ¹H NMR, LC MS/MS, FT-ICR-MS, as well as Matrix-Assisted Laser Desorption/Ionization Mass Spectrometry Imaging (MALDI-MSI), which enables visualization of metabolite spatial distribution within plant roots²⁸. Leveraging this approach, we investigated three key questions about the ERLE invasion system: (1) How do nitrogen acquisition and allocation strategies differ between invasive ERLE and native DICA at the molecular level, particularly given their contradictory responses to nitrogen availability? (2) What role do defense metabolites play in ERLE’s invasion success, considering the lack of evidence for enemy release effects? (3) How do metabolic responses to environmental variation differ between these species, potentially explaining ERLE’s success under varying conditions? We compared root and rhizosphere metabolomes of ERLE and DICA across open and closed canopy conditions of mesquite trees, where open canopy areas represent nutrient-poor, stressful environments, while closed canopy zones create nutrient-rich, less stressful microhabitats²⁹.

Based on current invasion theory and the observed contradictions in trait-based explanations, we hypothesized that ERLE would exhibit distinct root and rhizosphere metabolomic profiles compared to DICA, reflecting its invasion strategies: (1) enrichment in nitrogen-containing metabolites and upregulated nitrogen assimilation, with its rhizosphere showing enhanced microbial activity, reflecting alternative nitrogen acquisition strategies beyond what traditional C:N measurements reveal⁴. (2) ERLE would have fewer defensive metabolites in roots and rhizospheres, suggesting resource allocation strategies favoring growth over defense³⁰ and (3) ERLE would display greater phenotypic plasticity in its root and rhizosphere metabolomes across environmental conditions (open vs. closed canopy), potentially explaining its success despite showing no apparent drought adaptation advantage.

Our findings reveal molecular mechanisms behind ERLE’s invasion success, including efficient nitrogen allocation to shoots, reduced defense investment, and high metabolic plasticity. The findings demonstrate how molecular-level analysis can resolve apparent contradictions in trait-based approaches and suggest that climate-driven changes in resource distribution may be key drivers of ERLE’s invasion success. This study advances our understanding of invasion mechanisms beyond traditional trait-based approaches, offering new insights for predicting and managing plant invasions under environmental change.

Results

Multimodal metabolomics overview

Our complementary multimodal metabolomics approach (¹H NMR, LC MS/MS, MALDI-MSI, and FT-ICR-MS) revealed a complex chemical landscape in the context of a major invasive plant. Each technique captured a distinct range and number of metabolic features, emphasizing their complementarity (Fig. 1a, e; see supplementary note for exact mass range and number). Rhizosphere metabolomes exhibited more features and greater diversity than root metabolomes, reflecting the complexity of these environments.

**Fig. 1: Overview of results from different metabolomic platforms.**

Root metabolome shaped by plant identity, rhizosphere metabolome shaped by plant and canopy

Distinct metabolomic profiles were observed between invasive ERLE and native DICA using non-metric multidimensional scaling (NMDS) and permutational multivariate analysis of variance (PERMANOVA) analyses of root and rhizosphere samples. FT-ICR-MS analysis of root metabolomes revealed a clear separation between ERLE and DICA, with plant species as the primary driver (Fig. 1b; PERMANOVA: R² = 0.272, p = 0.001 for plant; R² = 0.033, p = 0.076 for canopy). This pattern was consistent across all modes of LC MS/MS but not ¹H NMR, where canopy conditions also significantly influenced ordination (Supplementary Fig. 1). This discrepancy may be attributed to ¹H NMR’s bias towards detecting abundant metabolites with a small molecular weight. The separation of DICA and ERLE in the NMDS plot of FT-ICR-MS result was associated with specific metabolite classes³¹. DICA roots were linked with recalcitrant metabolites such as lignin-like phenolic metabolites and condensed hydrocarbons, suggesting defensive strategies. ERLE roots showed higher levels of unsaturated hydrocarbons-, amino sugars-, and protein-like metabolites, indicating a focus on growth and nutrient acquisition.

In the rhizosphere, metabolomic profiles were primarily shaped by canopy conditions (Fig. 1f; PERMANOVA: R² = 0.164, p = 0.001 for canopy, R² = 0.041, p = 0.124 for plant). However, ¹H NMR and LC MS/MS results showed stronger plant influences (Supplementary Fig. 1). This discrepancy may be attributed to FT-ICR-MS capturing a wider range of metabolites beyond those produced by the plant. The separation of open and closed canopy in the NMDS plot of FT-ICR-MS result is only weakly associated with “other” metabolites. ERLE rhizospheres resembled DICA roots with recalcitrant metabolites, such as lignin- and condensed hydrocarbon-like metabolites (Fig. 1f). Conversely, DICA rhizospheres were enriched in proteins-, lipids-, unsaturated hydrocarbons-, and amino sugar-like metabolites (Fig. 1f), reflecting a pattern more akin to ERLE roots. This discrepancy in metabolite classes between root and rhizosphere samples suggests that DICA may invest differently in root exudates compared to ERLE, potentially leading to distinct rhizosphere microbial communities and metabolic processes.

DICA upregulated root nitrogen assimilation, ERLE upregulated nitrogen carriers

To investigate whether Lehmann lovegrass demonstrates superior nitrogen acquisition, we analyzed nitrogen-containing metabolites in roots using FT-ICR-MS, root metabolites involved in organic nitrogen assimilation pathways via LC MS/MS, and potential nitrogen carriers identified by LC MS/MS and MALDI-MS. FT-ICR-MS revealed a higher proportion of nitrogen-containing metabolites in ERLE roots compared to DICA roots (Fig. 2a; Kruskal-Wallis test: χ² = 11.036, p = 0.008 for plant; χ² = 2.027, p = 0.155 for canopy). However, LC MS/MS and ¹H NMR results indicated that most metabolites, particularly those involved in organic nitrogen assimilation pathways³², were upregulated in DICA, especially under closed canopy conditions (Fig. 2b, Supplementary Note 2, 3).

**Fig. 2: Root and rhizosphere metabolites involved in nitrogen use strategies.**

We focused on three major nitrogen carriers—asparagine³³, allantoin³⁴, and glutamine³⁵—to assess nitrogen uptake in roots. ERLE significantly upregulated two of these carriers, particularly showing notably higher levels of asparagine irrespective of canopy conditions (Fig. 2c1; Kruskal-Wallis test: χ² = 20.469, p < 0.0001 for plant; χ² = 0.144, p = 0.704 for canopy). ERLE also exhibited increased glutamine levels under open canopy conditions (Fig. 2c3; Dunn’s test: z-value = 2.77, p.adj = 0.033 for DICA-open vs ERLE-open; Fig. 2c3). In contrast, allantoin levels did not differ significantly between the species but were higher under closed canopy conditions (Fig. 2c2; Kruskal-Wallis test: χ² = 0.361, p = 0.548 for plant; χ² = 13.009, p = 0.0003 for canopy). MALDI-MSI revealed that glutamine and allantoin were uniformly distributed across the entire root section, whereas asparagine in ERLE roots was localized predominantly around the xylem, suggesting its role in nitrogen transport within the plant (Fig. 2d). Overall, DICA showed a higher abundance of nitrogen-containing metabolites associated with assimilation, while ERLE showed a greater presence of nitrogen carriers.

In the rhizosphere, nitrogen-containing metabolites revealed by FT-ICR-MS were more abundant under closed canopy conditions, regardless of plant species (Fig. 2e; Kruskal-Wallis test: χ² = 0.400, p = 0.527 for plant; χ² = 11.247, p = 0.0008 for canopy; Fig. 2e). This finding suggests that rhizospheres, while critical interfaces of plant-soil interactions, are substantially influenced by micro-environmental factors, such as the nitrogen fixed by mesquite trees.

The nominal oxidation state of carbon (NOSC) of FT-ICR-MS features, which reflects substrate quality and availability for rhizosphere microorganisms³⁶, provided further insights into microbial activities. Under open canopy conditions, features unique to ERLE rhizospheres had slightly but significantly lower NOSC than those unique to DICA (Dunn’s test: z-value = –2.28, p.adj = 0.0447). This lower NOSC suggests a depletion of easily accessible carbon sources, potentially indicating a more active rhizosphere microbial community around ERLE roots in these less favorable conditions. Conversely, under closed canopy, DICA-specific features exhibited a pronounced lower NOSC than ERLE-specific feature (Dunn’s test: z-value = 13.5, p.adj <0.0001), suggesting greater microbial consumption of labile carbon sources in the DICA rhizosphere under these improved environmental conditions (Fig. 2f). Furthermore, enrichment analysis revealed most metabolite subclasses were enriched in DICA rhizospheres, except for fatty amides and amines under closed canopy (Fig. 2g). These findings challenge our initial hypothesis that ERLE fosters a more active rhizosphere, suggesting that the native DICA supports greater microbial diversity and activity under favorable conditions.

DICA produces abundant defensive metabolites; ERLE introduces potential novel weapons

To test the hypothesis that invasive Lehmann lovegrass (ERLE) invests less in defensive metabolites than native Arizona cottontop (DICA), we investigated the presence, abundance, and spatial distribution of these metabolites in both root and rhizosphere metabolomes. Our chemical ecology investigation revealed distinct defensive strategies between ERLE and DICA that largely aligned with our hypothesis. FT-ICR-MS analysis revealed a higher presence of recalcitrant metabolites, such as lignin-like (Fig. 3a1; Kruskal Wallis test: χ² = 20.757, p < 0.001 for plant; χ² = 1.938, p = 0.164 for canopy) and condensed hydrocarbon-like metabolites (Fig. 3a2; Kruskal Wallis test: χ² = 25.629, p < 0.001 for plant; χ² = 0.006, p = 0.937 for canopy), in DICA roots compared to ERLE.

**Fig. 3: Root and rhizosphere secondary metabolites associated with plant defense.**

Further analysis using LC MS/MS and ¹H NMR data focused on specific chemical classes known for their defensive roles, such as benzene, phenols, flavonoids, cinnamic acids³⁷, azole³⁸, fatty acyls³⁹, prenol lipids⁴⁰ and indoles⁴¹. In most of these classes, more defensive features were upregulated in DICA compared to ERLE, with the exception of azoles and prenol lipids (Fig. 3b). Specifically, we identified four metabolites (hydroperoxyoctadecadienoic acid, gentisic acid, hesperetin, and 3−Indoleacrylic acid) within these classes with significant differential expressions, annotated chemical names, and confirmed roles in plant defense from literature^42,43,44,45. Among these, only hesperetin (a flavonoid; Fig. 3b3; Kruskal-Wallis test: χ2 = 18.243, p < 0.0001 for plant; χ2 = 0.064, p = 0.800 for canopy) was annotated with confirmed ID. All four metabolites were consistently upregulated in DICA roots compared to ERLE roots, regardless of canopy conditions (Fig. 3b1; Kruskal-Wallis test for hydroperoxyoctadeca-9,11-dienoic acid (HpODE; a fatty acyl): χ2 = 22.523, p < 0.0001 for plant; χ2 = 0.169, p = 0.681 for canopy. Fig. 3b2; Kruskal-Wallis test for gentisic acid (a benzene derivative): χ2 = 25.306, p < 0.0001 for plant; χ2 = 0.009, p = 0.924 for canopy. Fig. 3b4; Kruskal-Wallis test for 3-Indoleacrylic acid (an indole): χ2 = 25.626, p < 0.0001 for plant; χ2 = 0.196, p = 0.658 for canopy).

Untargeted MALDI-MSI revealed the spatial distribution of defensive metabolites within root sections. By matching candidate metabolites identified by MALDI-MSI with those identified by LC MS/MS and ¹H NMR, we annotated two defensive metabolites – hesperetin and 3-Indoleacrylic acid– that were enriched in the outer epidermal layer of both species (Fig. 2c), with significantly higher levels in DICA (Kruskal-Wallis for hesperetin test: χ² = 12.403, p < 0.001 for plant; χ² = 0.75, p = 0.387 for canopy; Kruskal-Wallis test for 3−Indoleacrylic acid: χ² = 5.357, p = 0.021 for plant; χ² = 0.003, p = 0.953 for canopy). Guided by this pattern, we targeted 65 additional MALDI-MSI features exhibiting similar epidermal localization (Supplementary Note 4). Interestingly, chelidonic acid was specifically enriched in ERLE root epidermis, hinting at a unique chemical defense mechanism in the invasive species. Additionally, non-gap-filled LC MS/MS results identified six potential defensive metabolites, including the confirmed annotation gamma-linolenic acid, among 19 features unique to ERLE (Supplementary Data 1), supporting the hypothesis that ERLE produces unique defensive metabolites that may contribute to its invasive success.

Rhizosphere metabolomes identified from LC MS/MS and ¹H NMR also aligned with our hypothesis. Although the number of significantly different features was lower than in roots, (Fig. 2f), we observed consistent upregulation of defensive metabolites, including HpODE, indole-3-aldehyde, and 3-indoleacrylic acid, in the DICA rhizosphere across both canopy conditions (Fig. 3d; Kruskal-Wallis test for HpODE: χ2 = 11.036, p = 0.0008 for plant, χ2 = 2.604, p = 0.107 for canopy; for indole-3-aldehyde: χ2 = 15.391, p < 0.0001 for plant, χ2 = 3.367, p = 0.067 for canopy; for 3-indoleacrylic acid: χ2 = 13.938, p = 0.0002 for plant, χ2 = 0.049, p = 0.825 for canopy). Given the role of fatty acyls (like HpODE) as signaling molecules in defensive responses (Svoboda & Boland, 2010), this pattern suggests that ERLE may experience less pressure from belowground enemies despite lower investment in root defenses. Moreover, the higher number of defensive metabolites in root than rhizosphere may suggest a preference for retaining most defensive metabolites within roots for both DICA and ERLE.

ERLE shows greater phenotypic plasticity in root and rhizosphere metabolomes

To investigate the hypothesis that Lehmann lovegrass exhibits greater phenotypic plasticity in its root and rhizosphere metabolomes in response to varying environmental conditions, we investigated the change of root and rhizosphere metabolomes between open and closed canopy conditions for both ERLE and DICA. Consistent with our hypothesis, ERLE demonstrated a greater shift in metabolome between different environmental conditions, particularly in the rhizosphere. Root metabolite diversity measured by FT-ICR-MS did not differ significantly between canopy conditions for either species (Fig. 4a; Kruskal Wallis test: χ² = 0.169, p = 0.681 for plant; χ² = 0.677, p = 0.411 for canopy). However, ERLE had fewer metabolites shared between open and closed canopy (Fig. 4b) and a larger proportion of the variance in root metabolomes ordination explained by canopy for most platforms compared to DICA (Fig. 4c), indicating a stronger influence of the environment on ERLE root metabolism. Enrichment analysis revealed that more subclasses of metabolites in ERLE roots were enriched under open canopy (Fig. 4d), further supporting a more significant shift in metabolomes in response to environmental changes than DICA.

**Fig. 4: Root and rhizosphere metabolome responses to canopy conditions.**

Greater phenotypic plasticity of ERLE than DICA was also evident in the rhizosphere metabolome. Rhizosphere metabolite diversity was significantly higher for ERLE under closed canopy compared to open canopy, while no significant difference was observed for DICA (Fig. 4e; Dunn’s test: z-value = –2.75, p.adj = 0.036 for ERLE-closed vs ERLE-open; Dunn’s test: z-value = –1.59, p.adj = 0.337 for DICA-closed vs DICA-open). Consistent with the root metabolome findings, ERLE had fewer shared metabolites between open and closed canopy compared to DICA (Fig. 4f), and canopy conditions explained more variance in ERLE rhizosphere metabolomes across most platforms (Fig. 4g).

To specifically assess phenotypic plasticity in root exudate profiles, we focused on a subset of 397 rhizosphere features detected in both root and rhizosphere metabolomes by LC MS/MS and ¹H NMR, under the assumption that true root exudates should be present in both locations. While DICA generally had more upregulated root exudates than ERLE (Fig. 4h1), ERLE exhibited a more substantial shift in root exudate profiles between open and closed canopy conditions (Fig. 4h2). This differential response suggests that ERLE may be more adept at adjusting its root exudate composition depending on environmental conditions. Further, this chemical flexibility may provide ERLE with a competitive advantage in varying environments, contributing to its invasive success.

Discussion

Our multimodal metabolomics approach reveals how different platforms enhance understanding of invasion mechanisms at the molecular level. Each method’s strengths—FT-ICR-MS for broad metabolite coverage, LC MS/MS for relative abundance quantification, and ¹H NMR for concentration of small molecules—provided comprehensive cross-validation. Across all platforms, we observed consistent compositional patterns between DICA and ERLE in the Santa Rita Experimental Range (SRER): root metabolomes were primarily shaped by plant species identity, while rhizosphere metabolomes were influenced by both canopy status and plant species, with canopy being the dominant factor. These consistent patterns across platforms strengthen confidence in our findings.

Each method provided distinct mechanistic insights. While LC MS/MS indicated higher relative abundance of nitrogen-containing metabolites in DICA roots, the superior coverage of FT-ICR-MS revealed a higher percentage of these metabolites in ERLE, aligning with previously observed physiological patterns¹⁵. Similarly, the integration of FT-ICR-MS, LC MS/MS, and MALDI-MSI spatial analysis provided a multifaceted view of defense investment, from structural metabolites to specialized metabolites and their tissue-specific distributions. For root exudates, the combination of highly reliable ¹H NMR with the broader coverage of LC MS/MS yielded a comprehensive profile that neither method alone could achieve.

Our analysis uncovered three key strategies potentially contributing to ERLE’s dominance over the native DICA: (1) more efficient nitrogen allocation to shoot growth, driven by higher intensity of shoot nitrogen carriers and lower intensity of root metabolites involved in nitrogen assimilation compared to DICA; (2) reduced investment in root defensive metabolites, with lower abundances of most defensive compounds but emergence of a potential novel defense metabolite; (3) greater metabolic plasticity in response to environmental stress, evidenced by more evident changes in its root and rhizosphere metabolome in response to canopy conditions. These findings extend beyond previous trait-based studies by providing a molecular-level understanding of how invasive plants can simultaneously optimize resource acquisition, defense, and stress response.

The first strategy is associated with the hypothesis that invasiveness is linked to superior nitrogen acquisition and utilization^46,47, including higher leaf nitrogen content, greater biomass under high nitrogen availability, and enhanced soil nutrient-releasing enzyme activities^48,49,50. Our findings align with previous reports of lower root C:N ratios in ERLE compared to DICA¹⁵, indicating potentially greater nitrogen uptake by ERLE. However, traditional metrics often fail to elucidate how nitrogen is utilized within the plant. We found that DICA invested nitrogen in defense and root exudates, while ERLE conserved nitrogen for shoot growth.

Our results suggest that ERLE might adopt a more conservative nitrogen use strategy, primarily allocating nitrogen to leaves and replenishing the soil nitrogen pool through defoliation of high-nitrogen leaves to maintain dominance in nutrient-poor environments. In contrast, DICA aggressively invests acquired nitrogen in growth and competitive advantage, explaining its prevalence in fertile islands, positive response to fertilizers⁵¹, and rapid regrowth after grazing¹⁶. The investment in root exudates in DICA may foster a potentially more diverse and active rhizosphere microbial community⁵². The substantial NOSC decrease for DICA-unique metabolites under closed canopies, likely reflects a complex interaction between mesquite’s N-fixing symbionts, DICA’s root exudates, and the overall enhanced nutrient status of fertile islands⁵³. This three-way interaction may enable rhizosphere microbes to utilize more complex but energetically rewarding metabolic pathways, potentially explaining DICA’s enhanced performance in the nutrient-rich microsites despite showing lower absolute nitrogen content than ERLE³⁶. It’s important to note, however, that metabolomics analyses remain largely descriptive, and our current approach has limitations that constrain interpretation of nitrogen fate within plants (see Supplementary Note 5 for a detailed discussion of limitations).

The second strategy is associated with the enemy release hypothesis, a frequently cited yet often debated explanation for the dominance of exotic plants⁵⁴, which posits that invasive species experience reduced impact from enemies in their introduced range⁵⁵. Traditional tests of this hypothesis often rely on comparing enemy presence or abundance⁵⁶. However, such approaches may not capture the actual impact of enemies, especially in the case of our study species, forage grasses known for their grazing tolerance¹⁴. DICA is known for its high palatability to livestock⁵⁷, leading to the hypothesis that preferential grazing might contribute to the dominance of ERLE at the SRER. Nevertheless, a 35-year monitoring at SRER found no significant difference between grazed and ungrazed pastures on the spread of ERLE¹², suggesting that livestock may not be the primary enemy In the case of DICA and ERLE. Although grasses have traditionally been considered to lack diverse secondary metabolites and thus unsuitable for biological control⁵⁸, our metabolomic analyses revealed rich root and rhizosphere secondary metabolite profiles, suggesting potential for belowground biocontrol approaches⁵⁹.

While DICA showed a higher investment in defensive metabolites than ERLE across all methods, they still exhibited higher levels of defensive signaling compounds, including octadecanoids³⁹, gentisic acid in roots⁴⁵, and indole-3-aldehyde in rhizospheres⁶⁰. This suggests that DICA may be experiencing greater pressure from belowground enemies, potentially rendering its increased investment in defense worthwhile. This aligns with the general trend of negative plant-soil feedbacks for grasses, which are often less negative for invasive species⁶¹. Conversely, ERLE’s apparent reduced enemy pressure may result from its production of novel defensive metabolites. For instance, chelidonic acid, known to increase in response to pathogen infection despite its unclear defensive role⁶², and gamma-linolenic acid, associated with pathogen resistance variant of cabbage⁶³, were exclusively present in ERLE. The exclusive presence of these less-documented defensive metabolites in ERLE suggests they may act as novel weapons in its introduced range and aiding its invasive success⁶⁴. Together, these findings underscore distinct defensive strategies between ERLE and DICA, with DICA relying on abundant defensive metabolites and ERLE potentially utilizing novel metabolites to enhance its invasive success.

It is important to note that defensive metabolites, particularly when examined as chemical classes, often serve multiple functions beyond pathogen defense, including roles as growth regulators and even primary metabolites⁶⁵. For example, hesperetin exhibits both anti-pathogenic properties⁴² and can induce beneficial soil microorganism associations in certain legume plants⁶⁶. Determining the explicit functions of these secondary metabolites requires further controlled experiments, such as direct assessment of defensive metabolite effects through pathogen cultivation assays, and measurement of defensive metabolite production in plants following pathogen inoculation (see Supplementary Note 5 for a detailed discussion of further validation experiment).

The third strategy is associated with phenotypic plasticity, or the ability to express different traits in response to environmental variation⁶⁷. This hypothesis suggests that high levels of phenotypic plasticity not only enable introduced plants to establish in the introduced environments⁶⁸ but also to outcompete native plants under fluctuating conditions⁶⁹. Traditionally, this hypothesis has been tested by comparing the variance of plant traits between populations of co-occurring native and invasive plants⁶ or by examining changes in plant traits in response to varying stress levels⁷. In our study, we investigated phenotypic plasticity by comparing metabolite diversity, richness, and composition between the favorable closed canopy and the more stressful open canopy. We found that ERLE exhibited greater plasticity in richness and composition of both root and rhizosphere metabolomes, aligning with observations of larger plasticity in other plant traits of invasive plants⁷.

While plant traits can often be directly linked to fitness, the relationship between fitness and metabolites is more nuanced. The larger difference in root metabolites between open and closed canopy conditions in ERLE, compared to DICA, suggests two possibilities: either DICA is inherently better adapted to open canopy conditions and requires fewer metabolic adjustments, or ERLE exhibits greater plasticity, altering its metabolism more extensively in response to adverse conditions. Many annotated root exudates upregulated in DICA or in ERLE under open canopy, such as organic acids like butyrate, malate, succinate, acetate, and citrate, are induced by drought or nutrient deficiency^70,71,72. This might suggest that ERLE employs a “Jack-of-all-trades” strategy (i.e., exhibit superior survival under stressful environments), conserving resources under favorable conditions and activating specific stress responses to increase its fitness under unfavorable conditions. In contrast, DICA appears to maintain a constitutive stress response, which may confer greater tolerance to environmental fluctuations but could also entail unnecessary resource expenditure under benign conditions.

Our multi-platform metabolomics approach has provided insights into the complex factors contributing to ERLE’s invasion success, offering a fresh perspective on plant competition and highlighting the importance of climate-driven changes in resource distribution. Our results do not support the notion that direct competition between DICA and ERLE is the primary driver of ERLE’s dominance and DICA’s decline at SRER since the 1960s. Contrary to expectations, DICA’s greater investment in root exudates and defensive metabolites, along with its more active rhizosphere microbiome, suggest a competitive advantage particularly under closed canopy^15,73. Instead, our findings suggest that ERLE’s dominance is more likely a byproduct of increasing aridity⁷⁴ with drier conditions creating increasingly patchy resource distributions that favor ERLE’s conservative nitrogen utilization and phenotypic plasticity^29,75. Under harsh open canopy conditions, ERLE increases root exudate production and prioritizes shoot nitrogen allocation, redirecting resources toward stress tolerance and rapid colonization of nutrient-poor soils. In contrast, under closed canopy conditions, ERLE reduces exudate production and defensive investments, conserving energy for environmental fluctuations and future stress events. This collective strategy—balancing resource allocation between harsh and favorable conditions—enhances ERLE’s ability to exploit heterogeneous nutrient availability in arid systems, whereas DICA’s inflexible high-resource strategy becomes maladaptive as aridity increases.

These insights demonstrate the value of multimodal metabolomics in offering additional insights into invasion dynamics, complementing trait-based approaches by comprehensively profiling plant inner biological processes, plant-microbe interactions, nutrient uptake and utilization, and stress responses. This approach allows for testing multiple hypotheses grounded in functional traits and offers the potential to examine multiple co-occurring species within the same system, providing a more comprehensive understanding of community-level interactions. This observed shift towards ERLE dominance and the concomitant decrease in forage value may represent a form of land degradation driven by climate change⁷⁶. Consequently, site-level management strategies may be limited in controlling ERLE, and interventions like adding labile carbon to stabilize soil nitrogen⁷⁷ could further disadvantage DICA. Given these findings, more effective management strategies might focus on maintaining or enhancing soil moisture retention and nutrient stability, potentially through practices that increase plant diversity and soil organic matter content. Conversely, manual removal and herbicide application might even exacerbate the problem by increasing barren ground and further depleting soil nutrients.

As climate change increases the aridity in drylands worldwide, which predicted to further increase the heterogeneity of resources and stress between islands and interspaces^29,78, our study predicts that species with traits similar to ERLE—conservative resource use and high phenotypic plasticity—may become increasingly dominant, potentially leading to widespread shifts in community composition and ecosystem function. While our multi-platform metabolomics approach has yielded valuable insights, there remain methodological challenges and areas for future development to fully harness the potential of metabolomics in invasion ecology and trait ecology (Supplementary Note 5). Key limitations include the difficulty in comparing metabolomes across species due to inherent variations in metabolic profiles and resource acquisition rates, as well as challenges in quantifying metabolites given variable ionization efficiencies and distinguishing plant-derived from microbial metabolites. Future validating experiments using isotope labeling and controlled soil conditions, coupled with advanced targeted metabolomics and imaging techniques, will be crucial for addressing these limitations and strengthening metabolomic applications in ecological research. Overall, our study highlights the immense potential of metabolomics in expanding the horizons of trait-based ecology and unraveling the complex mechanisms underlying plant invasions and informing future management strategies in a changing world.

Methods

Study site and experimental design

This study was conducted at an elevation of approximately 1125 meters within the Santa Rita Experimental Range (SRER) in Arizona (Fig. 5a). The SRER is situated at the northwestern base of the Santa Rita Mountains, approximately 55 kilometers south of Tucson, Arizona, with the soil series Sasabe-Baboquivari complex and the ecological site R041XC319AZ—Sandy Loam Upland 12–16⁷⁹. The elevation within the SRER ranges from 900 to 1450 m, with annual precipitation varying from 275 to 450 mm and the biome ranges from a desert scrub at the lowest elevations to savanna woodlands at the highest¹⁴. The landscape consists of open grasslands consist of perennial grasses, cacti, succulents, and seasonal herbaceous species, interspersed with woody patches dominated by velvet mesquite (Prosopis velutina Woot.), catclaw acacia (Acacia greggii Gray), and blue palo verde (Cercidium floridum Benth.)¹⁴. Established in 1902, the SRER was created to understand ecosystem dynamics and apply sustainable management practices, making it one of the oldest continuously operating rangeland research facilities globally⁸⁰. During the monitoring period, the most noticeable changes over the past 100 years have been the steady increase in mesquite tree density and the increasing dominance of nonnative Lehmann lovegrass¹⁴. In 2006, a new grazing scheme was implemented, emphasizing avoiding re-grazing of plants during the growing season and adjusting dormant season grazing capacity based on summer production⁸¹. Our study site, pasture 2S, was not grazed in 2022 but experienced intense grazing for one month in 2021, resulting in approximately 60% of grass biomass having been consumed.

Lehmann lovegrass (Eragrostis lehmanniana; ERLE) is a warm-season, perennial bunchgrass initially introduced to SRER in 1937 from South Africa for range restoration and high-quality forage grass. It quickly spread and outcompeted native perennial grasses and shrubs⁸². Lehmann lovegrass is less palatable than native perennial grasses during the summer growing season, though its ability to establish under adverse conditions, rapid recovery after disturbances such as fire and grazing, and high-quality winter forage compared to other native grasses make it a desirable forage grass⁸³. Arizona cottontop (Digitaria californica; DICA) is native, warm-season perennial bunchgrass and high-value forage grass due to its high palatability and tolerance to grazing⁵¹. It is a climax dominant in the semidesert grassland, but its density has decreased steadily since the invasion of Lehmann lovegrass began in 1970¹⁴.

Three mature velvet mesquite trees (Prosopis velutina) were selected from pasture 2S for our study. Each tree had ERLE and DICA occurring together both under the canopy and in nearby inter-canopy regions. Each focal mesquite tree, along with a surrounding area of 10 m in diameter, constituted a single site, with each tree spaced approximately 35 m apart (Fig. 5a). At each tree, four plot types (1 m × 1 m) were established to represent two dominant patch types (ERLE and DICA) and two canopy conditions (closed and open), resulting in a total of 12 plots. Three individual plants were randomly selected and uprooted from each plot, which sum up to total of 36 samples for each root or rhizosphere metabolome.

To ensure robust representation while controlling for spatial variability, we employed a nested sampling design. Each mesquite tree served as a block, with three individual plants sampled per species per canopy condition within each block. This design provides nine total replicates per treatment combination. We specifically sampled from established monoculture patches to minimize interference from neighboring species. Rhizosphere soil was collected from 2 to 8 cm depth, where root density was highest and plant’s influence on soil properties was most pronounced. Preliminary analyses showed that variation between blocks was minimal compared to treatment effects, supporting the effectiveness of this sampling approach in capturing consistent biological patterns while controlling for spatial heterogeneity.

Closed plots were located beneath the dripline of the mesquite canopy, while open plots were positioned outside the dripline but no more than 10 m radially from the focal mesquite. In arid ecosystems, the closed canopy conditions represent “fertile islands” characterized by higher fertility, microbial activity, and reduced abiotic stress, providing a contrasting environment to open canopy conditions^29,84. Open canopy areas are often infertile and support less biomass, experiencing not only lower litter input and higher temperature and water stress, but also lose surface soil through erosion, with these resources accumulating in nearby fertile islands⁸⁵. In our study site, despite relatively large variation between mesquite trees and grass patches, soil under open canopy often had lower NO₃⁻ and PO₄²⁻ concentrations compared closed canopy areas (Fig. 5c).

Sample collection and processing

Root and rhizosphere soil samples were collected on July 5, 2022. Rhizosphere soils were obtained by gently agitating the soil surrounding the roots from a depth of 2–8 cm below the soil surface. After the rhizosphere soil was extracted, roots were carefully separated from the plants. Both the soil and root samples were immediately placed in a cooler and transported to the laboratory within 4 h. Soil and root samples were first stored in a −20 °C freezer upon arrival at the laboratory. Prior to long-term storage in a −80 °C freezer, soil samples were processed through a 500-μm sieve, while root samples were cleaned and cut into 1-cm pieces within 3 days. Rhizosphere soil metabolite extractions were carried out between November 30 and December 5, 2022. Root metabolite extractions were conducted between January 27 and February 2, 2023. Samples were stored at −4 °C during and following extraction. Soil and root metabolite extraction protocols were adapted from Tfaily et al.^86,87 for FT-ICR-MS, Hildebrand et al.⁸⁸ for ¹H NMR, and Portman et al.⁸⁹ for LC MS/MS. Briefly, five grams of soil from each sample were separated into 50 ml Falcon tubes containing 20 ml of MilliQ water and briefly vortexed. Three grams of roots, consisting of approximately 20 pieces of 1-cm segments to account for metabolite variation within individual root systems, were homogenized in 10 ml of MilliQ water using an electric probe homogenizer (Fig. 5d). We employed water extraction to simulate natural conditions where rain events mobilize metabolites, focusing on compounds that would be biologically available in the soil-root interface. While organic solvent extractions might capture a broader range of semi-polar defensive compounds, water extraction better reflects the natural soil solution conditions and plant-soil interactions. Subsequently, all samples were sonicated in a water bath for 30 min at 20–30 °C (FisherBrand CPX3800) and centrifuged for 10 min at 500 rpm (Eppendorf Centrifuge 5430). The supernatant was then transferred to a new Falcon tube, centrifuged for a second time, and passed through a Minisart syringe filter with a 0.45 µm pore size. The filtered extracts were processed differently for each analytical platform. For LC MS/MS analysis, 3 ml of supernatant from root samples and 6 ml from rhizosphere samples were transferred to two 2 ml glass autosampler vials and dried using a vacuum centrifuge (Eppendorf Vacufuge plus). For FT-ICR-MS analysis, additional 1.5 ml and 3 ml of root and rhizosphere supernatants, respectively, were used for sample clean-up via solid-phase extraction (SPE) (Bond Elut PPL SPE Cartridges). This clean-up step was necessary to remove salts and concentrate samples for direct injection. Samples were acidified to a pH of 2–3 using HCl (Fisher Chemical) and loaded onto SPE cartridges under vacuum filtration. The cartridges were then washed with 15,000 µl of 0.01 M HCl solution, dried, and eluted into 2 ml autosampler vials with 1500 µl of methanol (MeOH, VWR chemicals, LLC). Vials were stored at −80 °C before sending for FT-ICR-MS analysis. The remaining aqueous extracts were used for ¹H NMR analysis. All processed samples were stored at −80 °C before shipping. All metabolite extracts were sent to the Environmental Molecular Sciences Laboratory in Richland, WA for FT-ICR-MS, ¹H NMR, LC MS/MS, and MALDI-MSI analysis.

Metabolomics

Different metabolomic techniques have unique strengths and limitations in their ability to detect various aspects of metabolites^90,91,92. For example, Fourier transform ion cyclotron resonance (FT-ICR-MS) is ideal for revealing the overall composition of the metabolome due to its high coverage (200–1200 Da) and the ability to separate peaks with similar m/z ratios²⁵, but it only provides molecular formula assignments, thus limiting its ability to directly predict specific biological processes. In contrast, liquid chromatography tandem mass spectrometry (LC MS/MS) provides structural elucidation on metabolites ranging from 50 to 1000 Da²⁶. Liquid state nuclear magnetic resonance (¹H NMR) offers structural and quantitative data on small metabolites that fall below the detection range of FT-ICR-MS and LC MS/MS, particularly root exudates²⁴. Additionally, matrix-assisted laser desorption/ionization mass spectrometry imaging (MALDI-MSI) allows visualization of the spatial distribution of metabolites within plant roots²⁸. By strategically integrating data from these diverse metabolomic techniques, we achieve a comprehensive and holistic understanding of the plant-soil metabolomes and generate complementary datasets that can be integrated and cross-validated, increasing confidence in the identified metabolites and metabolic patterns associated with plant invasion success.

The extracted metabolites were analyzed using Fourier transform ion cyclotron resonance mass spectrometry (FT-ICR-MS) following the protocol developed by Tfaily et al.^86,87. Root and rhizosphere metabolites were obtained using a 12 Tesla (12 T) Bruker SolariX Fourier transform ion cyclotron resonance mass spectrometer (FT-ICR-MS). Samples were directly injected into the instrument using a custom automated direct infusion cart, which performed two offline blanks between each sample to minimize contamination. Equipped with a standard electrospray ionization (ESI) source, the FT-ICR-MS operated in negative mode with the needle voltage set to +4.0 kV. Mass data were collected over a range of 150 m/z–1000 m/z at a resolution of 8 M. For each sample, three hundred scans were co-added and internally calibrated using OM homologous series separated by 14 Da (–CH2 groups). The mass measurement accuracy typically remained within 1 ppm for singly charged ions across a broad m/z range of 150 m/z–1100 m/z. Subsequently, Bruker Data Analysis (version 5.0) was utilized to convert raw spectra to a list of m/z values by applying the FTMS peak picker module with a signal-to-noise ratio (S/N) threshold set to 7 and an absolute intensity threshold to the default value of 100. Due to the non-linear nature of baseline noise and the need to process multiple samples in bulk, we applied a strict S/N threshold, using the lowest acceptable intensity value to ensure peak selection was based solely on the S/N ratio. This approach ensured consistency and reproducibility across datasets. Chemical formulae were then assigned using Formularity⁹³, an in-house software, following the Compound Identification Algorithm^87,94,95. Chemical formulae were assigned based on the following criteria: S/N ≥ 7, and mass measurement error ≤0.5 ppm, taking into consideration the presence of C, H, O, N, S, and P while excluding other elements.

Proton nuclear magnetic resonance (¹H NMR) spectroscopy was performed following the protocol described by Hildebrand et al.⁸⁸. Metabolite extracts (180 µL) were combined with 2,2-dimethyl-2-silapentane-5-sulfonate-d6 (DSS-d6, Chenomx Inc.) in D2O (20 µL, 5 mM) and thoroughly mixed prior to transfer to 3 mm NMR tubes. NMR spectra were acquired on a Bruker Neo spectrometer operating at 18.8 T (¹H ν0 of 800.30 MHz) equipped with a 5 mm Bruker TCI/CP HCN (inverse) cryoprobe with Z-gradient at a regulated temperature of 298.0 K. The 90° ¹H pulse was calibrated prior to the measurement of each sample. The one-dimensional ¹H spectra were acquired using a Nuclear Overhauser Effect Spectroscopy (NOESYpr1d) pulse sequence with a spectral width of 20.1 ppm and 2048 transients. The NOESY mixing time was 100 ms and the acquisition time was 4 s followed by a relaxation delay of 1.5 s during which presaturation of the water signal was applied. The 1D ¹H spectra were manually processed, assigned metabolite identifications, and quantified using Chenomx NMR Suite 9.0. Time domain free induction decays (72,114 total points) were zero filled to 131,072 total points before Fourier transform, followed by exponential multiplication (0.3 Hz line-broadening) and semi-automatic multipoint smooth segments baseline correction. Chemical shifts were referenced to the ¹H methyl signal in DSS-d6 at 0 ppm. Metabolite identification was based on matching the chemical shift, J-coupling and intensity of experimental signals to metabolite signals in the Chenomx library. Quantification was based on fitted metabolite signals relative to the internal standard (DSS-d6). Signal to noise ratios (S/N) were measured using MestReNova 14.1, with the limits of quantification and detection equal to an S/N of 10 and 3, respectively. Standard 2D experiments such as ¹H/¹³C—heteronuclear correlation (HSQC) or 2D ¹H/¹H Total Correlation spectroscopy (TOCSY) further aided corroboration of several metabolite identifications where there was sufficient S/N. Structure-based chemical taxonomy was assigned using ClassyFire⁹⁶.

Liquid chromatography electrospray ionization tandem mass spectrometry (LCESIMS/MS) used a Thermo Vanquish Flex UHPLC system interfaced with a Thermo QExactive HF-X mass spectrometer. Parameters and protocols were adapted from DiDonato et al.⁹⁷. Full MS scan data are acquired at a resolving power of 140,000 FWHM at m/z 200 and the scanning range of m/z 80–800. Samples are analyzed in both positive and negative ionization modes using higher-energy collision dissociation. Spray voltage 3.7 or 3.0 kV for positive and negative modes respectively; capillary temperature 350 °C. A pooled quality control (QC) sample, created by combining aliquots from each biological sample, was used to normalize ionization levels and correct systematic biases across the large-scale metabolomics dataset.

Reverse phase (RP) chromatography was used to analyze hydrophobic compounds. 5–10 µl of metabolite extracts were injected and compounds separated using a Thermo Hypersil GOLD column (2.1 × 150 mm, 3 μm particle size) with a column temperature of 40 °C and a flow rate of 400 μL min–1. Mobile phase A (Water with 0.1% formic acid, Fisher Scientific International, Inc.) and B (acetonitrile (ACN) with 0.1% formic acid, Fisher Scientific International, Inc.) were initially 90:10, respectively. The gradient method continued as follows: from 0–2 min, the gradient was held at 90% A, followed by a decrease to 10% A from 2 to 11 min. At 11 min, 10% A was held until 12 min, after which the flow rate was increased to 500 μL min⁻¹ from 12–12.5 min. The gradient was then returned to 90% A, maintaining a flow rate of 500 μL min⁻¹ from 12.5–13.5 min, held until 14 min. From 14–14.5 min, the flow rate was reduced to 400 μL min⁻¹, and the method concluded with a hold at 90% A from 14.5 to 16 min. The column was recalibrated during the chromatographic run between 12.5 and 16 min between samples. Hydrophilic interaction liquid chromatography (HILIC) was used to analyze polar compounds. 5–10 µl of metabolite extracts are injected and compounds separated using a using a ACQUITY UPLC BEH HILIC column (2.1 × 100 mm, 1.7 μm particle size) with a column temperature of 50 °C and a flow rate of 300 μL min^–1. Chromatography was adapted from clendinen et al.⁹⁸. Mobile phase A (5% ACN: 95% 10 mM NH₄OAc in H2O with 0.05% NH₄OH, Sigma Aldrich) and B (100% ACN with 0.05% NH₄OH, Sigma Aldrich) were initially 5% and 95%, respectively, at a flow rate of 300 μL min^–1. The gradient method continued as follows: from 0 to 6.0 min, A was increased to 63%. This was maintained until 7.0 min, after which the gradient quickly reduced to 5%A between 7.0 and 7.1 min. From 7.1 to 7.2 min, the flow rate increased to 500 μL min⁻¹ while holding at 5% A, which was then maintained from 7.2 to 9.5 min. Between 9.5 and 9.7 min, the flow rate decreased to 300 μL/min, and the gradient held steady at 5% A until 12.0 min. The column was recalibrated during the chromatographic run between 7–12 min between samples. Data-dependent acquisition MS/MS was performed only on pooled samples.

For LC MS/MS data processing, we used Thermo Compound Discoverer 3.3 to confidently identify metabolites. Spectra were aligned using an adaptive curve with a maximum RT shift of 0.3 (RP) or 0.6 (HILIC) minutes and a 3-ppm mass tolerance. Peaks were selected based on a minimum intensity of 1e6 and a chromatographic signal-to-noise ratio (S/N) of 3. These thresholds were determined through manual inspection of chromatograms to achieve optimal data clarity, using the typical S/N ratio of 3 to ensure clean and reliable peak selection. Detected features were grouped with a 3-ppm mass tolerance and a 0.25-min RT tolerance, and features that did not meet the sample count threshold were filtered out. Gap filling by default was applied to the remaining features to detect any low-abundance metabolites that might have been missed, either by redetecting at a lower threshold or filling with background noise when appropriate. We used gap filling to enhance downstream data analysis by reducing zero-inflated data and ensuring a more complete profile of metabolites, which is especially meaningful in rhizosphere samples due to their spatial proximity, where metabolites are likely to be consistently present. To assess any metabolites potentially unique to one species, we also performed an analysis without gap filling by removing gaps filled with background noise.

Compounds were assigned using isotopic patterns, RT, MS1, and/or MS2 data. All identifications and integrated peaks were manually validated. We employed a six-level confidence scale for identification: 1. Confirmed ID: MS1, MS2, and RT match an internal database entry; 2. Match with Mass and RT only: MS1 and RT match an internal database entry; 3. Not in internal database: MS2 matches external databases, but not internal; 4. Potential ID: Primarily matches MS1 and partial MS2 match external databases; 5. MS2: Has MS2 spectrum but no database match; 6. No MS2: Lacks MS2 spectrum. The chemical taxonomy of annotated metabolites was determined using ClassyFire with InCheKey based on the putative chemical names⁹⁶. InChIKeys were generated through the PubChem Identifier Exchange Service (https://pubchem.ncbi.nlm.nih.gov/docs/identifier-exchange-service). Batch processing for ClassyFire was facilitated using the ClassyFire Batch Tool by FiehnLab (https://cfb.fiehnlab.ucdavis.edu/). For features with MS2 but lacking annotation, we performed structural computations using SIRIUS graphical user interface (version 5.8.6), followed by ClassyFire assignment using Canopus⁹⁹.

Matrix-assisted laser desorption/ionization Fourier-transform ion cyclotron resonance mass spectrometry imaging (MALDI-MSI) randomly selected three root samples per treatment (12 total). Following removal of soil particles, roots were cut into 2 cm segments and transferred to 1.5 ml Eppendorf tubes before shipping to the Environmental Molecular Sciences Laboratory in Richland, WA. Roots were embedded within a mixture of 7.5% Hydroxypropyl methyl-cellulose (viscosity 40–60 cP, 2% in water at 20 °C, Sigma-Aldrich) and 2.5% polyvinylpyrrolidone (average mol weight 360,000, Sigma-Aldrich), cryosected at −10 °C using a CryoStar NX-70 Cryostat (Thermo Scientific, Runcorn, UK), and 12 µm sections were subsequently thaw-mounted onto ITO coated slides. Sections of twelve roots, where each root was represented by the two serial sections, were distributed over the two slides, and three such sets (for three different MALDI-MSI analyses) were prepared. Three different MALDI-MSI analyses were performed on replicate slides for maximizing molecular coverage: (1) positive polarity mode analysis using 2,5-dihydroxybenzoic acid (DHB, Sigma-Aldrich) as MALDI matrix (2) negative polarity mode analysis using N-(1-naphthyl) ethylenediamine dihydrochloride (NEDC, Sigma-Aldrich) as MALDI matrix, and (3) 4-(2-((4-bromophenethyl)dimethylammonio)ethoxy)benzenaminium bromide (4-APEBA, synthesized in-house) on-tissue chemical derivatization (OTCD) for detection of carbonyl molecules, MALDI matrices and OTCD were applied on the tissue sections using M5 Sprayer (HTX Technologies, Chapel Hill, NC) as described in Zemaitis et al.¹⁰⁰. Briefly, DHB (40 mg/mL in 70% MeOH) was applied using 50 uL/min flow rate and 1200 mm/min linear flow, at 70 °C, over 12 cycles at 3 mm track spacing. NEDC (7 mg/mL in 70% MeOH) was sprayed at 120uL/min flow rate and 1200 mm/min linear flow, at 70 °C, over 8 cycles at 3 mm track spacing. For OTCD, 1-ethyl-3-(3-(dimethylamino) propyl) carbodiimide (6 mg/mL in water, Sigma-Aldrich) and 4-APEBA (2 mg/mL in water) were sequentially deposited at a flow rate 25ul/min, 37.5 °C, 4 cycles.

All imaging analyses were performed on a Bruker Daltonics 12 T solariX FT-ICR-MS, equipped with a ParaCell, Apollo II dual ESI and MALDI source with a SmartBeam II frequency-tripled (355 nm) Nd: YAG laser (Bremen, Germany). Positive and negative ion mode acquisitions were acquired with broadband excitation from m/z 98.3 to 1,000, resulting in a detected transient of 0.5593 s—the observed mass resolution was ~110k at m/z 400. FlexImaging (Bruker Daltonics, v.5.0) was used for the imaging experiments, and analyses were performed with a 35 µm step size. FlexImaging sequences were directly imported into SCiLS Lab (Bruker Daltonics, v.2023.a Premium 3D) using automatic MRMS settings. Ion images were directly processed from the profile datasets within SCiLS Lab, and automated annotation of the centroided dataset was completed within METASPACE against KEGG-v1 database. In case of OTCD, a chemical modifier corresponding to the mass shift expected from 4-APEBA derivatization (+C₁₈H₂₂N₂Br, +345.0966 Da) was used. All annotations were imported back to SCiLS as a new peak list, and average intensity (after root mean square normalization) of each peak across the root section was used for further statistical comparison. In METASPACE, regions of Interest were defined and METASPACE calculated the average intensity of each metabolite within these ROIs. The results, including annotations and average intensities, were exported as CSV files for further statistical analysis.

Statistics and reproducibility

All statistical analyses were conducted using R version 4.2.1¹⁰¹. FT-ICR-MS data were pre-processed using the MetaboDirect pipeline (version 1.04) with binary normalization to target metabolite presence/absence rather than relative abundance³¹, as metabolites are ionized differently and preferential ionization can skew abundance measurements¹⁰². Putative molecular classes were assigned based on O/C and H/C ratios and visualized using Van Krevelen diagrams. The nominal oxidation state of carbon (NOSC) was calculated using the formula ({NOSC}=4-frac{4C+H-2O-3N-2S+5P-z}{C}) to reflect microbial substrate availability and predict microbial activities¹⁰³. The letters represent the stoichiometric number of corresponding elements, z is the charge of molecules.

¹H NMR, LC MS/MS, and MALDI-MSI data were Pareto scaled and FT-ICR-MS data was binary normalized prior to multivariant analysis¹⁰⁴. Although auto-scaling has been shown to yield more most biologically meaningful results of metabolite intensity¹⁰⁴, Pareto scaling can preserve large variations that could be relevant. Binary normalization for FT-ICR MS to be more robust because the data are compositional, and intensities may be influenced by ion competition rather than the actual metabolite concentrations. Permutational multivariate analysis of variance (PERMANOVA) using Jaccard distance (FT-ICR-MS) or Euclidean distance (¹H NMR, LC MS/MS) assessed the effects of plant species and canopy on root and rhizosphere metabolomes using Vegan package¹⁰⁵. To isolate the effects of canopy, PERMANOVA was also performed on subsets of each plant species for ¹H NMR, LC MS/MS, and FT-ICR-MS datasets. The results were visualized using non-matrix distance scaling (NMDS). For FT-ICR-MS data, molecular classes were fitted onto NMDS ordinations using the envfit function from the vegan package to identify those driving differences between plant species and canopy conditions¹⁰⁶.

Kruskal-Wallis tests and log2 fold change calculations on raw intensities were used to identify metabolites differentially regulated between plant species and canopy conditions. Features with |log2 fold change | >1.5 and p < 0.05 were considered significant. Dunn’s test was applied as a post-hoc test to find the differences between groups¹⁰⁷. Enrichment analysis was conducted using ChemRICH to identify ClassyFire subclasses significantly impacted by plant species or canopy effects¹⁰⁸. Briefly, p-values from the Kruskal-Wallis tests on raw metabolite intensities within each subgroup were used in a Kolmogorov-Smirnov test to determine if the p-value distribution within the subgroup differed from that of the whole dataset. Log2 fold changes calculated previously were used as effect sizes in ChemRICH to indicate the magnitude of change within each subgroup. Putative root exudates were determined by the shared features between the rhizosphere and root metabolome. For features with no annotated name, we considered the difference of m/z within 0.001 and difference of retention time smaller than 0.5 the same features. The m/z and retention time tolerance interval was estimated from the differences between those shared features and confirmed ID.

To ensure reproducibility, our study employed a robust sampling design with 36 total samples (18 root and 18 rhizosphere) collected from three individual plants from four monoculture plots around three mesquite trees, utilizing a nested approach to control spatial variability. All samples were processed following rigorous protocols to maintain metabolomic integrity to the best of our ability. Chemical reagents were obtained from reputable vendors meeting high-purity standards, with all sources, instrument configurations, and software versions meticulously documented to facilitate potential future replication. Data processing and statistical analyses followed standardized procedures, with all data processing and statistic analysis workflows thoroughly documented and openly accessible.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.