The Metabolic Response of Brachypodium Roots to the Interaction with Beneficial Bacteria Is Affected by the Plant Nutritional Status

The potential of plant growth promoting (PGP) bacteria in improving the performance of plants in suboptimal environments is increasingly acknowledged, but little information is available on the mechanisms underlying this interaction, particularly when plants are subjected to a combination of stresses. In this study, we investigated the effects of the inoculation with the PGP bacteria Azospirillum brasilense (Azospirillum) on the metabolism of the model cereal Brachypodium distachyon (Brachypodium) grown at low temperatures and supplied with insufficient phosphorus. Investigating polar metabolite and lipid fluctuations during early plant development, we found that the bacteria initially elicited a defense response in Brachypodium roots, while at later stages Azospirillum reduced the stress caused by phosphorus deficiency and improved root development of inoculated plants, particularly by stimulating the growth of branch roots. We propose that the interaction of the plant with Azospirillum was influenced by its nutritional status: bacteria were sensed as pathogens while plants were still phosphorus sufficient, but the interaction became increasingly beneficial for the plants as their phosphorus levels decreased. Our results provide new insights on the dynamics of the cereal-PGP bacteria interaction, and contribute to our understanding of the role of beneficial microorganisms in the growth of cereal crops in suboptimal environments.


Introduction
Low temperature and P deficiency are two of the most common abiotic stresses for crops, and their combination can be particularly detrimental during the early stages of the development of winter cereals in P-poor regions of the world [1].
Low temperature affects virtually all metabolic processes of plants by changing the enzyme kinetics of biochemical reactions, inhibiting photosynthesis, damaging intra-and intercellular structures, decreasing water and nutrient uptake capacity, and causing oxidative stress [2,3]. Plants use P in most of their metabolic reactions, and phosphate deficiency activates various metabolic rearrangements, which increase external inorganic P (Pi) uptake and optimize internal Pi use efficiency [4]. In the case of extended P deprivation, plants often start recycling P from P-containing biomolecules such as RNA, phospholipids and small phosphorylated metabolites, and the metabolic pathways these compounds are involved in are severely impacted by prolonged P deficiency [5].
Low temperatures can affect plant metabolic responses to low P availability and vice versa. By decreasing plant growth rate, low temperature might indirectly affect plant P requirements; furthermore, these stresses have an opposite effect on root developmentwhich is generally stimulated by P deficiency and hindered by low temperatures [6]. Low temperature hampers photosynthesis by inhibiting sucrose synthesis, and this causes an accumulation of phosphorylated intermediates and an imbalance in the ATP/ADP ratio [7]. As mentioned before, these compounds are also affected by low P availability, as they decrease in P deficient plants. Arabidopsis thaliana mutants with reduced shoot Pi adapt to low temperature better than mutants with increased shoot Pi, but this might be reversed in highly P deficient plants, as these plants might not be able to readjust their metabolism to face low temperature [8]. When comparing the metabolic response of maize leaves to low temperatures and P limitation, Schlüter et al. [9] reported that plants subjected to low temperatures increased their carbohydrate and phosphorylated metabolite levels, while the opposite was observed in P deficient plants.
Plant growth promoting (PGP) bacteria can improve plant adaptation to abiotic stresses by affecting their metabolism, particularly by producing PGP compounds, improving nutrient uptake and interfering with plant stress signaling pathways [10][11][12]. The metabolic response of cereal roots to inoculation with PGP bacteria has been increasingly investigated in recent years [13][14][15][16][17], and few of those studies have analyzed the role of PGP bacteria in shaping the root metabolome of cereals facing stress [18,19]. Gagné-Bourque et al. [18] studied the impact of PGP bacteria at various stages of the interaction with timothy grass (Phleum pratense) subjected to drought. Their results are of particular interest as inoculation affected plant phenotype and metabolome differently through the experiment, showing the plasticity of the plant-bacteria interaction in different conditions. There is instead a lack of time-resolved studies on the effects of the interaction with PGP bacteria on the metabolome of cereal roots subjected to low P and cold temperatures, applied separately or simultaneously.
This study is part of a project that investigated the effect of the interaction with the PGP bacteria Azospirillum brasilense (Azospirillum) on the development of the model cereal Brachypodium distachyon (Brachypodium) subjected to suboptimal temperature and low P availability. In our previous phenotypic analysis, we compared the development of plants grown at suboptimal temperatures and various P levels and inoculated or not with Azospirillum. At the lowest P availability, bacteria increased shoot biomass in the early stages of the interaction, while the shoot nutrient concentration was not affected. Additionally, inoculated plants featured an increase of their branch root length in the second half of the experiment, concurrently with plant P deficiency. We hypothesized that the inoculation improved plant resistance to P deficiency stress by optimizing their root architecture for the exploration of the substrate [1]. Thus the current study focuses on investigating the metabolism of the previously phenotyped Brachypodium roots, with plants grown in the conditions where Azospirillum had a positive effect on their development [1]. We resolved the changes occurring in the root polar metabolites and lipids at various stages of the interaction with Azospirillum, and linked them to the plant nutritional status and phenotype. The levels of phosphorylated compounds decreased throughout the experiment, causing relevant rearrangements in plant metabolism. Inoculated plants synthesized defense-related compounds during the early stage of the interaction, while at later stages inoculation decreased abiotic stress. We propose that the plant P deficiency shaped the interplay with the bacteria, which in turn affected the plant response to the abiotic stresses they were subjected to.

Results
Brachypodium distachyon Bd21-3 (Brachypodium) [20] with or without Azospirillum brasilense Sp245 (Azospirillum) [21] inoculation were grown in a hydroponic system for 21 d at suboptimal temperatures (20/10 • C day/night), and throughout the experiment plants were supplied with a low P (7 µM KH 2 PO 4 ) nutrient solution. At 7, 14, and 21 days after the inoculation (DAI), roots were harvested and polar metabolites and lipids were extracted and analyzed via untargeted GC-MS or LC-MS. The root metabolic profiles were compared between inoculated and non-inoculated plants at each time point, and between different time points of the same treatment.
The analysis of the plant phenotype at various time points in the interaction with Azospirillum has been reported by Schillaci et al. [1].

P-Containing Compounds
Five P-containing compounds were detected in the root samples: glycerol-3-phosphate (G3P), inositol-1-phosphate, mannose-6-phosphate, myo-inositol-2-phosphate, and phosphoric acid. The levels of all these compounds decreased between 7 and 21 DAI, with a similar magnitude in both treatments ( Figure 1). The only compound significantly affected by the inoculation was G3P, which was more abundant (fold change (FC) = 1.53) in Azospirillum-inoculated plants at 7 DAI (Figure 1a).

Results
Brachypodium distachyon Bd21-3 (Brachypodium) [20] with or without Azospirillum brasilense Sp245 (Azospirillum) [21] inoculation were grown in a hydroponic system for 21 d at suboptimal temperatures (20/10 °C day/night), and throughout the experiment plants were supplied with a low P (7 µM KH2PO4) nutrient solution. At 7, 14, and 21 days after the inoculation (DAI), roots were harvested and polar metabolites and lipids were extracted and analyzed via untargeted GC-MS or LC-MS. The root metabolic profiles were compared between inoculated and non-inoculated plants at each time point, and between different time points of the same treatment.
The analysis of the plant phenotype at various time points in the interaction with Azospirillum has been reported by Schillaci et al. [1].

P-Containing Compounds
Five P-containing compounds were detected in the root samples: glycerol-3phosphate (G3P), inositol-1-phosphate, mannose-6-phosphate, myo-inositol-2-phosphate, and phosphoric acid. The levels of all these compounds decreased between 7 and 21 DAI, with a similar magnitude in both treatments ( Figure 1). The only compound significantly affected by the inoculation was G3P, which was more abundant (fold change (FC) = 1.53) in Azospirillum-inoculated plants at 7 DAI (Figure 1a). Relative response of P-containing compounds in roots of plants with or without Azospirillum inoculation and harvested after 7, 14, and 21 DAI. Means ± standard error are presented. n = 6 for all samples except for non-inoculated at 21 DAI (n = 8). Asterisks in the graphs indicate a fold change (FC) <−1.5 or >1.5 and probability of significant difference between relative response of Azospirillum-inoculated and non-inoculated plants based on Student's t-test p < 0.05 (see also Table S1).

Figure 1.
Relative response of P-containing compounds in roots of plants with or without Azospirillum inoculation and harvested after 7, 14, and 21 DAI. Means ± standard error are presented. n = 6 for all samples except for non-inoculated at 21 DAI (n = 8). Asterisks in the graphs indicate a fold change (FC) <−1.5 or >1.5 and probability of significant difference between relative response of Azospirillum-inoculated and non-inoculated plants based on Student's t-test p < 0.05 (see also  Table S1).

Comparison between Polar Metabolites of Inoculated and Non-Inoculated Plants
The response of polar compounds detected in Brachypodium roots was assessed using relative quantification and compared between inoculated and non-inoculated samples at each time point. Hierarchical clustering paired with a heatmap of metabolite abundances and principal component analysis (PCA) allowed the behavior of root metabolic profiles to be visualized throughout the experiment.
Polar metabolite patterns firstly clustered based on time point, while the bacterial inoculation was a weaker discriminant ( Figure 2, Figure S1). Based on the levels of detected metabolites, samples from different treatments were often more similar than samples belonging to the same treatment. The clearest separation was observed at 14 DAI, when samples from the two treatments clustered separately ( Figure 2). The response of polar compounds detected in Brachypodium roots was assessed using relative quantification and compared between inoculated and non-inoculated samples at each time point. Hierarchical clustering paired with a heatmap of metabolite abundances and principal component analysis (PCA) allowed the behavior of root metabolic profiles to be visualized throughout the experiment.
Polar metabolite patterns firstly clustered based on time point, while the bacterial inoculation was a weaker discriminant ( Figure 2, Figure S1). Based on the levels of detected metabolites, samples from different treatments were often more similar than samples belonging to the same treatment. The clearest separation was observed at 14 DAI, when samples from the two treatments clustered separately ( Figure 2).  While most metabolites did not show significant fold changes between the two treatments at any time point, the response of a small group of metabolites varied dynamically throughout the time series (FC < −1.5 or FC > 1.5, p < 0.05, Table S1, Figure S2). A total of 18 compounds were significantly increased or decreased when the two treatments were compared at each time point. These included aminocaproic acid, campesterol, βcyano-alanine, diethylene glycol, G3P, hexacosanol, malic acid, α-ketoglutaric acid (α-KG), pantothenic acid, pentonic acid, 1,4-lactone, putrescine, pyroglutamic acid, quinic acid, ribonic acid, serotonin, sinapic acid, trehalose, and xylose.
Apart from pentonic acid, 1,4-lactone, trehalose and aminocaproic acid, no compounds were significantly affected at more than one time point.

Root Glycerophospholipids and Galactolipids
The comparison of specific lipid classes at different time points showed converse trends between GPs (PCs and PEs) and galactolipids (MGDGs and DGDGs) in both treatments (Table S2b)

Comparison between Lipids of Inoculated and Non-Inoculated Plants
Lipid and unidentified feature profiles of inoculated and non-inoculated plants were compared at 7, 14, and 21 DAI. The hierarchical clustering and PCA plot of lipid species and unidentified features detected at different time points show that features separated based on time of harvest. Lipid species did not separate based on the bacterial treatment ( Figure 3, Figure S3), while a more pronounced clustering was observed in unidentified features, particularly at 7 and 21 DAI ( Figure S4).

Comparison between Lipids of Inoculated and Non-Inoculated Plants
Lipid and unidentified feature profiles of inoculated and non-inoculated plants were compared at 7, 14, and 21 DAI. The hierarchical clustering and PCA plot of lipid species and unidentified features detected at different time points show that features separated based on time of harvest. Lipid species did not separate based on the bacterial treatment ( Figure 3, Figure S3), while a more pronounced clustering was observed in unidentified features, particularly at 7 and 21 DAI ( Figure S4).  To further investigate the effects of bacterial inoculation on the Brachypodium root lipidome, K-means clustering analysis was performed separately on the lipids identified in inoculated and non-inoculated samples, which were grouped based on their levels at 7, 14, and 21 DAI (Figure 4). Seven lipid clusters were fitted and grouped between the two treatments based on their metabolite composition. The degree of similarity between the clusters in each pair was generally good, varying between 48.7 and 75.5% for six out of the seven paired clusters (for the complete cluster composition and the similarity scores of each cluster pair, see Table S2c). The clusters of inoculated and non-inoculated samples differed particularly in the third and fourth pairs. In the third pair (Figure 4c), DGs and GPs identified in non-inoculated samples were generally stable between 7 and 14 DAI, and decreased at 21 DAI, while those from inoculated samples increased between 7 and 14 DAI and decreased again between 14 and 21 DAI. In the fourth pair (Figure 4d), TGs from non-inoculated samples decreased between 7 and 14 DAI and were generally stable between 14 and 21 DAI, while features from inoculated samples were stable between 7 and 14 DAI and increased between 14 and 21 DAI. The analysis of lipid species in the two treatments at different time points shows a limited effect of the inoculation on individual features. DGDG 36:1, LPC 14:0, and PE 33:2 were less abundant in inoculated plants at 7 DAI; PC 30:1 was more abundant in inoculated plants at 14 DAI. No lipid species differed significantly at 21 DAI ( Figure S5).
To further investigate the effects of bacterial inoculation on the Brachypodium root lipidome, K-means clustering analysis was performed separately on the lipids identified in inoculated and non-inoculated samples, which were grouped based on their levels at 7, 14, and 21 DAI (Figure 4). Seven lipid clusters were fitted and grouped between the two treatments based on their metabolite composition. The degree of similarity between the clusters in each pair was generally good, varying between 48.7 and 75.5% for six out of the seven paired clusters (for the complete cluster composition and the similarity scores of each cluster pair, see Table S2c). The clusters of inoculated and non-inoculated samples differed particularly in the third and fourth pairs. In the third pair (Figure 4c), DGs and GPs identified in non-inoculated samples were generally stable between 7 and 14 DAI, and decreased at 21 DAI, while those from inoculated samples increased between 7 and 14 DAI and decreased again between 14 and 21 DAI. In the fourth pair (Figure 4d), TGs from non-inoculated samples decreased between 7 and 14 DAI and were generally stable between 14 and 21 DAI, while features from inoculated samples were stable between 7 and 14 DAI and increased between 14 and 21 DAI.

Discussion
The comparison of polar metabolites and lipids extracted from the roots of Azospirilluminoculated and non-inoculated roots of Brachypodium at 7, 14, and 21 DAI revealed that the time of the harvest was a strong discriminant between samples, while the bacterial treatment affected specific sections of root metabolism at different time points.

Plants Consumed Phosphorylated Compounds to Face the Increasing P Deficiency
All detected P-containing compounds decreased in both treatments through the experiment ( Figure 1) as might be expected given the progressive decrease of shoot P concentration in those same plants [1]. The reduction in phosphoric acid indicates that the inorganic P reserves in the vacuole were declining (Figure 1e), a process associated with the reduction in phosphorylated metabolites as internal P recycling from senescing tissues occurred [22]. Plants progressively switched to metabolic pathways that require little or no phosphorylated compounds, allowing P to be remobilized for essential metabolism [5]. The large reduction in all detected P-containing compounds suggests that not only were the shoots of both treatments increasingly P deficient [1] but this was also the case in the roots. As previously reported, despite a similar metabolic response indicating P deficiency the growth of shoots and roots was stronger in plants inoculated with Azospirillum.

Azospirillum Elicited Different Responses in Plants at Different Stages
The profile of polar metabolites extracted from Brachypodium roots with and without Azospirillum inoculation changed through time, suggesting that bacterial interactions dynamically and specifically altered plant metabolism at low temperature and P conditions. The time point at which inoculation induced the largest number of metabolite changes was 14 DAI. Eleven metabolites had significant changes, of which 4 were increased and 7 were decreased by inoculation. Changes in polar metabolites were more limited at 7 and 21 DAI, with four compounds increased and one decreased in inoculated plants at both time points (Table S1, Figure S2).
The metabolic pathways associated with significant changes in this study are shown in Figure 5. Compounds significantly affected by the inoculation are scattered across primary metabolism rather than belonging to specific pathways, with the partial exception of sugar metabolism (trehalose, xylose) and the citric acid cycle (malic acid, α-ketoglutaric acid(α-KG)). It was not possible to map hexacosanol, diethylene glycol, pentonic acid, 1,4-lactone, and ribonic acid onto Figure 5, because a pathway for those compounds has not been characterized in plants yet. Further information on affected compounds in relation to stresses and/or interaction with PGP bacteria can be found in Appendix A Table A1.
Dynamic metabolite fluctuations in association with bacterial inoculation are evident by mapping affected compounds at different time points ( Figure 5). This reveals metabolites which were most significantly changed in response to growing conditions (time, inoculation).
We cannot exclude that compounds increased in inoculated plant roots were produced by the bacteria or by both bacteria and the plant but, due to the great disproportion between host plant and bacterial biomass, it is unlikely that bacterial metabolites had a significant impact on the observed results.  It is noteworthy that all compounds whose levels were increased in inoculated plants at 7 DAI-hexacosanol, campesterol, trehalose, and G3P-can be involved in plant response to biotic stresses and this suggests that the plant was initially producing antimicrobial compounds to limit bacterial colonization of the roots. Hexacosanol stimulates wax barriers formation, and can also be one component of the wax layer suberin, a structure involved in plant response against biotic and abiotic stresses [23,24]. Campesterol is a precursor for the synthesis of brassinosteroids, a class of hormones involved in plant responses to environmental stresses, both biotic and abiotic [25][26][27]. Brusamarello-Santos et al. [14] report that trehalose levels increased in maize inoculated with Azospirillum brasilense or Herbaspirillum seropedicae, and the authors hypothesize a role of this compound in the regulation of the plant-bacteria interaction, possibly in defense related mechanisms. G3P plays an essential role as a signaling molecule in plant systemic acquired resistance, a defensive mechanism against a broad variety of pathogens which is activated during early infection stages, preventing pathogen propagation in distal tissues [28,29]. This type of defense response is not unusual; PGP bacteria can be sensed as pathogens during early stages of the interaction, as both beneficial and detrimental microorganisms can trigger the plant immune system [17]. Once different recognition mechanisms of plant and microorganisms are established, the plant response to the PGP bacteria changes as was seen at 14 and 21 DAI.
Brassinosteroids, the phytohormones synthesized from campesterol, are normally involved in cell elongation [30], and specific concentrations can stimulate branch root formation [31]. The untargeted analysis we performed only allowed us to compare the levels of detected compounds and did not provide information on their absolute quantity, but it is possible that the higher campesterol levels in inoculated Brachypodium at 7 DAI ( Figure 5, Table S1) caused the earlier development of branch roots observed in those plants [1]. As branch roots are of great importance for plants growing in a P-poor environment [32] this could have produced a significant advantage for Azospirillum-treated plants.
Trehalose and sucrose are sugars separated by two enzymatic reactions only ( Figure 5) and can accumulate in plants subjected to cold stress [33]. Compared to sucrose, trehalose biosynthesis requires much smaller amounts of carbon, and allows the uridine diphosphate glucose and glycerol-6-phosphate pools to be conserved [34]. In stressed plants, the metabolism of sucrose and trehalose is typically related [33], but in our study the two compounds had analogous trends only in non-inoculated plants, with a decrease during the second week of the experiment and an increase during the third. In contrast, the trehalose content was relatively stable in inoculated plants ( Figure 5), while the sucrose content followed the same trend as in non-inoculated plants ( Figure S6). This suggests that the bacteria caused a shift in sugar metabolism of inoculated plants, redirecting it towards less P-consuming pathways. Finally, genes involved in trehalose biosynthesis differed in the roots of two rice cultivars subjected to various stresses, including P deficiency [35]. Plants with higher expression of such genes had shorter primary root and better developed shoots when subjected to P deprivation, which is consistent with the phenotype of Brachypodium plants inoculated with Azospirillum [1]. We hypothesize a role for trehalose in shaping the development in cereals facing P deficiency, in particular preventing the redirection of resources from the shoot to the root growth, which is a common response of plants to this stress [36].
Compared to samples harvested at 7 DAI, those from 14 DAI show a higher number of affected metabolites ( Figure 5, Table S1). The characteristics of most compounds whose levels were changed at 14 DAI-quinic acid [37], α-KG [38], malic acid [5,39,40], β-cyanoalanine [41,42], and sinapic acid [43]-suggest inoculation with Azospirillum provided protection against both nutrient deficiency and low temperature stress (see also Appendix A Table A1). This hypothesis is supported by the reduction in compounds with antimicrobial activity compared 7 DAI in inoculated plants, as only sinapic acid is reported to be involved in plant defense mechanisms against pathogens [44].
Since plants undergoing P deficiency exudate organic acids to increase substrate P availability [39], they need to use amino acids as alternative carbon sources for various metabolic processes. One of the by-products of amino acid degradation is NH 4 + , which can reach toxic concentrations in the plant. To reduce NH 4 + levels, plants can use amination of the organic acid α-KG to synthesize specific amino acids and amines (asparagine, glutamine, putrescine). This leads to reduced abundance of α-KG in plants facing severe P starvation [5] as was seen in non-inoculated plants between 7 and 14 DAI in our study ( Figure 5). In inoculated plants the abundance of α-KG increased at 14 DAI. As the production of this organic acid does not seem to be part of the A. brasilense strategy to increase P availability [45], this result suggests that the bacterial interaction with Brachypodium preserved the plant TCA cycle, preventing the use of α-KG in reducing ammonium levels at this stage. This hypothesis is also supported by the fact that at 14 DAI, inoculated plants had lower levels of putrescine (Table S1, Figure 5), one of the compounds used by plants to incorporate NH 4 + . The higher levels of β-cyano-alanine in non-inoculated Brachypodium roots suggests that they were facing higher levels of stress compared to inoculated plants. β-cyano-alanine is a product of the degradation of cyanide which, in turn, is the biproduct of the conversion of 1-amino-cyclopropane-1-carboxylic acid (ACC) to ethylene, a hormone whose synthesis is enhanced in plants facing stress [42]. Thus, stressed plants can often produce high levels of cyanide, which needs to be quickly metabolized as it inhibits electron transport and metalloenzymes [46]. The main pathway for this process is the β-cyano-alanine synthase pathway, which incorporates the cyanide into a cysteine molecule forming β-cyano-alanine, subsequently converted into asparagine, aspartate, and ammonia, and these compounds can be used for plant primary metabolism [47]. In this study, inoculated Brachypodium plants had significantly lower β-cyano-alanine content at 14 DAI (Table S1, Figure 5), which suggests a lower production of ethylene at this specific time point and hence lower stress levels. Non-inoculated plants had higher levels of asparagine (FC = 2.63) compared to the inoculated ones. Unfortunately, asparagine quantification is often problematic using GC-MS and LC-MS is preferred [48]. The detected levels in our samples varied substantially within each treatment resulting in asparagine FC at 14 DAI being not statistically significant (p-value = 0.68, see Table S1). However, this result at least partially supports the hypothesis of cyanide detoxification occurring in non-inoculated plants, as asparagine is one of the products of the β-cyano-alanine synthase pathway ( Figure 5).
A smaller number of metabolites differed significantly between the two treatments at 21 DAI, and the identity of most compounds changed compared to the previous time point (Table S1, Figure 5). This suggests a change in the plants' strategy to face the increasing P deficiency in their tissues. As at 14 DAI, some affected compounds-ribonic acid [49] and serotonin [50]-may play a role in plant adaptation to suboptimal conditions, rather than limiting bacterial growth as observed at 7 DAI.
Ribonic acid levels were generally steady through time in non-inoculated plants, while there was a great increase in inoculated plants between 14 and 21 DAI (Figure 5). In a recent study, the foliar application of γ-aminobutyric acid (GABA) increased the leaf concentration of ribonic acid in the leaves of the grass Agrostis stolonifera subjected to drought and improved their adaptation to this stress [51]. Azospirillum has been reported to increase the GABA concentration in the tissues of inoculated sorghum plants [52], and similar results were found in our study, where the levels of GABA increased in inoculated plants between 14 and 21 DAI ( Figure S7). This suggests a connection between ribonic acid and GABA in Brachypodium's responses to abiotic stresses.
Serotonin, also known as 5-hydroxytryptamine, is an indoleamine mainly synthesized in plant roots [53] and various studies suggest that serotonin can play a primary role in plant resistance against biotic [54] and abiotic [50] stresses. Most importantly, serotonin can shape root architecture of plants, modulating the development of different root types. Arabidopsis primary and lateral root growth is inhibited by high exogenous levels of serotonin, which affects cell division and elongation, particularly of the stem cell niche. Conversely, lower serotonin levels, while not affecting primary root, increased root branching by stimulating the maturation of lateral root primordia [55,56]. The difference in serotonin levels of the two treatments at 21 DAI (Table S1, Figure 5) could suggest a role of serotonin in shaping the root architecture of Brachypodium, which differed in inoculated plants particularly in the last week of the experiment [1]. Our study did not differentiate between the metabolomes of primary and branch roots, and future studies should determine whether serotonin concentrations differ between inoculated root types, possibly inhibiting primary root growth while at the same time stimulating branch roots.

Lipids
Untargeted analysis of lipid compounds extracted from Azospirillum-inoculated and non-inoculated Brachypodium roots at various time points of their growth provided three major insights. First, lipids changed in roots grown under both conditions similarly through time, thus indicating a developmental effect on root lipid composition. Secondly, both treatments reacted in a similar way to the increasing P deficiency, with a strong increase in galactolipids and a strong decrease in GP levels. Finally, the bacterial treatment did not significantly affect the root lipid profile, at any of the time points; however, specific lipid classes behaved differently through time in the two treatments.

P Deficiency Strongly Remodeled Root Lipid Profiles of Both Treatments
The levels of both GPs and galactolipids fluctuated similarly in both treatments, with most GPs decreasing and most galactolipids increasing between 7 and 21 DAI (Table S2b). As mentioned previously, plants subjected to P deficiency can react by recycling internal P for essential metabolic processes. GPs are second only to RNA as the most abundant phosphorylated compounds in plants, and their replacement with non-phosphorous lipids in plasma membranes, mitochondrial membranes, and tonoplasts is an established strategy in roots of various cereal crops to deal with P deficiency [57].
In particular, GPs can be converted to MGDGs and DGDGs, and the latter can then form a bilayer with similar characteristics to the GP bilayer in cell membranes [58,59]. When oat plants were supplied with varying P levels for 30 days, low P caused an increase in the DGDG/PC ratio in plant roots through time, while this ratio remained stable in plants supplied with optimal P levels [57]. As our study did not include a P-sufficient treatment, we cannot exclude that the decrease of GPs and concurrent increase of galactolipids in both the treatments through time was caused by the aging of the plant. Nevertheless these trends, together with the decreased shoot P concentration [1] and the reduction in Pcontaining polar compounds ( Figure 1) observed in both treatments, are further indications of the increasing P deficiency that plants were facing.

Azospirillum Inoculation Had a Limited Effect on Brachypodium's Root Lipids
The comparison of K-means clusters of identified features through time revealed that when treated as a single class of mean responses the GP, DG, and TG lipid classes behaved differently in the two treatments through time ( Figure 4).
Conversely, the comparison of root lipid profiles and of the abundances of single lipid species between the two treatments revealed a weak effect of the bacterial inoculation. No separation between the Azospirillum-inoculated and non-inoculated samples was observed at any time point (Figure 3) and lipids significantly affected by the inoculation during the experiment were limited: three lipid species were less abundant in inoculated plants at 7 DAI and one was more abundant in inoculated plants at 14 DAI ( Figure S5). Inoculation effects were more profound on unidentified features, which displayed a clearer clustering at different time points ( Figure S4).
Very few studies have so far focused on the effect of beneficial bacteria on the root lipid composition of plants, either grown in optimal conditions or facing stresses. Due to the importance of GPs in cell membranes, most studies focus on the effects of bacterial inoculation on this lipid class. The K-means clustering suggests that, through time, GPs, might have behaved differently in the two treatments ( Figure 4c). As DGs can be a product of GP hydrolysis, the fact that their levels changed similarly through time is noteworthy, but further studies on the importance of this lipid class in interactions between plants and PGP bacteria are required.
Zhang et al. [60] report that in the interaction between soybean and rhizobia, TG synthesis was decreased in nodules compared to roots. Conversely, the halophyte Bassia indica grown under salt stress showed improved growth and increased TG content in shoots when inoculated with the PGP bacteria Bacillus subtilis [61]. Calderon-Vazquez et al. [62] hypothesize that TG degradation could be part of the maize response to P deficiency, to provide acyl groups for the synthesis of glycolipi ds (e.g., galactolipids) as GPs substitutes. Overall, the available studies on the link between plant TGs, PGP bacteria, and P deficiency are still scarce and results are inconsistent, but the different behavior of this class of lipids between inoculated and non-inoculated plants in our study (Figure 4d) suggests that they may be involved in plant response to inoculation and could have been involved in plant adaptation to low P.
Although in the earlier stages of the experiment (7 DAI) a number of defense-related polar metabolites were affected, the same was not observed for the identified lipid species, despite membrane lipids being pivotal in plant interactions with microorganisms [63]. In general, the limited effect of inoculation on Brachypodium root lipid profile could have different explanations. Most of the cited studies on plant-microbe interactions showed lipid rearrangements in the first stages of the interaction (12-48 h) [64,65]. It is possible that, while the immune response involving polar metabolites continued at least until 7 DAI, the lipid response happened mainly in the first hours/days of the interaction. Since many of the membrane lipids involved in plant response to biotic interactions are GPs, it is also possible that Brachypodium response to Azospirillum was heavily affected by the increasing P deficiency they suffered through the experiment, and this would also explain why the number of lipids significantly affected by the inoculation gradually decreased during the experiment. The third possible explanation is that some affected lipids lie among the features that were not identified as lipids using the currently available libraries, as reliable and comprehensive databases for the identification of detected lipid species are still under construction [66]. Finally, the fact that the levels of identified lipid species did not vary between the two treatments does not necessarily mean that plants were not synthetizing more or less of certain lipid species and those changes were hidden by the bacterial production or consumption of those same compounds. Transcriptomic analysis of inoculated and non-inoculated plant roots could help clarify this hypothesis, by linking the abundance detected metabolites with the expression of the genes involved in their synthesis and/or degradation.
Differences between inoculated and non-inoculated plants were observed in the fluctuations of specific lipid classes by K-means cluster analysis. Compared to other unsupervised methods for the graphical representations of high-dimensional data such as PCA and hierarchical clustering, K-means clustering can be more effective in identifying subtle changes in metabolite abundance [67], and allowed the effect of inoculation on DGs, GPs and TGs to be detected but this will need to be further investigated in future studies.

Plant Growth Conditions
The experimental setup and conditions for plant inoculation and growth are described in detail in Schillaci et al. [1]. Briefly, 48-h-old Brachypodium distachyon Bd21-3 (Brachypodium) seedlings were inoculated by submerging them in a bacterial solution of Azospirillum brasilense Sp245, while sterile PBS was used as mock-inoculum for noninoculated seedlings. Plants were then grown for 21 d in the GrowScreen-PaGe platform [68], a hydroponic system where they were supplied weekly with a modified Hoagland solution containing extremely low (7 µM KH 2 PO 4 ) P levels. For the whole experiment, plants were grown at 20/10 • C day/night temperature. Plants were harvested at three time points: at 7 DAI (84 inoculated and 87 non-inoculated); at 14 DAI (24 inoculated and 36 non-inoculated) and at 21 DAI (25 inoculated and 34 non-inoculated). In order to minimize variations caused by circadian changes in the plant metabolism [69], all tissue harvesting was done at the same time of the day.
Roots were separated from shoots using a scalpel and immediately snap-frozen in liquid nitrogen. Frozen root tissues were ground to a fine powder using a Retsch Mixer Mill MM400 (Retsch GmbH, Haan, Germany) and plants from the same treatment and time point were pooled into samples with approximately the same number of individuals and fresh weight, measured using a digital scale with a 0.1 mg resolution (Mettler Toledo, Columbus, OH, USA). Six pools were formed for each treatment × time point, except for non-inoculated plants at 21 DAI which had eight groups. Samples were kept frozen during the grinding and pooling steps, to maintain quenched metabolism and thus prevent degradation of metabolites.

Polar Metabolites Extraction from Roots
For polar metabolites analysis, a variation of the method described in Cheong et al. [70] was used. 250 µL 100% LC-MS grade MeOH (Roth, Karlsruhe, Germany), containing 4% 13 C 6 sorbitol/valine (Sigma-Aldrich, Castle Hill, Australia) was added to approximately 25 mg of frozen root tissue from each pooled group. Tubes were shaken at 800 rpm for 15 min at 30 • C, centrifuged at 15,700× g for 15 min at room temperature and the supernatant collected. The pellet was re-extracted with 250 µL of Milli-Q H 2 O as above and the supernatants combined. In case of cloudy supernatant or precipitate observed during supernatant transfer, the pooled supernatant was centrifuged at 15,700× g for 10 min at room temperature before transferring as much supernatant as possible to a new tube. From each sample, 50 µL aliquots of the supernatant were transferred into glass insert (Agilent, Santa Clara, CA, USA) and dried under vacuum at 30 • C for 90 min.

Lipid Extraction from Roots
Lipids were extracted using a single step extraction protocol [71] as described by Kehelpannala et al. [72]. Briefly, 10 mg of frozen root tissue from each pooled group was homogenized with 400 µL of isopropyl alcohol (Roth, Karlsruhe, Germany) containing 25 µM deuterated cholesterol (internal standard, Sigma-Aldrich, Munich, Germany) and 0.01% butylated hydroxytoluene (BHT, Sigma-Aldrich, Munich, Germany). The tubes were then shaken at 1400 rpm for 15 min at 75 • C. Once cooled to room temperature, 1.2 mL of a mixture of chloroform: MeOH: water (30:41.5:3.5, v:v:v) mixture were added to each sample (all reagents were LC-MS grade and purchased from Roth, Karlsruhe, Germany). The tubes were then shaken at 300 rpm for 24 h at 25 • C and centrifuged at 15,700× g for 15 min at room temperature. The supernatant was then dried down with a N 2 gas stream for approximately 50 min.

GC-MS Analysis of Polar Metabolites
Polar metabolite extracts were derivatized as described by Dias et al. [48], and 1 µL of sample was injected onto the GC column using a hot needle technique. Each sample was injected pure (splitless) and 1:30 diluted (split), in order to detect both lowly and highly abundant metabolites within the dynamic range of the instrument.
The GC-MS system used for the analysis was composed of Gerstel 2.5.2 autosampler, a 7890A Agilent gas chromatograph and a 5975C Agilent quadrupole mass spectrometer (Agilent, Santa Clara, CA, USA).
Gas chromatography was performed following the method described by Hillyer et al. [73], with the only differences that samples were run through a 30 m Agilent J & W VF-5MS column with 0.25 µm film thickness and 0.25 mm internal diameter with a 10 m Integra guard column.
The chromatograms and mass spectra obtained from the split and splitless injections were analyzed using the Agilent MassHunter Qualitative Analysis software (v B.07.00, Agilent Technologies, Santa Clara, CA, USA) and the AMDIS software (v 2.73, (National Institute of Standards and Technology, Gaithersburg, MD, USA) for identification of the compounds, and the Agilent MassHunter Quantitative Analysis software (v B.08.00) for their quantification. The commercial mass spectra library NIST (http://www.nist.gov, accessed on 16 September 2020) and the in-house Metabolomics Australia mass spectral library were used as references for the compound identification. During the compound identification step, chromatograms with high peak intensity analysis for each treatment/time point set of samples were selected to build the list of detected compounds in the samples.
The resulting data were normalized to the 13 C 6 sorbitol/valine internal standard value, and to the exact weight of tissue used for metabolite extraction. When integrating the data from compounds detected in both the splitless and the split samples run, data from the splitless run were preferred, unless the chromatogram showed poor quality due to the overloading of the column or the overlap of different compounds.

LC-MS Analysis of Lipids
Dried lipid extracts were resuspended in 200 µL of a butanol: MeOH (1:1, v/v) mixture containing 10 mM ammonium formate. Eight pooled biological quality control (PBQCs) samples were prepared by pooling 10 µL from each resuspended sample, and aliquoting this sample into eight vials.
Untargeted analysis of samples extracted from Brachypodium roots was then performed using the liquid chromatographic conditions derived from those described by Kehelpannala et al. [72], but using an injection volume of 10 µL. Lipids were analyzed using a Sciex Triple TOF TM 6600 QqTOF mass spectrometer (Sciex, Framingham, MA, USA) equipped with a Turbo VTM dual-ion source (electro-spray ionization (ESI) and atmospheric pressure chemical ionization (APCI)) and an automated calibrant delivery system (CDS) using Information Dependent Acquisition method (IDA) in positive ion mode. The parameters were set as follows: MS1 mass range, 80- Three root PBQC extract samples were further analyzed using a different IDA method and a Sequential Window Acquisition of All Theoretical Fragment-ion Spectra (SWATH) method, both in positive ion mode. For the SWATH analysis, the HPLC settings were the same as previously described, except that 15 µL of sample was injected in the column. The same mass spectrometer employed previously was set as described by Kehelpannala et al. [74]. The compiled list of lipids was used to annotate the mass spectrometric features extracted from the IDA of the root samples. The same settings described previously were used, except for the alignment parameter settings, where a retention time tolerance of 1 min and MS1 tolerance of 0.02 Da were used.
MS-DIAL output data consisting of peak areas of annotated lipids and unidentified features was preprocessed prior to comparing them among the different treatments and time points. Peak area of all detected features were normalized to the fresh weight of each sample, and features with coefficient of variation CV = SD(PBQCs) mean(PBQCs) above 20% were discarded, to ensure the reproducibility of the results [76].

Statistical Analyses and Data Visualization
The abundances of polar metabolites and lipids were compared using the web-based tool Metaboanalyst 4.0 (https://www.metaboanalyst.ca/, accessed on 13 November 2020). Lipid data were normalized by median, and both polar metabolite and lipid data were log transformed and auto scaled.
The metabolite profiles of the two treatments at each time point were displayed using 2D principal component analysis (PCA) and hierarchical cluster analysis of the root metabolome paired with a heatmap of metabolite abundances, with samples and features clustered using Euclidean distance function and Ward clustering algorithm.
The log transformed and auto scaled level of the compounds detected in the two treatments was compared at each time point, and the statistical significances between the measured abundances were analyzed using the Student's t-test (group variance assessed using the Levene's test). Compounds with a fold change (FC) < −1.5 or >1.5 and with raw p-value < 0.05 were considered differently abundant between inoculated and noninoculated treatments.
K-means clustering of the relative responses from identified lipid species through time was performed in RStudio v. 1.2.5019 (Rstudio, Boston, MA, USA) [77]. A function named "KmeansPlus" was written and compiled in the publicly available R package RandoDiStats. The built function uses and depends on the R package ComplexHeatmap [78] as is detailed in the GitHub repository (https://github.com/MSeidelFed/RandodiStats_package, accessed on 22 December 2020). The function call defaults were used. Briefly, Pearson correlation was taken as a similarity distance between lipid variables, average was taken as the cluster construction method and a clustering consensus from 1000 runs was compiled to prevent misrepresentation from an outlier run. Finally, seven partitions or clusters were formed based on the encountered patterns between conditions. The R function returned in the environment a data frame containing the lipid variables that belong to each of the produced clusters, easing the interpretation of results. Additionally, in the working directory, plots dependent on ggplot2 [79] were generated. The plots featured the mean intensity of lipids that belong to each cluster as solid-colored lines and the standard error as shade around the mean.

Conclusions
In summary, the comparison between the root metabolic profiles of inoculated and non-inoculated plants at different stages of their growth tells us that the interaction with Azospirillum affected Brachypodium metabolism dynamically through time, particularly with regard to polar metabolites.
During the early stages of the experiment, we hypothesize that bacteria were sensed as pathogens by the host, and Brachypodium produced polar compounds with antimicrobial properties in order to limit the bacterial colonization of roots. This trend changed at later stages of the experiment, when fewer and different defense compounds were produced, and most affected compounds suggest an improved response to P deficiency by inoculated plants. This hypothesis is supported by their enhanced root development compared to non-inoculated plants during the second half of the experiment [1].
We cannot exclude that some of the compounds detected in inoculated plants might be synthesized by Azospirillum when interacting with Brachypodium, and further research is required to determine with certainty which organism produced those metabolites. In order to resolve the mechanisms for promotion of Brachypodium plant growth by Azospirillum when plants were subjected to simultaneous abiotic stresses, our experimental setup was highly controlled and simplified. However this setup allowed detailed phenotypic analysis that would have been almost impossible in a soil based system. This study lays the foundation for future studies where cereals and bacteria interact in conditions closer to an agricultural setting, such as soil substrates characterized by fluctuating humidity and temperatures, and by the native microbiota. Further studies of the interaction between Brachypodium and Azospirillum in hydroponic systems will allow the accurate characterization of the root exudates, an essential component of plant-bacteria association [80,81].
In conclusion, we suggest that the increasing P deficiency changed the interaction between Brachypodium and Azospirillum, as the beneficial effects of bacterial inoculation were most evident at 14 and 21 DAI, and at those time points inoculated roots did not display higher content of antimicrobial compounds. Numerous studies which subjected inoculated plants to varying degree of stresses report that the beneficial effect of PGP bacteria can increase as the growing conditions worsen, thanks to their capacity of decreasing the damage in plants [82][83][84]. As soil degradation and climate change are likely to aggravate the growing conditions of crops in many areas of the world, our study confirms that PGP bacteria can be a valuable resource to still achieve profitable and sustainable agriculture in the future.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/metabo11060358/s1, Figure S1: Principal component analysis of the polar metabolite profiles of Azospirillum-inoculated (B) and non-inoculated roots (C) of Brachypodium harvested at 7, 14, and 21 DAI. Ellipses around samples display 95% confidence areas, Figure S2: Volcano plots of significantly affected polar metabolites in the comparison between bacteria inoculated and noninoculated samples at 7 DAI (a), 14 DAI (b) and 21 DAI (c). Red dots represent compounds with FC <−1.5 or >1.5 and probability of significant difference between mean of Azospirillum-inoculated and non-inoculated plants based on Student's t-test p < 0.05, Figure S3: Principal component analysis of the lipid profiles of Azospirillum-inoculated (B) and non-inoculated roots (C) of Brachypodium harvested at 7, 14, and 21 DAI. Ellipses around samples display 95% confidence areas, Figure S4: Hierarchical clustering coupled with heatmap of the unidentified feature profiles of Azospirilluminoculated (B) and non-inoculated roots (C) of Brachypodium harvested at 7, 14, and 21 DAI. The lettering at the bottom of the heatmap indicates the replicates: for each replicate, blue and red colors indicate the lower or higher abundance of specific lipids compared to the other replicates, with darker colors indicating more pronounced differences, Figure S5: Volcano plots of significantly affected lipids in the comparison between bacteria inoculated and non-inoculated samples at 7 DAI (a), 14 DAI (b), and 21 DAI (c). Red dots represent compounds with FC < −1.5 or >1.5 and probability of significant difference between mean of Azospirillum-inoculated and non-inoculated plants based on Student's t-test p < 0.05, Figure S6: Relative response of Sucrose in roots of plants with or without Azospirillum inoculation harvested at 7, 14, and 21 DAI. Means are presented. n = 6 for all pooled samples except for non-inoculated at 21 DAI (n = 8), Figure S7: Relative response of γ-aminobutyric acid in roots of plants with or without Azospirillum inoculation harvested at 7, 14, and 21 DAI. Means are presented. n = 6 for all pooled samples except for non-inoculated at 21 DAI (n = 8), Figure S8: Relative response of Pi-Cer in roots of plants with or without Azospirillum inoculation harvested at 7, 14, and 21 DAI. Means are presented. n = 6 for all pooled samples except for non-inoculated at 21 DAI (n = 8), Table S1: List of polar metabolites detected in Brachypodium roots harvested at 7, 14, and 21 DAI. Columns display in ascending order the class of each compound, its name, the relative response fold changes (FC) (Azospirillum inoculated/non-inoculated plants) and the related p-value for each time point. Compounds are divided into amino acids (AA), organic acids (OA) sugars (SU) or "Others" if they did not belong to any of the previous classes. Compounds with a FC < −1.5 or FC > 1.5 and p < 0.05 are presented in red font, Table S2: (a) Features detected with LC-MS analysis of Azospirillum-inoculated (B) and non-inoculated (C) Brachypodium roots harvested at 7, 14, and 21 DAI. Average retention time (Rt) and mass to charge ratio (m/z) of each compound are presented. Columns D to AO report the relative level of detected features, with darker colors indicating higher abundance compared to the other replicates. ADGGA: acyl diacylglyceryl glucuronide, Cer-AP: ceramide alpha-hydroxy fatty acid-phytospingosine, DG: diacylglycerol, DGDG: digalactosyldiacylglycerol, DGGA: diacylglyceryl glucuronide, HexCer: hexosylceramide, HexCer-AP: hexosylceramide alpha-hydroxy fatty acid-phytospingosine, LPC: lysophophatidylcholine, MGDG: monogalactosyldiacylglycerol, PC: phophatidylcholine, PE: phosphatidylethanolamine, PG: phosphatidylglycerol, PI: phosphatidylinositol, Pi-Cer: ceramide phosphoinositol, TG: triacylglycerol. (b) List of glycerophospholipids and galactolipids detected in Brachypodium roots harvested at 7 and 21 DAI. Columns display the relative response fold changes (FC) (21 DAI/7 DAI) and the related p-value for each treatment. Red and blue font indicate a FC < −1.5 or >1.5 respectively, and probability of significant difference between time points based on Student's t-test p < 0.05. DGDG: digalactosyldiacylglycerol, LPC: lysophosphatidylcholine, MGDG: monogalactosyldiacylglycerol, PC: phosphatidylcholine, PE: phosphatidylethanolamine, PG: phosphatidylglycerol, PI: phosphatidylinositol, Pi-Cer: ceramide phosphoinositol, (c) Composition and similarity of the K-means clusters of lipid species extracted from Azospirillum-inoculated and non-inoculated Brachypodium roots harvested at 7, 14, and 21 DAI. ADGGA: acyl diacylglyceryl glucuronide, Cer-AP: ceramide alphahydroxy fatty acid-phytospingosine, DG: diacylglycerol, DGDG: digalactosyldiacylglycerol, DGGA: diacylglyceryl glucuronide, HexCer: hexosylceramide, HexCer-AP: hexosylceramide alpha-hydroxy fatty acid-phytospingosine, LPC: lysophophatidylcholine, MGDG: monogalactosyldiacylglycerol, PC: phophatidylcholine, PE: phosphatidylethanolamine, PG: phosphatidylglycerol, PI: phosphatidylinositol, Pi-Cer: ceramide phosphoinositol, TG: triacylglycerol.     [108] no unknown α-ketoglutaric acid P deficiency rice [109] no unknown