The Metabolic Profile of Anchusa officinalis L. Differs According to Its Associated Arbuscular Mycorrhizal Fungi

Anchusa officinalis (L.) interacts with various microorganisms including arbuscular mycorrhizal fungi (AMF). Recently, the AMF Rhizophagus irregularis MUCL 41833 has been shown to modulate the metabolome of A. officinalis. However, little information is available on the impact that different AMF species may have on primary and secondary plant metabolites. In this study, four AMF species belonging to the genus Rhizophagus (R. irregularis MUCL 41833, R. intraradices MUCL 49410, R. clarus MUCL 46238, R. aggregatus MUCL 49408), were evaluated for their potential to modulate A. officinalis metabolome under controlled semi-hydroponic cultivation conditions. An untargeted metabolomic analysis was performed using UHPLC-HRMS followed by a multivariate data analysis. Forty-two compounds were reported to be highly modulated in relation to the different AMF associations. Among them, six new secondary metabolites were tentatively identified including two acetyl- and four malonyl- phenylpropanoid and saponin derivatives, all presenting a common substitution at position C-6 of the glycosidic moiety. In addition, an enhanced accumulation of primary and secondary metabolites was observed for R. irregularis and R. intraradices, showing a stronger effect on A. officinalis metabolome compared to R. clarus and R. aggregatus. Therefore, our data suggest that different AMF species may specifically modulate A. officinalis metabolite production.


Introduction
The interactions between plants and beneficial microorganisms improve not only growth and health, but also modify the metabolome considerably [1,2]. Among these microorganisms, arbuscular mycorrhizal fungi (AMF) are known for their intimate symbiotic relationship with circa 72% of land plants [3]. They facilitate phosphate and nitrogen uptake in exchange for carbohydrates [4,5], offer enhanced resistance to pests and diseases, drought and salinity, and are involved in the heavy metal detoxification process [6].
An increasing number of studies (reviewed extensively by Zeng et al. [7] and by Kaur and Suseela [2]) have reported the impact of AMF on the production of primary and secondary metabolites in different plant species, including crops. For instance, a higher accumulation of amino acids, fatty acids, isoflavonoids and phenylpropanoid derivatives was reported in roots of Medicago truncatula colonized by Rhizophagus irregularis [8]. Similarly, a significant upregulation of secondary metabolites belonging mainly to blumenol derivates and pyrrolizidine alkaloids, was detected in the roots of Senecio jacobaea associated with R. irregularis [9].
Not necessarily all plant-AMF associations result in increased metabolite production. The productivity of the host plant-i.e., the outcome of the plant-microbe interactionis highly dependent on the identity of the AMF symbiont [2]. Indeed, different AMF species can differently affect the production of specific phytochemicals on a specific plant species. For instance, Rivero et al. [10] reported that, although the metabolic pathways of Solanum lycopersicum altered by mycorrhizal symbiosis were common for both Funneliformis mosseae and R. irregularis, the compound's biosynthesis was altered depending on the particular AMF species involved. Jasmonic acid (JA), methyl-JA and JA-Isoleucine conjugates were accumulated in significantly higher amounts only in the plants associated with F. mosseae. This suggests that specific AMFs have a fine-tuned regulation role in biosynthetic pathways.
Besides the increasing interest in plant-AMF associations, those studies applying a metabolomic approach were limited to only two AMF species, R. irregularis (formerly Glomus intraradices) and F. mosseae (formerly Glomus mosseae) [8][9][10][11][12][13]. In a recent study, Cartabia et al. [14] described the effects of R. irregularis MUCL 41833 on shoots and roots of A. officinalis, growing under a semi-hydroponic cultivation system for a period of 9 days. The untargeted metabolomic approach showed an upregulation of 35 primary and secondary metabolites (e.g., organic acids, phenolic compounds, oleanane-types glycosides) in mycorrhized plants as compared to the non-mycorrhized ones. Interestingly, A. officinalis is always found associated with AMF in the wild [15]. Thus, considering non-mycorrhized plants as control does not truly reflect the natural conditions. Therefore, in the present study, a metabolomic analysis of roots and shoots of A. officinalis associated with four different AMF species belonging to the same genus (R. irregularis MUCL 41833, R. intraradices MUCL 49410, R. clarus MUCL 46238 and R. aggregatus MUCL 49408) was conducted to test the hypothesis that different AMF strains impact the plant metabolome differently. R. irregularis MUCL 41833 was chosen as the reference strain since it is considered as a model organism in AMF research and the most used AMF strain in commercial inoculants [10,16]. It is a generalist colonizer, present in almost all soils and climatic zones [16], readily grown in vitro on root organs [17] and its whole genome sequence has been published recently [18,19]. Moreover, its effects on the metabolome of A. officinalis were recently published [14].
To test our hypothesis, the plants were grown for 9 days in the semi-hydroponic cultivation system described by Cartabia et al. [14] and an ultrahigh-performance liquid chromatography high-resolution mass spectrometry (UHPLC-HRMS) analysis was performed. An untargeted metabolomics approach was further conducted to shed light on the overall effects of the different AMF species on A. officinalis root and shoot parts.

Root Colonization by AMF and Plant Total Fresh Weight
Intraradical AMF structures were assessed in A. officinalis plants treated with four different AMF species (R. irregularis, R. intraradices, R. aggregatus and R. clarus) both at the plants' transfer in the containers (T0) and at the end of the experiment (i.e., after 9 days of growth in the system-T1). The overall data for AMF colonization and the total fresh weight (TFW, e.g., the sum of roots and shoots) at the two sampling times are reported in Table S1. To correctly analyze the data, a mixed model for repeated measurements was conducted. For the AMF total colonization (TC%), no significant interaction (p-value = 0.0853) was found between "treatments" and "time", while a significant effect (p-value < 0.001) was noticed for both single factors. The pairwise multiple comparison test revealed a general significant decrease in TC% (i.e., averaged of the four AMF treatments) between T0 and T1 (data not shown). A significantly lower TC% (i.e., average of the two sampling times for each AMF treatment) was also observed in plants associated with R. irregularis as compared to R. intraradices, while the plants associated with R. clarus and R. aggregatus had intermediate values (Table 1). Similarly, for the percentage of arbuscules (AC%), no significant interaction (p-value = 0.7479) was found between "treatments" and "time", while a significant effect (p-value < 0.001) was noticed only for the single factor "time". The pairwise multiple comparison test revealed a general decrease in AC% (i.e., average of the four AMF treatments) between T0 and T1 (data not shown). Table 1. AMF-root colonization (TC% and AC%) and total fresh weight (TFW) averaged between the two-time samplings (T0 and T1) of A. officinalis associated with each AMF species (R. irregularis, R. intraradices, R. clarus and R. aggregatus) growing for 9 days in the semi-hydroponic cultivation system. The total fresh weight was measured on the same plants associated with the four AMF species at T0 and T1. No significant effect (p-value = 0.5699) was reported between "treatments" and "time", or for the single factors (p-value = 0.3570 and p = 0.0897, respectively) ( Table 1).

Metabolic Profiles and Metabolomic Analysis of A. officinalis Plants
Metabolic profiles of root and shoot samples of plants associated with four different AMF species were analyzed using multivariate analyses. The results of the unsupervised principal component analyses (PCA) highlighted the presence of two major and distinguished clusters (R. irregularis/R. intraradices on one side and R. clarus/R. aggregatus on the other side) in both the root and shoot samples (Figure 1a,b).
Shoot samples showed a higher proximity of individuals in cluster 2 (R. clarus/ R. aggregatus) as compared to the corresponding root samples, which were more widespread in the model. Both root and shoot parts of the plants showed an outlier in cluster 1, in R. irregularis and R. intraradices treatments, respectively, and were subsequently removed from the analysis.
The UHPLC-HRMS metabolic analysis of root and shoot samples associated with the four AMF species gave, in total, 201 and 265 different mass signals, respectively, during the peak picking process. Based on the Volcano-plot analysis performed, only 96 and 101 compounds passed the defined criteria (fold change > 1.5 and p-value < 0.05) in all the possible AMF combinations tested (Figures 2 and 3).

Figure 2.
Volcano-plot analysis-Identification of up-and down-regulated compounds (p-value < 0.05 and fold change > 1.5) between A. officinalis root samples associated with four AMF species (R. irregularis, R. intraradices, R. clarus and R. aggregatus) after 9 days of growth in the semi-hydroponic cultivation system.  No differentiation was observed in treated root and shoot samples belonging to the same PCA cluster. However, significant up and down-regulation of mass signals were highlighted when comparing R. irregularis and R. intraradices treatments, from cluster 1 to R. clarus and R. aggregatus treatments from cluster 2 (Figures 2 and 3). Among root and shoot samples, 22 and 35 mass signals, respectively, showed significant differences in all possible AMF-plants combinations between the two generated clusters. Thus, special attention was given and a detailed dereplication process was followed.

Identification of Primary and Secondary Metabolites Affected by AMF
The dereplication process of root and shoot extracts obtained from A. officinalis plants associated with the four different AMF species highlighted 42 differently modulated primary and secondary metabolites, strictly related to the specific AMF association.
As shown in Table 2a, six compounds were characterized as primary metabolites. Among them, compounds 1, 2, 3, 6 were tentatively identified as amino acids, while compounds 4 and 5 as organic acids. Compounds 3 and 4 (glutamic and threonic acids, respectively) were overall increased in the shoots of the AMF-treated plants. Indeed, in the R. intraradices treatment, a significant induction of compound 3 was noticed and, together with R. irregularis treatment, of compound 4 in the shoots. Similarly, compound 1 (glutamine) was up-regulated in the roots associated with R. irregularis and R. intraradices.    Among the common compounds, which were affected in both parts of the plant and felt within our selection criteria (fold change > 1.5 and p-value < 0.05), compounds 5 and 6 showed the highest increment in plants associated with R. irregularis. In particular, compound 6 showed a maximum fold change of 8.66× and 7.37× in roots and shoots, respectively, while regarding compound 5, the increment was slightly higher in the shoots (7.67×) as compared to the roots (6.85×). Conversely, compound 2 showed a 6.22× fold change in roots and a 5.87× fold change in shoots of A. officinalis associated with R. intraradices. Regarding the increment of compounds 1, 3 and 4, only those related to the roots, for compound 1 (1.7× fold change), and to the shoots, for compounds 3 and 4 (1.58× and 5.85× fold change, respectively) of the plants associated with R. intraradices, felt within our selection criteria (fold change > 1.5 and p-value < 0.05) ( Figure 4).
Syringin derivative Benzoic acid derivative  (Table 2b). Overall, thirty-six secondary metabolites were characterized, belonging to different chemi-cal classes of natural products. The majority of the identified compounds were assigned to phenylpropanoid derivatives (20) and to glycosidic triterpenes (9, saponins), while a less pronounced impact was observed in other chemical classes such as benzoic acids (1), secoiridoids (1), coumarins (2) and imidazolidines (1)  Among the C6-C3 derivatives, sixteen compounds were identified as mono-, di-, tri and tetrameric derivatives of caffeic acid (compounds 8, 9, 12, 15-17, 21, 23, 26, 28-30, 32-35) and four of them were found in the corresponding glycosidic form (compounds 9, 12, 15, 21). Besides the hydroxycinnamates, compounds 13, 24 and 31 were identified as derivatives of syringin and showed characteristic fragment ions at m/z 191, resulting from the cleavage of the glucose moiety, and fragment ions at m/z 176, 161, and 121 from the fragmentation of the remaining sinapyl alcohol [14]. In particular, compounds 13 and 24 were identified as the methyl derivatives of syringin and of syringinoside, already described in A. officinalis [14] The phenolic C 6 -C 3 derivatives represented the chemical group with the most important variation. In total, twenty phenylpropanoid derivatives were annotated. Fifteen of them were significantly increased in roots, nine in shoots, while three were affected in both roots and shoots (compounds 29, 31 and 33) (Figures 4 and 5). Compound 29 was equally affected in both parts, while compounds 31 and 33 showed a stronger accumulation in shoots and roots, respectively. In all cases, a higher impact was observed in the R. irregularis and R. intraradices treatments (Figures 4 and 5).
Among the C 6 -C 3 derivatives, sixteen compounds were identified as mono-, di-, tri and tetrameric derivatives of caffeic acid (compounds 8, 9, 12, 15-17, 21, 23, 26, 28-30, 32-35) and four of them were found in the corresponding glycosidic form (compounds 9, 12, 15, 21). Besides the hydroxycinnamates, compounds 13, 24 and 31 were identified as derivatives of syringin and showed characteristic fragment ions at m/z 191, resulting from the cleavage of the glucose moiety, and fragment ions at m/z 176, 161, and 121 from the fragmentation of the remaining sinapyl alcohol [14]. In particular, compounds 13 and 24 were identified as the methyl derivatives of syringin and of syringinoside, already described in A. officinalis [14]. This was evident by the mass difference of 14 (208,191,176,121) suggesting their structural similarities to syringin. An additional diagnostic fragment ion at m/z 384 suggested the presence of an additional acetyl group (CH 3 CO) in the structure, which was assigned to the position C-6 of glucose moiety (Table 3). This leads to the identification of compound 31 as the acetyl derivative of compound 24 which represents a previously undescribed molecule in the literature. Compound 10 showed a mass difference of 30 Da with respect to syringin (sinapyl alcohol glucoside). The presence of an additional fragment at m/z 208, characteristic to a methyl ester group in the structure, and of the ion at m/z 193 allowed us to hypothesize the presence of a dihydrosinapic acid glucoside moiety. Further investigation of the MS/MS data and by comparison with bibliographic references [26], compound 10 was identified as a glycosylated methyl ester derivative of dihydrosinapic acid. syringin. An additional diagnostic fragment ion at m/z 384 suggested the presence of an additional acetyl group (CH3CO) in the structure, which was assigned to the position C-6 of glucose moiety (Table 3). This leads to the identification of compound 31 as the acetyl derivative of compound 24 which represents a previously undescribed molecule in the literature. Compound 10 showed a mass difference of 30 Da with respect to syringin (sinapyl alcohol glucoside). The presence of an additional fragment at m/z 208, characteristic to a methyl ester group in the structure, and of the ion at m/z 193 allowed us to hypothesize the presence of a dihydrosinapic acid glucoside moiety. Further investigation of the MS/MS data and by comparison with bibliographic references [26], compound 10 was identified as a glycosylated methyl ester derivative of dihydrosinapic acid.  The above-mentioned fragments suggested the presence of two additional hydroxyl groups in the aglycone with respect to bayogenin and leads to the tentative identification of compound 27 as a diglycosidic derivative of dihydroxybayogenin [35].
An increased production of compounds 37 to 40 was also noticed in shoots of A. officinalis associated with R. irregularis and R. intraradicens. Both compounds, 37 and 40, showed a common molecular ion at m/z 1028.5187 and a common MS/MS fragmentation pattern, suggesting their structural similarity. Indeed, they both shared fragments at 779, 659, 617 and 455 (Table 3) (Table 3). Based on the above-mentioned HRMS/MS data and by comparison with previously reported data [14,38], compounds 37 and 40 were tentatively assigned as isomers of anchusoside-2 and anchusoside-7.

Discussion
The association between plants and arbuscular mycorrhizal fungi (AMF) is one of the most widespread symbioses [5]. These fungi provide numerous benefits to the host plants, especially in terms of nutritional assistance and resistance to a/biotic stresses [5,40]. A growing body of studies has reported the beneficial effects of these root symbionts on the modulation of specific biosynthetic pathways increasing/modifying the production of primary and secondary metabolites [2,8,9,14].
Herein, an untargeted metabolomic study on several AMF species belonging to the same genus was conducted under the highly-controlled semi-hydroponic cultivation system developed by Cartabia et al. [14], to evaluate the effects of four different AMF species (R. irregularis, R. intraradices, R. clarus and R. aggregatus) on the metabolome of A. officinalis.
During the experimental period, the root colonization was high for the four AMF species, even if a general significant decrease was noticed after 9 days of growth in the semi-hydroponic cultivation system (i.e., total colonization mean values varying between 85% at T0 to 66% at T1, and the percentage of arbuscules above 10% and close to 10% at T0 and T1, respectively).

Impact of AMF Species on Primary and Secondary Metabolites in Roots and Shoots of Anchusa officinalis
The major impact of AMF on A. officinalis metabolome was detected in the primary metabolism, mainly in the amino acid and organic acid content (compounds 1-6), but also in some specific secondary metabolites, derived from the phenylpropanoid (compounds  8-10, 12, 13, 15-17, 19, 21, 23, 24, 26, 28, 29-35) and the mevalonate (compounds 19, 20, 22, 25, 27, 37-40) pathways. The untargeted metabolomic approach performed on root and shoot tissues of A. officinalis evidenced forty-two compounds that fulfilled the defined threshold (fold change > 1.5 and p-value < 0.05) applied in the Volcano-plot analysis in at least one of the associations between the AMF species and A. officinalis.

Impact on Primary Metabolism
Six primary metabolites (compounds 1-6) were significantly affected in root and shoot samples in relation to the four AMF species. While compound 1 (glutamine) was produced in significantly higher amounts in the roots, compounds 2-6 were mainly accumulated in the shoots. All compounds (1-6) showed a similar accumulation increment in plants associated with R. irregularis and R. intraradices (cluster 1). R. irregularis and, especially, R. intraradices are strong elicitors of amino acids, such as aspartic (compound 2), glutamic acid (compound 3), glutamine (compound 1) and its derivative, pyroglutamic acid (compound 6), as well as of organic acids, such as threonic acid (compound 4) and malic acid (compound 5), as compared to R. clarus and R. aggregatus. Based on these results, R. clarus and R. aggregatus, belonging to cluster 2 (Figure 1), affect the primary metabolism less than the two other AMF strains belonging to cluster 1 (Figure 1).
Among the primary metabolites, compounds 2 and 3 have been reported as N precursors and donors, reflecting the AMF's ability to enhance ammonium assimilation in mycorrhized plants by the GS/GOGAT enzymatic pathway [41]. Together with compound 5, they are also important intermediate of the tricarboxylic acid cycle and they act as precursors for the synthesis of key amino acids, such as asparagine, threonine, lysine, isoleucine and glutamine (compound 1), which are building blocks for the production of macromolecules [14,41]. Pyroglutamic acid (compound 6), reported as the lactam of glutamic acid, is considered an important reservoir of glutamate [42] while compound 4 (threonic acid) is linked to ascorbic acid metabolism and catabolism, involved in anti-oxidant activities and correlated to the well-maintaining of plant fitness [43]. Threonic acid was already mentioned in previous studies as a naturally occurring constituent of shoots [44,45], and its modulation in AMF-plant symbiosis was mentioned by Schweiger et al. [11]. These results are in accordance with Cartabia et al. [14] who showed an important accumulation of the above-mentioned compounds in A. officinalis plants associated with R. irregularis.

Impact on Secondary Metabolism
Thirty-six secondary metabolites emerged from the Volcano-plot analysis of root and shoot samples, as the most affected compounds by the colonization of A. officinalis with the four different AMF species. The annotated compounds were divided into four major categories: (1) twenty C 6 -C 3 derivatives, from which four were classified as syringin derivatives and sixteen as caffeic acid derivatives; (2) one C 6 -C 2 derivative; (3) nine glycosylated triterpenoids; (4) six compounds belonging to other chemical classes, such as secoiridoids, coumarins, and imidazolidines. Among them, eight compounds were affected both in roots and shoots of the AMF-colonized plants in one, at least, symbiotic association (compounds 14, 20, 29, 31, 33, 34, 41, 42).
Phenylpropanoids were the most impacted secondary metabolites class by the AMF treatments. They were characterized as mono, di, tri or tetrameric derivatives of caffeic acid and of syringin according to their characteristic MS and MS/MS fragments (see Section 2.3). Monomers and dimers of caffeic acid derivatives share characteristic fragment ions at m/z 179, 161, 135, 121, deriving from the cleavage of a single C 6 -C 3 unit, while the tri-and tetrameric forms show additional ions at m/z 339, 295 and 185 derived from the cleavage of multiple units. Our results are in line with reported data suggesting a discernible enhancement of phenylpropanoid pathway in the roots of mycorrhized plants [2,10,14]. In addition, these compounds were already reported in A. officinalis [14,46].
All the identified C 6 -C 3 and C 6 -C 2 compounds were found in significantly higher amounts in plants associated with R. irregularis and R. intraradices (cluster 1) as compared to those associated with R. clarus and R. aggregatus (cluster 2). This suggests that R. irregularis and R. intraradices could influence A. officinalis metabolome in a similar way by activating common biosynthetic pathways. However, in a few cases, such as compounds 11 and 16 in shoots, their amount was exclusively affected (p-value < 0.05) in plants associated with R. intraradices. Besides the close response on metabolite productions by R. irregularis and R. intraradices species, minor differences can be observed in triggering specific compounds. In fact, AMF species belonging to the same cluster in the principal component analysis (PCA) affected A. officinalis plants in a similar way, however, not strictly identical. On the other hand, AMF-plant associations, which are differently clustered and present major differences in affecting A. officinalis metabolome, could equally affect the accumulation of specific secondary metabolites. This is the case of compounds 12, 23, 31 and 32 in roots, and compounds 11, 15 and 16 in shoots, which did not show a significant accumulation (p < 0.05) among the AMF species belonging to different clusters.
The AMF species also affect the production of oleanane-type saponins. In our analysis, and in contrast to the phenylpropanoids derivatives, these compounds were essentially affected/modulated in the shoots of A. officinalis. Saponins are involved in plant defense mechanisms against biotic constraints, such as pest or herbivores attack, and their content is strongly influenced by plant-AMF symbiosis [47][48][49]. Nine significantly modulated compounds from the different AMF treatments were tentatively identified. Our analysis led to the identification of four oleanane-type derivatives, which possessed two glycosidic units in their configuration (compounds 19, 22, 25 and 27) and five tri-glycosylated derivatives of bayogenin (compounds 20) and of oleanolic acid (compounds 37-40). The tri-glycosylated compound 20 and 25 presented a similar glycosylation, with the presence of a glucuronic methyl ester group, observed for the first time in the study of Cartabia et al. [14]. The ability of A. officinalis to produce saponins was already reported in previous studies [31,38,50] and their strong accumulation in shoot parts of mycorrhized A. officinalis is in line with the study by Cartabia et al. [14].
Seven compounds (20,25,27,(37)(38)(39)(40) were essentially modulated in the shoots of plants associated with R. irregularis and R. intraradices (cluster 1) as compared to the plants associated with R. clarus and R. aggregatus (cluster 2), from which six are exclusively identified in the shoots. On the other hand, Volcano-plot analysis of root parts showed the accumulation of two saponins (compounds 19 and 22). Similarly, to the phenylpropanoid derivatives, triterpenoids are mainly affected during the association of A. officinalis with R. irregularis and R. intraradices (cluster 1). However, within cluster 1, compounds 37 and 40 appeared to be more affected in the shoots of plants associated with R. intraradices  (fold change of 1.68 and 1.60, respectively). When compared to the AMF of cluster 2, R. irregularis failed to exert any significant upregulation of compound 25 in the shoots, while R. intraradices failed to induce any significant effect of compound 19 in the roots. In all the other cases, AMF belonging to cluster 1 significantly affected triterpenoids accumulation.
Special attention was given to compounds 37-40, which represent the most affected saponins of plants associated with R. irregularis and R. intraradices (fold change ranking from 1.34 to 4.55). They are all undescribed molecules in the literature, characterized by a conjugated malonyl-sugar moiety and they represent derivatives of anchusoside-7 and anchusoside-2 [14,31,38]. Our results are in line with similar conjugated structures identified from the association of different plants with R. irregularis [2,9]. Kobayashi et al. [51] reported the presence of FAS II gene, responsible for the synthesis of lipoic acid through the mitochondrial pathway in bacteria, in both R. irregularis and R. clarus from which some subunits are encoding for enzymes such as malonyl-CoA ACP transacylase. This could explain the potential ability of AMF to upregulate and synthesize malonyl conjugated compounds.
The four AMF species also promoted, differently, the accumulation of acetylated compounds. This is the case of compound 19, characterized as the acetyl derivative of anchusoside-9 (compound 22), already reported in A. officinalis [14]. The mass difference between those two compounds, equal to the presence of an additional acetyl group, as well as to the presence of the diagnostic MS/MS fragment ions at m/z 707, deduced an acetylated hexose unit at the C-21 of the structure [28]. Compound 19 represents a previously undescribed molecule in the literature. Besides saponins, two additional compounds from the phenylpropanoid pathway, both induced in plants associated with R. irregularis and R. intraradices (compounds 12 and 31) presented a similar acetyl-sugar conjugation (see Section 2.3). Compound 31 was tentatively identified as a new acetylated derivative of methylsyringin, while compound 12 was the 3-feruloyl-6 -acetyl sucrose [28]. To the best of our knowledge, these results pointed out, for the first time, the ability of specific AMF species to enhance the production of acetylated secondary metabolites.
The accumulation of methylated compounds was also reported in our analysis, essentially in the root parts of A. officinalis associated with R. irregularis and R. intraradices. Indeed, a significant increment of methylated syringin derivatives (compounds 13, 24, 31) and of the methylated phenylpropanoids, methyl dihydrosinapic acid glucoside (compound 10) and methylrosmarinic acid (compound 35), was observed. This result is consistent with our previous study [14] pointing out the methylation potential of AMF R. irregularis.
Despite the widely accepted fact that different AMF genera could affect plant metabolome differently [52], one of the main observations of the present study is that AMF species belonging to the same genus may induce similar, but not strictly identical, metabolomic responses in A. officinalis plants, without being strongly related phylogenetically. Indeed, the latest updates regarding the phylogenetic classification of under-investigated AMF strains showed that R. irregularis is phylogenetically more closely related to R. clarus than to R. intraradices [53]. Therefore, the outcome of the association in terms of, e.g., plant growth promotion and metabolites enhancement, is highly specific to the identity of the AMF symbiont [2,54,55].

Biological Material
Seeds of Anchusa officinalis L. were provided by Rühlemann's herbs and aromatic plants (Germany). They were surface-disinfected by immersion in sodium hypochlorite (8% active chloride) for 5 min and rinsed three times with sterilized (121 • C for 15 min) deionized water. The seeds were then germinated in plastic seed trays (37.5 × 23 × 6 cm) filled with a mix (w/w, 1:2) of sterilized (121

Anchusa officinalis Colonization
Two-week-old A. officinalis plants were transferred in 10 L pots containing a sterilized (121 • C for 15 min) mix of lava and perlite (w/w, 2:1). The substrate was half mixed with the AMF-inoculum substrate above (final-ratio lava:perlite w/w, 5:1). The plants were grown under the same greenhouse conditions as above.

Experimental Setup
Two-month-old plants (7 replicates per AMF treatment) were gently removed from the 10 L pots above and their root systems were rinsed with deionized water to eliminate lava and perlite debris. They were subsequently transferred to the semi-hydroponic cultivation system as detailed in Cartabia et al. [14] (Figure 6). Briefly, the plants were placed in a 500 mL plastic bottle (VWR INTERNATIONAL, Leuven, Belgium), cut at the base and with a 100 µm size pore nylon mesh (Prosep B.V.B.A., Zaventem, Belgium) glued on the top. The bottles (called containers thereafter) were used bottom-up, filled with 32 g of perlite (KNAUF GMBH, Dortmund, Germany), covered with a superficial layer of black lava rock (1-3 mm) and wrapped in aluminum foil to avoid algae development. The containers were transferred randomly in holes made in flex-foam supports and were maintained in the greenhouse set at the same conditions as described above. A 90% P-impoverished modified Hoagland solution (see [56]) was used at two different concentrations: diluted by 200× (referred to as Hoagland dil200× ) during the acclimatization phase (7 days) and diluted by 100× (referred to as Hoagland dil100× ) during the circulation phase (42 h). Both phases were presented in detail in Cartabia et al. [14]. Before starting with the circulation, initial flushing was performed. Each plant container received 200 mL of Hoagland dil100× solution, which was circulated at a velocity of 44 mL min −1 through the containers and then discarded. After this initial flushing, regular circulation was initiated and maintained at 7.5 mL min −1 for 42 h (T1). and diluted by 100× (referred to as Hoagland dil100× ) during the circulation phase (42 h). Both phases were presented in detail in Cartabia et al. [14]. Before starting with the circulation, initial flushing was performed. Each plant container received 200 mL of Hoagland dil100× solution, which was circulated at a velocity of 44 mL min −1 through the containers and then discarded. After this initial flushing, regular circulation was initiated and maintained at 7.5 mL min −1 for 42 h (T1).

Plant Harvest and AMF Roots Colonization
Total fresh weight (TFW), as well as root colonization, were assessed on the same 7 replicates per AMF treatment at the start (T0) and the end of the experiment (T1). The root colonization was evaluated by McGonigle et al. [57], on one-third of the root system (i.e., 1/6 of the root system at T0 and 1/6 at T1). The root fragments were first stained, following the method developed by Walker [58], and subsequently placed on microscope slides and covered with a 40 × 22 mm coverslip before observation under a bright field light microscope (Olympus BH2-RFCA, Japan) at ×10 magnification. Around 200 intersections were observed per plant to evaluate the total colonization (TC%) of roots (e.g., hyphae, arbuscules, and vesicles) and the percentage of arbuscules (AC%) (Figure 6). At the end of the circulation above (T1), the plants were harvested to proceed with the primary and secondary metabolites analysis on the other two-thirds of the root system (see below).

Plant Harvest and AMF Roots Colonization
Total fresh weight (TFW), as well as root colonization, were assessed on the same 7 replicates per AMF treatment at the start (T0) and the end of the experiment (T1). The root colonization was evaluated by McGonigle et al. [57], on one-third of the root system (i.e., 1/6 of the root system at T0 and 1/6 at T1). The root fragments were first stained, following the method developed by Walker [58], and subsequently placed on microscope slides and covered with a 40 × 22 mm coverslip before observation under a bright field light microscope (Olympus BH2-RFCA, Japan) at ×10 magnification. Around 200 intersections were observed per plant to evaluate the total colonization (TC%) of roots (e.g., hyphae, arbuscules, and vesicles) and the percentage of arbuscules (AC%) (Figure 6). At the end of the circulation above (T1), the plants were harvested to proceed with the primary and secondary metabolites analysis on the other two-thirds of the root system (see below).

Analysis of Primary and Secondary Metabolites in Roots and Shoots of A. officinalis 4.6.1. Samples Preparation
The remaining two-thirds of the root systems, as well as the shoot parts of each plant, were prepared according to Cartabia et al. [14]. Briefly, 20 mg of powdered and freeze-dried plant material was subjected to an exhaustive ultrasound-assisted extraction with 1 mL of EtOAc/MeOH ( In the following steps, an untargeted metabolomics workflow, including the normalization of the dataset (deconvolution, deisotoping, retention time -RT-alignment and gap-filling procedures), was developed for the detection of known and unknown compounds. All the raw data obtained from the high-resolution metabolomic profiling were uploaded in Compound Discoverer 3.2.0.421 software (Thermo Fisher Scientific, San Jose, CA, USA). Briefly, the peak alignment of the selected spectra was performed from 1 to 12 min with a mass tolerance (MT) of 5 ppm and a maximum shift of 2 min. The spectrum properties applied for the peak picking and detection of compounds were the following: S/N > 3, Min. peak intensity of 7.5 × 10 5  The structural elucidation of the metabolites of interest was performed by comparison of the chromatographic and spectrometric features of each respective peak with data from the literature. The high resolving power for both full scan experiments and the MS/MS fragments of the Q-Exactive Orbitrap analyzer in correlation to the accurate mass measurements assured the identification of the very important variables (VIP) compounds with high confidence. The suggested EC (Elemental Composition) for molecular ions and MS/MS fragments, as well as the respective RDBeq (Ring Double Bond Equivalents) further assisted the safe identification process. More prediction of compounds was performed by comparing data with in-house and online libraries, and fragment ions were correlated with spectra online databases. The pre-processed data of ESI − ionization were exported as a .csv file to a Microsoft Excel spreadsheet and manipulated accordingly for the data filtering and the multivariate statistical analysis followed.

Statistical Analysis
AMF colonization parameters (TC% and AC%) were subjected to a mixed model for repeated measurements fit by restricted maximum likelihood estimation. The normal distribution of residuals was checked for each dependent variable. The model took into account the heterogeneity of the variance (only for TC%) modeling the within-group errors variance structure with the "varIdent" matrix and assuming the different AMF treatments as a stratification variable [59]. Moreover, the repeated measurements (i.e., two sampling times conducted on the same replicates) were modeled through an autoregressive correlation structure of order 1 ("corAR1") for all the dependent variables (TC% and AC%). An ANOVA test of the model was provided and the interaction between "time" (T0 and T1) and "treatments" (R. irregularis, R. intraradices, R. clarus and R. aggregatus) was checked as well as the significance of every single factor separately. A pairwise multiple comparison test (with the Bonferroni correction) was computed to separate means (p-value < 0.05). Similarly, the total fresh weight of the plants was subjected (as a dependent variable) to the same mixed model as described above. Normal distribution of residuals and homogeneity of variance was checked. The model took into account the repeated measurements (i.e., two sampling times conducted on the same replicates) through an autoregressive correlation structure of order 1 ("corAR1"). Data analyses were performed by R [60] using the "nlme" package [61].
Multivariate analyses of HRMS data were carried out using SIMCA 14.1 software (Soft Independent Modelling of Class Analogy) (Umetric, Malmo, Sweden) to assign the discriminant metabolic changes between the different AMF treatments (R. irregularis, R. intraradices, R. clarus and R. aggregatus) after 9 days of the experiment. The interpretation of imported data was performed through a principal component analysis (PCA) and partial least squares discriminant analyses (PLS-DA) according to Pareto correlation. In addition, a permutation test with n = 100 was performed to exclude any overfitting of the aforementioned PLS-DA models (Figures S4-S7). Volcano-plot analyses were carried out using Compound Discoverer 3.2.0.421 (Thermo Fisher Scientific, San Jose, CA, USA) on the basis of filtering criteria such as p-value < 0.05 and fold change > 1.5, while the graphical representation (bar charts) of discriminant variables (targeted compounds) were generated with GraphPad Prism 7 (GraphPad Software, San Diego, CA, USA). One-way ANOVA and Tukey's test were provided in order to reveal significant differences (p-value < 0.05) in the discriminant metabolites between the four AMF treatments (R. irregularis, R. intraradices, R. clarus and R. aggregatus). Data analyses were performed by R (R Core Team, 2018) using "ggplots 2" [62] and "agricolae" [63] packages. Finally, the identification and matching of discriminant variables were performed by comparing the MS and MS/MS spectra with bibliographic data as well as with commercial and in-house libraries.

Conclusions
In the present study, we demonstrated that the association between Anchusa officinalis and different AMF species (Rhizophagus irregularis, R. intraradices, R. clarus and R. aggregatus) belonging to the same genus resulted in a different modulation of several metabolites. Based on our data, primary and secondary metabolites production was significantly affected especially in the plants associated with R. irregularis and R. intraradices. Indeed, a higher accumulation of phenolic compounds and of saponins was detected in roots and shoots of A. officinalis plants associated with these two AMF species. Additionally, an increased production of malonyl, acetyl and methyl derivatives of phenylpropanoids (e.g., 3-feruloyl-6 acetyl sucrose, methylsyringinoside, methylsyringin, 6"-acetylmethylsyringin) and of oleananetype saponins (e.g., acetylanchusoside-9, malonylanchusoside-2, malonylanchusoside-7) was observed. Among them, six compounds (acetylanchusoside-9, 6"-acetylmethylsyringin, malonylanchusoside-2, hydroxy-malonylanchusoside-7, hydroxy-dimalonylanchusoside2/7 and malonylanchusoside-7) were tentatively characterized as new secondary metabolites. Within this study, evidence leads to the AMF species-specific metabolic response of A. officinalis. However, some AMF may be more closely related to each other in modulating the plant metabolome of their host. These observations may open the door to the selection of the most adequate AMF species and/or strains for the production of desirable active compounds.