Targeted and Untargeted Metabolomic Analyses Reveal Organ Specificity of Specialized Metabolites in the Model Grass Brachypodium distachyon

Brachypodium distachyon, because of its fully sequenced genome, is frequently used as a model grass species. However, its metabolome, which constitutes an indispensable element of complex biological systems, remains poorly characterized. In this study, we conducted comprehensive, liquid chromatography-mass spectrometry (LC-MS)-based metabolomic examination of roots, leaves and spikes of Brachypodium Bd21 and Bd3-1 lines. Our pathway enrichment analysis emphasised the accumulation of specialized metabolites representing the flavonoid biosynthetic pathway in parallel with processes related to nucleotide, sugar and amino acid metabolism. Similarities in metabolite profiles between both lines were relatively high in roots and leaves while spikes showed higher metabolic variance within both accessions. In roots, differences between Bd21 and Bd3-1 lines were manifested primarily in diterpenoid metabolism, while differences within spikes and leaves concerned nucleotide metabolism and nitrogen management. Additionally, sulphate-containing metabolites differentiated Bd21 and Bd3-1 lines in spikes. Structural analysis based on MS fragmentation spectra enabled identification of 93 specialized metabolites. Among them phenylpropanoids and flavonoids derivatives were mainly determined. As compared with closely related barley and wheat species, metabolic profile of Brachypodium is characterized with presence of threonate derivatives of hydroxycinnamic acids.


Introduction
Purple false brome (Brachypodium distachyon (L.) P. Beauv.; hereafter Brachypodium) is closely related to wheat and barley, making it potentially useful for functional genomics of these crops. Its main advantage as a model plant is the smallest genome found in the Poaceae family comprising five chromosomes spanning over 272 Mbp, in which about 25,000 protein-coding sequences are predicted [1]. In addition, Brachypodium is self-fertile and has a rapid life cycle of 8-10 weeks, depending on the environmental growth conditions [2]. A breakthrough point in Brachypodium research was the genome sequencing of accession Bd21 [3], which contributed to several genetic and genomic resources including Phytozome [4] and Gramene [5] and gave rise to initiatives like BrachyPan (Brachypodium pan-genome) [6]. Consequently, Brachypodium became an object of intense research in many fields serving in understanding interaction of grasses with viruses [7], bacteria [8], fungi [9] and invertebrates [10] as well as their responses to abiotic stresses [11].
insight into the global metabolite profile in analysed organs of both Brachypodium lines we performed principal component analysis (PCA) with all m/z signals detected during analyses performed with high resolution MS system in positive and negative ionization mode. The obtained PC3 plot revealed clear metabolic discrimination among tested organs and relative similarity between both lines (Figure 1). The highest consistency of metabolic profiles was observed within roots of both lines whereas the biggest interline differences were visible for spikes.

Comparison of Metabolomics Profiles in Analyzed Brachypodium Organs and Lines
Metabolic diversity within the studied Brachypodium lines and organs was represented in our LC/MS data sets by 22,307 individual signals detected in 48 analysed samples (three organs, two lines, two experiments and four biological replicates). To have a better insight into the global metabolite profile in analysed organs of both Brachypodium lines we performed principal component analysis (PCA) with all m/z signals detected during analyses performed with high resolution MS system in positive and negative ionization mode. The obtained PC3 plot revealed clear metabolic discrimination among tested organs and relative similarity between both lines ( Figure 1). The highest consistency of metabolic profiles was observed within roots of both lines whereas the biggest interline differences were visible for spikes. To corroborate our observations from the PCA plot, we used univariate two-way ANOVA analysis for each signal to classify signals into three following groups: (i) signals differentiating organs (O: comparison of the mean values of signal intensities from roots, spikes and leaves), (ii) signals differentiating lines (L: comparison of the mean values of To corroborate our observations from the PCA plot, we used univariate two-way ANOVA analysis for each signal to classify signals into three following groups: (i) signals differentiating organs (O: comparison of the mean values of signal intensities from roots, spikes and leaves), (ii) signals differentiating lines (L: comparison of the mean values of signal intensities from Bd21 and Bd3-1) and (iii) signals revealing significant interaction between organ and line factors (L×O: comparison of the mean values of signal intensities from Bd21 roots, spikes and leaves, and Bd3-1 roots, spikes and leaves) ( Figure 2A). As already indicated by the PCA plot (Figure 1), there was a relatively low number of signals discriminating lines whereas the majority of signals were organ specific. Nevertheless, the PCA plot was created on the basis of all detected signals while only fraction of them was filtered for O effect after ANOVA. This indicated a high impact of the differentiating signals on the entire metabolomic profiles in Brachypodium plants. Overall, these results are convergent with previous unbiased metabolom analyses conducted on the leaves and seeds of the Bd21 and Bd3-1 lines, which also revealed a stronger impact on organs than the genotype on metabolome [12].
Molecules 2022, 27, x FOR PEER REVIEW 4 of 35 signal intensities from Bd21 and Bd3-1) and (iii) signals revealing significant interaction between organ and line factors (L×O: comparison of the mean values of signal intensities from Bd21 roots, spikes and leaves, and Bd3-1 roots, spikes and leaves) ( Figure 2A). As already indicated by the PCA plot (Figure 1), there was a relatively low number of signals discriminating lines whereas the majority of signals were organ specific. Nevertheless, the PCA plot was created on the basis of all detected signals while only fraction of them was filtered for O effect after ANOVA. This indicated a high impact of the differentiating signals on the entire metabolomic profiles in Brachypodium plants. Overall, these results are convergent with previous unbiased metabolom analyses conducted on the leaves and seeds of the Bd21 and Bd3-1 lines, which also revealed a stronger impact on organs than the genotype on metabolome [12]. Differences in the metabolite set might contribute to the phenotypic differences between both studied lines, which have proven variation in many phenotypic traits [12,23,25], resistance to particular pathogens [26][27][28][29] or drought adaptation [17,32]. To obtain a better insight into the metabolic differences between Bd21 and Bd3-1 lines we selected signals corresponding to differentially accumulating metabolites (DAMs). We defined DAMs as signals significantly distinguishing both accessions (p-value < 0.01 for factors L or L × O) and differing at least two times (fold change; FC > 2) in their abundance in any of the tested organs of Bd3-1 and Bd21 accessions ( Figure 2B). Out of all signals, 2295 met these conditions for each organ suggesting a prevalent role for widely occurring elements in line differentiation. As suggested by the PCA plot, the proportion of DAMs indicated the lowest differences between Brachypodium lines in the roots and highest in the spikes. We selected 30 DAMs with the highest diversification among the studied groups to annotate respective m/z values and compare in detail differences in their abundances between particular lines and organs ( Figure 3). Despite good genetic characterization of Brachypodium, the metabolic pathways of this species are fragmentary in all dedicated metabolic platforms. Therefore, annotation of m/z values was per- (B) number of shared and organ-specific differentially accumulating metabolites (DAMs) defined as signals meeting the conditions: p-value < 0.05 for factor L or O × L; |log 2 (fold change)| > 1.5, where fold change was Bd3-1/Bd21 signal intensities.
Differences in the metabolite set might contribute to the phenotypic differences between both studied lines, which have proven variation in many phenotypic traits [12,23,25], resistance to particular pathogens [26][27][28][29] or drought adaptation [17,32]. To obtain a better insight into the metabolic differences between Bd21 and Bd3-1 lines we selected signals corresponding to differentially accumulating metabolites (DAMs). We defined DAMs as signals significantly distinguishing both accessions (p-value < 0.01 for factors L or L × O) and differing at least two times (fold change; FC > 2) in their abundance in any of the tested organs of Bd3-1 and Bd21 accessions ( Figure 2B). Out of all signals, 2295 met these conditions for each organ suggesting a prevalent role for widely occurring elements in line differentiation. As suggested by the PCA plot, the proportion of DAMs indicated the lowest differences between Brachypodium lines in the roots and highest in the spikes. We selected 30 DAMs with the highest diversification among the studied groups to annotate respective m/z values and compare in detail differences in their abundances between particular lines and organs ( Figure 3). Despite good genetic characterization of Brachypodium, the metabolic pathways of this species are fragmentary in all dedicated metabolic platforms. Therefore, annotation of m/z values was performed with a database created with Oryza sativa subsp. japonica (japonica rice), described at the metabolome level model plant from Poaceae family [35].
Among the annotated compounds putative derivatives of hydroxycinnamic acids (feruloylhydroxycitric acid, caffeoylpyruvylhexose, isomers of caffeoylthreonic acid and cinnamic acid ethyl ester) were highly represented ( Figure 3). These included conjugates of hydroxycitrate with hydroxycinnamic acids known from Zea mays as compounds with high variation in accumulation profile in different inbred lines [36]. The signal corresponding to feruloylhydroxycitric acid had the highest abundance in Bd3-1 roots compared with the Bd21 roots. The relatively high level of caffeoylthreonic acid isomers in all organs of the Bd21 line, as compared with Bd3-1, is noteworthy. The same trend of high abundance in Bd21 line was observed for hydroxybenzoic acid derivatives (N-salicyloylaspartic acid and N-pyruvoyl-methoxy-hydroxyanthranilic acid).
Putative derivatives of the flavone apigenin (isovitexin pentose-deoxyhexoside, pentahydroxy-dimethoxyflavone hexoside and apigenin hydroxy-methylglutaryl-hexoside) together with proanthocyanidin B and cyanidin acylated glycoside were representatives of differentiating flavonoids. Interestingly, cyanidin 3-O-glucoside (chrysantemin) has been already reported as differentially accumulating metabolite in Brachypodium spikes [13]. These correlative findings suggest that biosynthesis of cyanidin glycosides clearly discriminate both Brachypodium lines. Despite the common biosynthetic origin, differences in abundance of particular flavonoids and hydroxycinnamic acids were not correlated. However, it should be noted that most of the distinctive compounds were complex structures that were relatively distant from the common precursors in the metabolic pathway. This in turn indicated that the activities of the enzymes responsible for particular modifications of the core structures, including hydroxylation, acylation, methylation and glycosylation, were responsible for the observed differences in phenylpropanoid metabolism among compared organs and lines.

Most Represented Metabolic Pathways
For further global settling of Brachypodium metabolome into the biological context, functional analysis of obtained result on the basis of MetaboAnalyst was implemented. Firstly, pathway-level enrichment has been performed with all m/z signals detected in positive and negative ionization modes for overall picture of metabolites in Brachypodium plants ( Table 1). The same as for metabolite annotation, our analysis was performed on O. sativa database [35]. Direct tentative annotation of m/z values obtained during our analysis to rice metabolites enabled further calculation of pathway-level enrichment.
Our analysis performed with all signals indicated significant enrichment of flavonoidrelated pathways (Table 1). Forty signals have been matched to metabolites from flavonoids biosynthesis (47 metabolites in total). This was accompanied by signals matching 12 metabolites (all from the same pathway) from flavone and flavonol biosynthesis. This indicated a significant contribution of a specialized metabolism to the overall profile of Brachypodium metabolites.  [35,37]. l-leaves; r-roots; s-spikes.

Most Represented Metabolic Pathways
For further global settling of Brachypodium metabolome into the biological context, functional analysis of obtained result on the basis of MetaboAnalyst was implemented. Firstly, pathway-level enrichment has been performed with all m/z signals detected in positive and negative ionization modes for overall picture of metabolites in Brachypodium plants ( Table 1). The same as for metabolite annotation, our analysis was performed on O. sativa database [35]. Direct tentative annotation of m/z values obtained during our analysis to rice metabolites enabled further calculation of pathway-level en-  [35,37]. l-leaves; r-roots; s-spikes.  [37]; total-number of compounds included in biological pathway in the database; hits-number of compounds matched in our analysis; FDR-false discovery rate; impactpathway impact value related to the number of links occurred upon a node in pathway topology graph. Full set of data including other pathways as well as lists of annotated metabolites is available as Supplementary The remaining identified pathways represented primary metabolic processes, mainly nucleotide, sugar and amino acid metabolism. Twenty-five metabolites from the galactose metabolism matched with signals from our analysis, including galactinol and raffinose, which were previously described as involved cold and drought stress response in Brachypodium [38]. Significant annotation of pentose and glucuronate interconversions was mainly related to metabolites from modules of pectin degradation and glucuronate pathway. Key components of plant metabolism from the pentose phosphate pathway as source of substrates for synthesis of purine nucleotides, followed by purine metabolism, were highly matched. In the second mentioned pathway, annotation focused on modules of inosine monophosphate biosynthesis, adenine ribonucleotide biosynthesis and purine degradation. Pyridoxal and pyridoxine, as well as their phosphorylated derivatives from vitamin B6 metabolism, were also significantly enriched. Antioxidant and defence related metabolites from amino sugar and nucleotide sugar metabolism and ascorbate and aldarate metabolism, including phosphorus containing structures (D-Glucosamine phosphate and their derivatives and UDP-glycosides) were significantly annotated.
Amino acid metabolism branched-chain amino acids, valine, leucine and isoleucine biosynthesis were matched. Related to them 2-Oxocarboxylic acid metabolism was also highly scored in both statistical significance and pathway impact. Compounds annotated to this pathway focused on a branch of the 2-Oxocarboxylic acid chain extension by tricarboxylic acid module, which in Brassicaceae species leads to glucosinolate biosynthesis [39]. However, in Poaceae, glucosinolates are absent, therefore, signals annotated to sulphide compounds such as 2-(5'-methylthio)pentylmalic acid and isomer 3-(5'methylthio)pentylmalic acid), 2-(6'-methylthio)hexylmalic acid and isomer 3-(6'-methylthio) hexylmalic acid from this branch can be components of different pathways. Interestingly, elements of "2-Oxocarboxylic acid" were previously reported in grasses as factor involved in the response to salinity and drought stress [40,41]. Further inspection of LC-MS and MS/MS spectra showed the presence of such S-containing compounds in Brachypodium plants not previously identified, which confirms the validity of enrichment analysis in plant metabolite profiling (Supplementary Figure S1). Another amino acid related pathway, tyrosine metabolism, was selected based on the annotation of tyrosine and 3,4-dihydroxyphenylalanine, which, in grasses, plays a key role in lignin biosynthesis (Maeda, 2016).

Metabolic Pathways Distinguishing Bd21 and Bd3-1 Lines
In order to identify metabolic pathways discriminating on both analysed lines in particular organs, pathway enrichment analysis was only performed for signals representing DAMs selected based on the above-described ANOVA analysis (p-value ≤ 0.01 for factor L or O × L; FC > 2) ( Figure 2B, Table 2). Housekeeping and general metabolism-related biological pathways (galactose metabolism, pentose phosphate pathway, valine, leucine and isoleucine biosynthesis) highly varied among the six compared groups. Moreover, specialized metabolism of flavonoid biosynthesis, was related to flavone and flavonol biosynthesis and 2-Oxocarboxylic acid metabolism were commonly differentiated these groups.
Metabolic differences between the roots of Bd21 and Bd3-1 are manifested primarily by "Diterpenoid biosynthesis". In monocots, diterpenoids are known from large structural diversity and species-specificity. In Brachypodium, no specialized diterpenoids have been identified, however, a few Brachypodium genes are homologous to rice genes related to momilactone phytoalexin production [42]. In our analysis, the main differentiating module of diterpenoid biosynthesis was gibberellin production, which is the key factor in root elongation in monocots [43]. The differences in gibberellin levels in the roots of both Brachypodium lines could be related to observed differences in root morphology of both tested lines [25]. In this context it was also of interest that tetrahydrofolate from differentiating one carbon pool by folate pathway has been reported as key regulators of root development [44].
The caffeine metabolism pathway including purine alkaloids was also significantly different among the roots of both Brachypodium lines. Matched intermediates of this pathway included xantosine, 7-methyluric acid and their derivatives. However, caffeine itself has been not reported in grasses while at least some of the matched metabolites can be linked with purine salvage or degradation [45].
A pathway that specifically differed between spikes of both Brachypodium lines was the cysteine and methionine metabolism indicating possible differences in sulphate assimilation. S-Adenosyl-L-methionine, a key metabolite from this pathway, is a donor of methyl group in numerous transmethylation reaction influencing physical and chemical properties of lignin polymers, as well as hydroxycinnamic acids synthesis in Brachypodium plants [46].

Metabolite Identification with LC-MS Systems
In the next step of our study, we tried to identify a subset of detected metabolites based on their spectra obtained during the HPLC-ESI-MS n and UPLC-HR-MS/MS analyses. MS n spectra are helpful in the identification of complex metabolites, for example flavonoids glycoconjugates, where they can enable the determination of the place and character of the glycosidic bond. In addition, the order of detachment of individual fragments from complex structures with a simultaneous observation of the intensities of particular product ions enables the differentiation of isomeric and isobaric structures, unlike MS/MS, which, in many cases hampers, isomers differentiation. However, accurate measurement of m/z values obtained during HR MS/MS analysis allowed confirmation of tentative structures predicted by the MS n analysis. Overall, this analysis enabled us to identify 93 metabolites at levels 1-3 according to the Metabolomic Standards Initiative [47] (Table 3). This manual metabolite identification enabled us to describe the structures specific to Brachypodium, which are absent in metabolomic databases and, therefore, cannot be annotated with automated bioinformatics approaches. Table 2. Metabolic pathways discriminating Bd21 and Bd3-1 lines in particular organs selected based on pathway enrichment analysis performed with MS signals representing differentially accumulating metabolites (DAMs). KEGG-Kyoto Encyclopedia of Genes and Genomes [37]; total-number of compounds included in biological pathway in the database; hits-number of compounds matched in our analysis; FDR-false discovery rate; impact-pathway impact value related to the number of links occurred upon a node in pathway topology graph. Full sets of data including other pathways as well as lists of annotated metabolites are available as Supplementary

Hydroxycinnamoyl-Quinic Acids
MS signals corresponding to phenolic and hydroxycinnamic acids and their derivatives were detected in all studied organs at high abundances. p-coumaric, caffeic and ferulic acids have been identified in Brachypodium in conjugation with (methyl)quinic acids, sugar acids, or polyamines; while caffeic and ferulic acids (Table 3; compounds 21, 38) have been additionally observed as free molecules that have been identified by comparison to their standards.
Different isomers of hydroxycinnamoyl-quinic acids (HQA) were reported and can be found in metabolite databases, but proper annotation of these isomeric structures should be supported by fragmentation schemes in MS/MS or MS n spectra. In our analysis, HQAs (8, 26, 29, 43 and 48) were identified according to Clifford et al. [43] and Piasecka et al. [54]. The main product ion from deprotonated molecules of compounds 29, 43 and 48 at m/z = 173 Da corresponded to quinic acid, which is distinctive for 4HQA isomers ( Figure 4A). On the other hand, 5HQAs are characterized by the main product ion representing the respective hydroxycinnamic acid molecule [56] as we found for 26 where product ion at m/z = 193 Da corresponded to ferulic acid ( Figure 4B). Deprotonated ions of compounds 8 had the same exact m/z value as compounds 26 and 43 (367.10425 Da) with the same chemical formula C 17 H 20 O 9 calculated from exact mass. The UV spectra of 8 with a maximum of absorption at 305 nm (Table S2)

Esters of Hydroxycinnamic and Threonic Acids
In addition to quinic acid conjugates, the most abundant signals detected in all studied organs corresponded to sugar acid, mainly threonic acid, derivatives.  Overall, we observed three different hydroxycinnamic acids conjugated with quinic acid as 4HQA isomers (29, 43 and 48), while only ferulic acid conjugate has been detected as 5HQA isomer (26). However, the dominant presence of 4HQAs might result from the isomerization of 5HQAs that could occur either in vivo or in the extract solution [63]. Widely occurring grass 3HQA isomers [56] have been not detected in our study.

Esters of Hydroxycinnamic and Threonic Acids
In addition to quinic acid conjugates, the most abundant signals detected in all studied organs corresponded to sugar acid, mainly threonic acid, derivatives.  Figure 5B). Interestingly, according with our pathway enrichment analysis present in these highly abundant conjugates (18,42,63), threonic acid could be generated as a degradation product of ascorbate in plants [64].

Hydroxycinnamic Acid Amides
Our analysis of MS n and MS/MS spectra also enabled identification of conjugates of p-coumaric, caffeic and ferulic acid with putrescine and agmatine. Among these, compounds 10 and 12 have similar masses of protonated [M+H] + ion at m/z = 235 Da and calculated chemical formula from exact mass as C13H17O2N2 ( Figure 6A). Compounds 10 and 12 Also share product ion at m/z = 147 Da that correspond to dehydrated p-coumaric acid. In both compounds losses of -NH4 group and of entire putrescine moiety can be

Hydroxycinnamic Acid Amides
Our analysis of MS n and MS/MS spectra also enabled identification of conjugates of pcoumaric, caffeic and ferulic acid with putrescine and agmatine. Among these, compounds 10 and 12 have similar masses of protonated [M+H] + ion at m/z = 235 Da and calculated chemical formula from exact mass as C 13 H 17 O 2 N 2 ( Figure 6A). Compounds 10 and 12 Also share product ion at m/z = 147 Da that correspond to dehydrated p-coumaric acid. In both compounds losses of -NH 4 group and of entire putrescine moiety can be also observed ( Figure 6BC). Fragmentation of compound 12 was characterized by additional product ion at m/z = 119 Da, which is typical for p-coumaric acid with preserved peptide bond between acidic and polyamine substituents. also observed ( Figure 6BC). Fragmentation of compound 12 was characterized by additional product ion at m/z = 119 Da, which is typical for p-coumaric acid with preserved peptide bond between acidic and polyamine substituents.  Regarding this, we tentatively assigned compounds 10 and 12 as p-coumaroyl-Nputrescine isomeric structures. The same scheme of fragmentation and differences in product ion intensities were observed for compounds 11 and 23. The protonated [M+H] + ions of both metabolites yielded losses of fragment 88 Da, corresponding to putrescine, whereas the main product ions at m/z = 248 and 177 Da in an MS2 scan indicated a ferulic acid molecule. In an MS3 scan of compound 23, an additional product ion at m/z = 145 Da indicated peptide bond preservation during fragmentation steps. The calculated chemical formula of compounds 11 and 23 and fragmentation spectra corresponded to C 14 H 21 O 3 N 2 , which complied with feruloyl-N-putrescine isomers, similar to 10 and 12. In an MS3 scan of compound 23, an additional product ion at m/z = 145 Da indicated peptide bond preservation during fragmentation steps. In grasses, two geometric isomers of hydroxycinnamic acids cis and trans were previously reported [71]. Trans isomers constitute the predominantly occurring form in plants, but they can be transformed to corresponding cis isomers by UV-radiation. Concerning this, we assumed that the isomeric hydroxycinnamic acid amides revealing differences in their fragmentation pattern represented cis-trans conformers. were determined as p-coumaroylagmatine and feruloylagmatine. The latter compound has already been reported in Brachypodium leaves [72]. Agmatine derivatives were also observed in barley, including complex hordatine A, B and C structures that possess antifungal activities [58]. Agmatine conjugates with hydroxycinnamic acids were also reported in wheat [73] and species representing other plant families like African shrub Maerua edulis in which cis-trans conformers of p-coumaroylagmatine were distinct [74].

Flavonoid Glycosides
The highest number of compounds identified in this study belonged to flavonoids. MS n fragmentation schemes of the same or similar molecules were previously described in Poaceae plants including barley and wheat [58,62,75]. Therefore, structural similarities and conservation of structures can be deduced for those closely related species. Our analysis resulted in identification of isomeric structures of glycosylated flavonols (quercetin, isorhamnetin), flavones (apigenin, luteolin and chrysoeriol) and proantocyanidins. The aglycone type of molecular mass 270 and 286 Da were further confirmed as apigenin and luteolin, respectively, referring to pseudo-MS3 spectrum described previously [76]. The characteristic product ions at m/z 199, 175, 151 and 133 determined the aglycone as luteolin whereas the characteristic product ions at m/z 153, 145, 121 and 119 determined the aglycone as apigenin. In roots, derivatives of O-methylated flavonoids like isorhamnetin, tricin and chrysoeriol were mainly identified. Glycosides of all these aglycones were present in Poaceae species in diverse structural configuration [58,62,75]. We mainly found hexose(s), deoxyhexose(s) and pentose(s) sugar substituents of analyzed flavonoid aglycones. However, a precise sugar structure could not be determined by MS analysis. The presence of glucose and galactose could be suggested in Brachypodium flavonoid glycoconjugates as both hexoses were found in flavonoids of closely related species barley and wheat [77,78]. Among other sugar only rhamnose as deoxyhexose and arabinose as pentose were described in flavonoids of Poaceae family.
Among the different aglycone-sugar conjugates, we identified O-glycosides and C-glycosides as well as O-,C-glycosides. In addition, LC-MS n analysis of flavonoid diglycosides enabled distinguishing interglycosidic bonds in glycosyl(1→6)glycosides and glycosyl(1→2)glycosides [58,62]. Unfortunately, the exact positions of aglycone substitutions in flavonoid glycosides were difficult to establish without NMR analysis. However, on the basis of similarities in fragmentation scheme with mass spectra of flavone derivatives reported in species closely related to Brachypodium, including barley and wheat, 4-OHand 7-OH-groups could be suggested as positions of glycosylation of flavonoid derivatives from Brachypodium [58,62,75]. ions were detected at a lower abundance, suggesting a deoxyhexosyl(1→6)hexosidic bond in these structures according to [79]. The major protonated product ions at m/z = 331.08045 and m/z = 301.07037 Da corresponded to tricin and chrysoeriol aglycones, respectively. A relatively high, intense [M+H-164] + ion in MS2 spectrum of metabolite 90 proved that deoxyhexose is external sugar in this molecule. Furthermore, the detachment of an entire moiety of deoxyhexose was typical for deoxyhexosyl(1→2)glycosyl interglycosidic bond.  Figure 7A). The identification of compounds 57, 68, 69 and 76 highlighted the problem with the automatic annotation of metabolomic signals resulting from the diversity of isomeric or isobaric structures available in mass spectra databases. In the case of these four compounds (m/z = 593.15155), we found 110 entries (with error ppm = 5) in the Metlin database [81]. For this reason, identification of such isomeric structures could be only done by MS/MS or MS n fragmentation scheme analysis.

O,C-Glycosides
Structural similarity in C-glycosides of apigenin and luteolin between Brachypodium, barley wheat and maize were significant [58,75,82]. However, C-glycosylation could be catalysed by different enzymes in Poaceae plants, as for example, by UDP-glucose-dependent C-glucosyltransferase in rice and wheat [83] or by bifunctional C-/O-glycosyltransferase in maize [82]. In these cereals as well as in Brachypodium, C-glycosides of the flavones apigenin and luteolin were dominant metabolites, with glycosylation occurring singly or doubly at the 8-C and 6-C positions. Structural isomerism related to glycoconjugation hindered proper C-glycosides identification in Brachypodium plant by simple annotation to dedicated databases. Furthermore, significant changes in content of saponarin isomers can be observed during plant development as was detected for barley [84]; therefore, detailed flavonoids studies in Brachypodium should be further extended in developmental context.

Acylated Flavonoids
Acylated flavones, in which glycosides are substituted with hydroxycinnamic acids, i.e., p-coumaric, caffeic, ferulic and sinapic acids, present in barley and wheat have been studied by NMR and mass spectrometry and were identified as 7-O-[6 -acyl]-glucosides and 7-O-[6 -acyl]-glucosyl-4 -O-glycosides [85][86][87]. Surprisingly, in Brachypodium we identified only one acylated flavonoid (60) that was detected in leaves while in Poaceae family different isomeric and isobaric structures are present. According to literature data from other related species the structure was established as isoorientin 2 -O-hexoside 7-O-[6 -sinapoyl]-hexoside [58]. The absence of acylated flavonoids in our analysis may eventually arise from our plant growth conditions and/or age of plants used in our analysis, therefore presence of acylated flavonoids in Brachypodium cannot be excluded.

Flavan-3-ols
Our LC/MS analysis revealed that spikes of Brachypodium were rich in flavan-3ol derivatives, mainly from the proanthocyanidin group. Stereoisomers catechin and epicatechin as well as gallocatechin and epigallocatechin showed the same product ions and very similar ratios on the MS/MS or MS n spectra, thus even preliminary identification of these isomers was not possible based on the obtained MS spectra. The position and the stereochemistry of the interflavan linkage could not be elucidated by MS. Procyanidins can be classified into A-type and B-type depending on the stereo configuration and linkage between monomers. B-type procyanidins possess a single C-C interflavan bond while A-type procyanidins have additional ether bond. The (epi)catechin, (epi)gallocatechin and galloyl subunits were identified in Brachypodium as procyanidin A-and B-type on the basis of typical fragment detachments in MS/MS of deprotonated molecules according to [88,89]. High resolution mass spectrometry as well as MS n enabled to trace the way of fragmentation of dimeric and trimeric structures of (epi)catechin and (epi)gallocatechin proanthocyanidin. An example of the fragmentation of compound 15 is the hRetro Diels-Adler reaction (RDA), which resulted in the detachment of fragment 168 Da, typical for (epi)gallocatechin oligomers (Figure 8). The most abundant product ion at m/z = 407 Da resulted from water elimination from RDA product and [M−H-126] − ions from heterocyclic ring fission (HRF). In the parallel fragmentation, a Quinone-Methide cleavage of interflavan bond occurred and the remaining deprotonated ions corresponded to (epi)catechin monomer. Therefore, 15 was identified as a procyanidin B-type dimer. Compounds 7 and 27 possessed a similar fragmentation scheme as compound 15; however, their m/z ratios indicated A-type interflavan bonds in those proanthocyanidis. Trimeric proanthocyanidins (6, 7 and 20) were also observed. The typical RDA and QM fragments corresponded to (epi)catechin and (epi)gallocatechin subunits in those structures. fragments corresponded to (epi)catechin and (epi)gallocatechin subunits in those structures.
Proanthocyanidins are well described in several cereal and grass species, including sorghum and barley, due to their nutritional and technological importance [90]. For instance, barley grains have been shown to be especially rich in proanthocyanidins, which contribute to haze formation in barley beer [91]. However, despite well-characterized Brachypodium genes responsible for catechin and epicatechin biosynthesis [92] the number of studies on the accumulation of proanthocyanidins in Brachypodium plants is very limited. Only (-)-epicatechin was detected by widely targeted metabolome analysis in seeds and leaves of Bd21 and Bd3-1 [12]. Overall, our study revealed differential accumulation of proanthocyanidins among Brachypodium organs. We observed high accumulation levels of these compounds in seeds and intermediate accumulation level in leaves while only compounds 6 ((epi)gallocatechin trimer) and 24 ((epi)gallocatechin O-hydroxybenzoate) were detected in roots. Di-to penta-mers and galloylated proanthocyanidins were found in barley. Epiafzelechin derivatives, which are abundant in buckwheat (Fagopyrum esculentum) [93], were not detected during our analysis.  Proanthocyanidins are well described in several cereal and grass species, including sorghum and barley, due to their nutritional and technological importance [90]. For instance, barley grains have been shown to be especially rich in proanthocyanidins, which contribute to haze formation in barley beer [91]. However, despite well-characterized Brachypodium genes responsible for catechin and epicatechin biosynthesis [92] the number of studies on the accumulation of proanthocyanidins in Brachypodium plants is very limited. Only (-)-epicatechin was detected by widely targeted metabolome analysis in seeds and leaves of Bd21 and Bd3-1 [12]. Overall, our study revealed differential accumulation of proanthocyanidins among Brachypodium organs. We observed high accumulation levels of these compounds in seeds and intermediate accumulation level in leaves while only compounds 6 ((epi)gallocatechin trimer) and 24 ((epi)gallocatechin O-hydroxybenzoate) were detected in roots. Di-to penta-mers and galloylated proanthocyanidins were found in barley. Epiafzelechin derivatives, which are abundant in buckwheat (Fagopyrum esculentum) [93], were not detected during our analysis.

Plant Material
Seeds of Brachypodium lines Bd21 and Bd3-1 were obtained from Robert Hasterok (University of Silesia, Katowice, Poland). Plants were grown in soil in a controlled growth chamber at 23 or 20 • C (day or night) under short day conditions (8 h light and 16 h darkness) for 4 weeks at 50-60% relative humidity and irradiance of 100 µmol m −2 s −1 . Seeds were sown in 9 × 9 cm square pots filled with soil mixed with peat. In the third week of cultivation a supplemental liquid fertilizer (N total 12% w/w, P 2 O 5 4% w/w, K 2 O 6% w/w, B 0.01% w/w, Cu 0.0007% w/w, Fe 0.015% w/w, Mn 0.012% w/w, Mo 0.001% w/w, Zn 0.005% w/w) was applied in a concentration of 1 mL/1 L of water. Subsequently, light conditions were changed to a long day (16 h light and 8 h darkness) and plants were grown until the developmental stage of spikes, i.e, watery ripe: first grains have reached half their final size (71-73 in BBCH scale according to [94]). Samples of the spikes, leaves and roots were collected, weighted and placed in 2ml tubes containing zirconia beads, then immediately frozen in liquid nitrogen and stored at −80 • C until further processing. An extraction buffer containing 0.5 mM lidocaine and 0.5 mM camphorsulfonic acid in DMSO was added (2.5 µL /1 mg of fresh plant weight) to each tube. Samples were homogenized with a Precellys Evolution (Bertin, France) tissue grinder and centrifuged in 4 • C at 15,000 g. Supernatants were collected and directly subjected to LC-MS analysis.

Metabolite Profiling
Analysis and identification of metabolites was performed using two complementary LC-MS systems using the previously published approach [58]. First of them (low resolution HPLC-DAD-MS n ) consisted of 1100 HPLC system with a photodiode-array detector (Agilent, Santa Clara, CA, USA) equipped with an XBridge Shield C18 column (150 × 2.1 mm, 3.5 µm particle size; (Waters, Milford, CT, USA) coupled to an Esquire 3000 ion trap mass spectrometer (Bruker Daltonics, Billerica, MA, USA). Chromatographic separations were conducted with injection volume 10 µL using water with 0.1% formic acid (solvent A) and acetonitrile (solvent B) at the flow rate 0.2 mL/min and the following gradient: 0-6 min from 8% to 10% B, 6-40 min to 20% B, 40-46 min to 98% B maintained for 5 min. The MS n spectra were separately recorded in the negative and positive ion modes. The second system (high resolution UPLC-MS/MS) consisted of UPLC equipped with a photodiode-array detector (Acquity System; Waters) hyphenated to a high-resolution Q-Exactive hybrid MS/MS quadrupole Orbitrap mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). Chromatographic separation was performed on an Acquity UPLC HSS T3 C18 chromatographic column (2.1 × 50 mm, 1.8 µm particle size; Waters) at 22 • C using water containing 0.1% formic acid (solvent A) and acetonitrile (solvent B). The injection volume was 5 µL, gradient elution started at 100% of A and linearly changed to 20% of B over 2 min, then to 30% of B over 8 min and to 95% of B over 1 min maintained for 2 min. UV absorbance was recorded in the 230-450 nm wavelength range with a resolution of 2 nm. Q-Exactive MS operated in Xcalibur version 3.0.63 with the following settings: heated electrospray ionization ion source voltage −3 kV or 3 kV; sheath gas flow 30 L/min; auxiliary gas flow 13 L/min; ion source capillary temperature 250 • C; auxiliary gas heater temperature 380 • C. MS/MS mode (data-dependent acquisition) was recorded in negative and positive ionization, at resolution 70,000 and AGC (ion population) target 3e6, scan range 80 to 1000 m/z.

Metabolite Identification
Individual compounds were tentatively identified on the basis of low resolution LC-MS n mass spectra if corresponding fragmentation spectra and m/z signals were confirmed in high resolution UPLC-MS/MS. Particular structures were suggested via comparison of the exact molecular masses with ∆ less than 5 ppm, mass spectra and retention times to those of standard compounds, spectra in available databases (PubChem, ChEBI, Metlin, Reaxys, DynLib and KNApSAck) [48,66,81,[95][96][97] and literature data. Confirmation of isomeric aglycone type was based on available standard compounds and methods described previously [76]. Pseudo-MS3 spectra of flavonoids O-glycosides with in-source CID 80 eV enabled to confirm fragmentation typical for luteolin based on the presence of product ions at m/z 117 and 135 for and excluded presence of other isomers e.g., kaempferol.

Bioinformatic Processing
High-resolution raw UPLC-MS/MS data were separately processed by MZmine 2.53 [98] for negative and positive ionization. In first step, lists of masses were generated in each scan of the raw data files (Supplementary File S1). Chromatograms for each exact mass detected over the scans were built by a Chromatogram Builder algorithm. These chromatograms were deconvoluted using an ADAP Wavelets algorithm and subsequently subjected to isotope elimination, adduct and complex searching, followed by retention times normalization among peak lists. Such transformed peaks were aligned across all samples by a Join aligner module. The resulting peak table was completed by supplemental peak detection with a peak finder algorithm prior to missing value imputation (gap-filling). The obtained result table was subjected to further statistical analysis and visualizations.

Statistical Analysis
Statistical analysis was performed with a Genstat 21 (VSN International, Hempstead United Kingdom). Observations below the detection limit were substituted with half of the minimum non-zero observation for each metabolite and then observations were transformed by log 2 (x). Two-way analysis of variance (ANOVA) was performed with the experiments as a block (random effects) and organ line as two fixed factors. Analysis was performed together with positive and negative ionization. Significant changes in the accumulation of metabolites was indicated by the effect on an organ, line or by the interaction an organ x line with p-value < 0.05. To select metabolites with significant differences between lines in each organ, the definition of DAM was introduced. Only signals with a significant interaction of line and organ (LxO) or with significant effects on a line (L) were classified as DAMs. Additionally, for each signal, we calculated fold change (FC) in each organ as a ratio of signal intensities in the Bd3-1 and Bd21 lines (Bd3-1/Bd21) to restrict analysed DAMs to those with FC > 2 (|log 2 (Bd3-1/Bd21)| > 1). Visualizations including PCA 3D plot (generated with data after log 2 transformation), heatmaps and Venn diagrams were created in R (R Foundation for Statistical Computing).

Pathway Enrichment Analyses
Pathway enrichment analyses were conducted with all m/z signals from the combined positive and negative ionization modes and only with m/z signals representing DAMs. Data was imported to a functional analysis module in MetaboAnalyst 5.0 [40]. This enabled direct m/z value annotation to metabolic data base for O. sativa subsp. Japanese from the Kyoto Encyclopedia of Genes and Genomes (KEGG) [35,37]. Signals with significant annotation were further subjected to pathway-level enrichment on the basis of mummichog algorithms in a pathway analysis module of MetaboAnalyst 5.0. This module filtered metabolites overrepresented on the pathway level in addition to pathway topology analysis. Enrichment was ranked on the basis of mummichog algorithms followed by Benjamini-Hochberg false discovery rate (FDR) correction. Pathway topology was scored on the basis of relativebetweenness centrality, to estimate the relative importance of individual nodes to the overall pathway network. The node impact values were normalized by the sum of the importance of the pathway to estimate maximum impact of each pathway as 1. Significantly enriched metabolic pathways upon differentiating factors were selected if FDR < 0.03 and pathway impact > 0.3 were consistent across multiple comparisons.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/molecules27185956/s1, Figure S1: High resolution MS/MS spectra of sulphur-containing dicarboxylic acids; Figure S2: Representative UV chromatograms; Table S1: Scoring of pathway mapping; Table S2: Pathway mapping conducted on m/z signals significantly differentiating roots of Bd21 and Bd3-1 lines; Table S3: Pathway mapping conducted on m/z signals significantly differentiating leaves of Bd21 and Bd3-1 lines; Table S4: Pathway mapping conducted on m/z signals significantly differentiating spikes of Bd21 and Bd3-1 lines; File S1. Parameters of raw MS/MS data processing by MZmine software.

Data Availability Statement:
Most of the data used in the current study is contained within the article. Remaining data that support the findings of this study is available from the corresponding authors upon request.