The field experiment was conducted by the Università Cattolica del Sacro Cuore in the vineyard of the farm Res Uvae in Castell’Arquato (Piacenza, northern Italy, 44°51′26.0″N 9°51′23.8″E), and grapevines (Vitis vinifera L.) of the white variety Malvasia di Candia Aromatica, grafted onto Kober 5BB rootstocks, were used. Plants were planted in 2019, spaced 1.2 m within and 2.5 m between rows, and pruned using the Guyot training system.

The characteristics in the top 0–30-cm soil layer were as follows: sand 242 g/kg, silt 356 g/kg, clay 402 g/kg, pH H₂O 7.5, organic matter 8.8 g/kg, total N 0.71 g/kg, available P (Olsen) 32.07 mg/kg, exchangeable Ca 3,819.64 g/kg, exchangeable Mg 546.91 meq/100 g, exchangeable K 168.16 mg/kg, exchangeable Na 0.27 meq/100 g, and cation exchange capacity 26.1 meq/100 g.

Standard vineyard floor management practices—detailed in Table S1 —included phytosanitary treatments according to Integrated Pest Management (IPM) principles and supported by DSS (vite.net from Horta s.r.l.) to estimate the proper timing of interventions.

The experimental design consisted of thirty-one 85-m-long adjacent rows of grapevines arranged to test mycorrhizal-based treatments (listed in Table 1 ). In detail, seven different AMF treatments were used: Micoflow + Radimix advance (T1), MycoApply DR (T2), MycoUp (T3), Overground ST-Plus (T4), Team mix (T5), Tricoveg (T6), and Vici Rhyzoteam WG (T7). A control condition (Control) was also added to compare the no-treatment effect. All treatments included four rows, except for T7, which had only three. The treatment dosage and time of application (described in Table S2 ) followed the label instructions. The biostimulants were applied once per year for three consecutive years, starting in 2019.

Weather data were collected from the Università Cattolica del Sacro Cuore weather station (PESSL®). The Datalogger collected real-time data for in-field measurement of temperature, rain, flow rate, leaf wetness, and relative humidity, among other parameters. The data for the experimental period were downloaded from a web user interface.

In the same experimental vineyard, the study of fungal ecology has been carried out with periodic sampling of roots and soil of the different treatments. The samplings were carried out every February for 3 years to determine fungal colonization. The samples were sent to the Biome Makers laboratory (Valladolid, Spain) for molecular analysis. In order to characterize fungal microbial communities, the internal transcribed spacer (ITS) marker regions were selected, and libraries were prepared following the two-step PCR Illumina protocol using custom primers amplifying the ITS1 region described previously (Becares and Fernandez, 2017). Then, they were sequenced by an Illumina MiSeq device using pair-end sequencing (Illumina®, San Diego, CA, USA). These data were used to check if the application of the different biostimulants worked and was able to increase the abundance of AMF in the roots of inoculated plants. Bioinformatic processing of the reads obtained and taxonomic annotation were performed using the Vsearch function in UCHIME and SINTAX algorithms (Edgar, n.d.; Edgar et al., 2011; Rognes et al., 2016). The operational taxonomic unit (OTU) tables generated for each year were used to evaluate the relative abundance within the fungal community from the database of sequenced species using the dplyr R package (Yarberry, 2021). Once detected, their relative abundance was evaluated over the years for each treatment in comparison to the relative abundance of the same genus population detected in Control.

The MultispeQ 2.0 device (PhotosynQ, East Lansing, MI, USA) was used for photosynthetic analysis of grapevines. Starting July 14, 2022, field sampling was conducted once a week. Four plants per treatment were randomly chosen, and a lateral shoot’s fourth and sixth fully expanded leaves were used for data acquisition. The analyzed parameters included the quantum yield of photochemical energy conversion in PSII (Phi2), the quantum yield of non-regulated non-photochemical energy loss in PSII (PhiNO), and the quantum yield of regulated non-photochemical energy loss in PSII (PhiNPQ), which were calculated by applying the equations derived by Kramer et al. (2004) as modified by Tietz et al. (2017). Accordingly, the three light-adapted parameters add up to 1 (Phi2 + PhiNPQ + PhiNO = 1) (Kramer et al., 2004).

On July 26, 2022, four replicates per each mycorrhizal treatment were collected and freeze-dried. Specifically, each replicate was represented by a pool of a side shoot’s fourth and sixth leaves. The collection date (BBCH = 83) was defined on the basis of photosynthetic data, choosing the day on which the differences between the different treatments were most marked. Samples were ground with mortar and pestle using liquid nitrogen, and the leaves from the same plant were pooled. Then, extraction was performed according to Buffagni et al. (2021). Briefly, 1 g of plant material was extracted in 10 mL of 80% (v/v) methanol acidified with 0.1% formic acid. Mechanical homogenization was performed using an Ultra-Turrax (Polytron PT, Stansstad, Switzerland); then, extracts were centrifuged (12,000 × g) for 10 min and filtered on 0.22-µm cellulose filters into glass vials. The same extraction procedure was replicated twice from two independent sub-samples from each replication.

Untargeted metabolomic analysis was performed using a 1290 ultra-high pressure liquid chromatographic system coupled with a G6550 quadrupole-time-of-flight mass spectrometry (UHPLC/QTOF-MS) from Agilent Technologies (Santa Clara, CA, USA) as previously reported (Buffagni et al., 2021). Reversed-phase chromatography was performed using an Agilent Poroshell 120 PFP column (100 mm × 2.1 mm i.d., 1.9-μm particle size) and a mobile phase of water and acetonitrile (6%–94% acetonitrile in 33 min) both acidified with formic acid (0.1% v/v). The mass analyzer operated in positive mode (ESI+), and the mass spectrometer worked in full scan mode (range 100–1,200 m/z, 35,000 full width at half maximum (FWHM) nominal resolution).

The MassHunter Profinder 10.0 software (Agilent Technologies) was used to analyze raw data for alignments and annotation by applying the “find-by-formula” algorithm. In detail, mass (with a 5-ppm tolerance) and retention time (0.05 min maximum shift) alignment were combined with monoisotopic mass, isotopic accurate spacing, and ratio for compound annotation using the database PlantCyc 15.5, same as Hawkins et al. (2021). A filter reduction was then applied to annotated features, retaining compounds detected in at least 75% of replications within at least one treatment. According to COSMOS Metabolomics Standard Initiative (Salek et al., 2015), a Level 2 confidence in annotation, corresponding to putatively annotated compounds, was achieved (Schymanski et al., 2014).

To profile phenolic compounds in grapes, five bunches from five different plants per treatment were randomly collected at maturity on August 25, 2022, and stored at −20°C until analysis. Fruits of the same treatment were homogenized together, and the resulting pool was used to select four samples (2 g). Each sample was extracted in 20 mL of 80% methanol with 0.1% formic acid, centrifuged, filtered, and analyzed as previously described for leaf metabolomics.

The Agilent Profinder 10.0 software was used for data processing, and annotation of MS compounds was accomplished according to the Phenol-Explorer 3.6 database (Rothwell et al., 2013), using the algorithm above and achieving the same COSMOS level of confidence in annotation.

Finally, a semi-quantification of the phenolic compounds annotated was carried out. Phenolic compounds were grouped into sub-classes according to Phenol-Explorer annotations, and then a representative compound for each sub-class was used for semi-quantification. To this aim, a calibration curve was prepared for each class representative, based on standards (>98% purity, from Extrasynthèse, Genay, France) solutions in methanol 80% (v/v) analyzed under the same MS conditions. Specifically, the representative standards chosen were as follows: cyanidin (for anthocyanin equivalents), catechin (flavanols), quercetin (flavonols), luteolin (flavones, flavanones, dihydrochalcones, and isoflavonoids), sesamin (lignans), resveratrol (stilbenes), ferulic acid (phenolic acids), and tyrosol (tyrosol and other low molecular weight polyphenols). The results were expressed as equivalents (E)/g dry weight (dw) starting from the cumulate abundance of each compound within a phenolic sub-class.

Yield was determined at harvest (August 25). To this aim, 25 vines per treatment were randomly chosen, and all plant bunches were weighted together to calculate the average yield component (kg). Moreover, the weight of 50 randomly selected berries per treatment was used for the weight component (g).

Technological parameters were evaluated at harvest by randomly selecting three bunches of three different plants per treatment. Bunches were put in a plastic sample bag, squeezed to obtain juice, and filtered for analysis. The total soluble sugar (TSS; °Brix) was determined by a digital refractometer (HHTEC). Must pH and titratable acidity were measured with an automatic titrator (Crison Instruments, Barcelona, Spain). Acidity was titrated with 0.1 N of NaOH to a pH 7.0 endpoint and expressed as g/L of tartaric acid equivalents, according to International Organisation of Vine and Wine (OIV) standard methods (OIV, 2021). The concentration of l-(+) tartaric and l-(−) malic were assessed using Hyperlab analyzers (Steroglass S.r.l., Perugia, Italy) by enzymatic methods (Baiano et al., 2015).

Statistical analysis of photosynthetic data, phenolic profiles, and fruits’ qualitative traits was completed on RStudio software (4.2.1 version). One-way analysis of variance (ANOVA) combined with post hoc Tukey’s honestly significant difference (HSD) test (p-value ≤ 0.05) was carried out to highlight differences among treatments.

The software Mass Profiler Professional 15.1 (Agilent Technologies) was used to elaborate metabolomics data. The abundance of compounds was Log2 transformed, normalized at the 75th percentile, and baselined against the median in Control. Fold change (FC)-based unsupervised hierarchical cluster analysis (HCA) based on Ward’s linkage method and Euclidean distance was used to naively point out metabolomic patterns. Then, ANOVA (Benjamini multiple testing correction) was performed to establish statistically significant differences compared to Control (p-value ≤ 0.05). Starting from the MS dataset, multivariate data analysis was completed in SIMCA® 17 (Sartorius, Umeå, Sweden). Therein, metabolomics data were Pareto-scaled and log-transformed, and then supervised modeling was carried out by orthogonal projection to latent structures discriminant analysis (OPLS-DA). Calculation of R2Y (goodness-of-fit) and Q2Y (goodness-of-prediction) was performed, the model was validated through CV-ANOVA (considering a p-value ≤ 0.05), permutation testing was used to exclude overfitting, and Hotelling’s T1 distance was applied to investigate potential outliers. Thereafter, the variable importance in projection (VIP) was used to select discriminant compounds, using a score threshold of ≥1.2. Discriminant markers were finally uploaded into the Pathway Tool of PlantCyc 15.1.0 software to facilitate biochemical interpretations (Hawkins et al., 2021).

After that, the AMF-driven changes in qualitative and quantitative traits of fruits were jointly evaluated using multivariate analysis of variance (MANOVA), and the Pillai trace test was used to calculate the p-value. Subsequently, canonical discriminant analysis (CDA) was conducted to identify differences among treatments and understand the relationships among the variables measured within those groups.

The presence of mycorrhiza was determined from the data obtained by the sequencing analysis of root samples. An increase in the population of these arbuscular symbiotic fungi has been observed over the years ( Figure S1 ). Despite the limitations of the approach used, our sequencing data supported the success of the biostimulant treatments over the 3 years. The data refer to general AMF colonization, not limited to the applied strains, to account for the ecological processes driving mycorrhiza community dynamics post-inoculation. In fact, it has been suggested that knowledge of the past and present soil AMF community in the target field, and the local adaptation are essential to investigate colonization processes (Basiru and Hijiri, 2022). The whole AMF population may better represent community dynamics, as observed by Bender et al. (2019) in AMF-inoculated fields with highly diverse indigenous AMF communities. This is even more relevant when trials are performed for a long period of time following inoculation (Paz et al., 2021). Irrespective of the dynamics of the AMF community, improving AMF diversity and functions by commercial inoculants can contribute to the sustainability of the agroecosystem (van der Heyde et al., 2017).

On this basis, the AMF community has been investigated through ITS sequencing. Under our experimental conditions, after 3 years of application using a commercial strain, the abundance of AMF was double in the treated samples compared to untreated plants (Control).

The maximum and average daily temperature (°C) and the daily precipitation (mm) were recorded for the entire duration of the experiment. The data produced during the whole period of investigation (July 14 to August 25, 2022) are presented in Figure 1 , together with Phi2, PhiNO, and PhiNPQ data, which were determined on July 14, 22, and 26 and on August 2, 10, 17, and 25. Despite the hot and dry climates that characterized the entire period, the highest temperatures were remarked in the first part of our time frame (July 14 to August 6). In particular, July 15 and 22 and August 5 showed the highest average and maximum values, respectively, higher than 37°C and 29°C. Contrarily, the highest precipitation levels were observed in correspondence with sudden temperature drops (July 27 and August 18, respectively, 35.8 mm and 75.6 mm).

ANOVA and Tukey’s test on chlorophyll fluorescence parameters were separately conducted for each data collection ( Figure 1 ). Interestingly, most of the differences between treatments were observed under the highest temperatures registered (July 14, 22, and 26 and August 2). The data obtained on the hottest days (July 22 and 26) showed significant differences in photosynthetic performance. On both dates, the recorded Phi2 values registered with T3 treatment (0.558 ± 0.065 on July 22 and 0.471 ± 0.141 on July 26) were statically higher than those with Control (0.299 ± 0.033 on July 22 and 0.301 ± 0.031 on July 26).

Initially, UHPLC/QTOF-MS untargeted metabolomics approach was carried out on grape leaves to obtain a comprehensive view of biochemical plant response to the different biostimulant treatments. Overall, more than 2,400 compounds were putatively annotated and subsequently listed in Table S3 , together with their respective abundances, composite mass spectra, and retention time.

First, the fold change-based unsupervised hierarchical cluster analysis was performed to point out clusters on the metabolic signatures of samples ( Figure 2A ). The resulting dendrogram showed two main clusters, the first of which formed by T3 and T6, while the second one by Control, T1, T2, T4, and T7. Within the second cluster, two sub-groups were observed, formed by Control and T7 and by T1, T2, and T4. Interestingly, T5 samples were not merged with any of these clusters, highlighting a different modulation of metabolite profile compared to the other treatments.

Subsequently, unsupervised HCA results were further elaborated by the score plot of the OPLS-DA supervised approach, where sample projection in a two-dimensional space confirmed a clear separation as a function of the treatment factor ( Figure 2B ). The OPLS-DA model was validated by CV-ANOVA (p-value = 8E−21) and by inspecting the goodness-of-fit (Q2Y = 0.76) and prediction (R2Y = 0.98) parameters. At the same time, permutation testing (N > 100) excluded model overfitting. Finally, the VIP approach allowed us to select the compounds having the highest discrimination potential (VIP score ≥1.2) in the prediction model. Therein, isoprenoids—mainly diterpenoids, sesquiterpenes, and carotenoids—and alkaloids were mostly represented. The complete list of VIP marker compounds is provided in Table S4 .

After that, ANOVA (p-value ≤ 0.05) allowed us to annotate 246 differential compounds, which appeared significantly affected by the biostimulant treatments, compared to Control ( Table S5 ). These compounds were interpreted by the PlantCyc Pathway Tool software as provided in Figure 3 . Concerning primary metabolism, lipid biosynthesis was mainly impaired by biostimulant application, showing enhanced concentrations of glycolipids and phospholipids in the leaves of each treatment. Notably, this upregulation was observed in T1 and T2, which displayed higher levels of unsaturated lipids. Among these, markedly increased contents of 1-oleoyl-2-palmitoyl-phosphatidylglycerol, 1-palmitoyl-2-linoleoyl-phosphatidylcholine, and 1-18:2-2-18:2-monogalactosyldiacylglycerol were noticed. Moreover, secondary metabolism was also affected by the treatments; in particular, N-containing metabolites, phenylpropanoids, and terpene compounds exhibited the most substantial modulation. Treatments T1, T2, and T5 elicited N-containing metabolites, while T3 and T6 mainly showed an opposite trend. Most of these compounds belonged to alkaloids like sarpagine and norajmaline. The same accumulation trend was observed for terpene biosynthesis, even if the highest accumulation of these metabolites was registered in T1-treated samples. Among others, triterpenoids (3-methyl-1,2-didehydro-2,3-dihydrosqualene, arabidiol, and 11-deoxoglycyrrhetinate), tetraterpenoids (3,4,3′,4′-tetradehydroisozeaxanthin, prolycopene, and canthaxanthin), and sesquiterpenoids were up-modulated. Differently, phenylpropanoid concentrations were markedly decreased by all treatments, except for T2 and T4, which displayed a very slight accumulation of these compounds. Notably, cinnamate, (3R,4R)-7,2′-dihydroxy-4′-methoxyisoflavanol, 2-hydroxynaringenin, and apigeninidin contents decreased following treatments. Finally, except for the T6 application, a distinctive metabolic reprogramming of the phytohormone profile was achieved in all the treated plants, showing high concentrations of brassinosteroid (6-deoxoteasterone, 3-dehydroteasterone, and campestanol), auxin ((indol-3-yl)acetate), and cytokinin (trans-zeatin N ⁶-dimethylallyladenine) compounds.

The untargeted investigation of phenolic profiles in berries highlighted 156 putatively annotated compounds, as provided in the Supplementary Material ( Table S6 ), with flavonoids and phenolic acids being the most represented (43 and 42 metabolites, respectively). In detail, flavonoids mainly included anthocyanins (cyanidin 3-O-galactoside and pelargonidin 3-O-rutinoside) and flavonols (quercetin 3-O-glucuronide and isorhamnetin 3-O-glucoside), while hydroxycinnamic included p-coumaroyl glycolic acid and gallic acid 4-O-glucoside, and protocatechuic acid 4-O-glucoside was the most represented among hydroxybenzoic acids.

Polyphenols were then quantified as phenolic class equivalents, and the resulting values were provided as mg/g (FW). Treatment-specific profiles could be observed; as reported in Table 2 , significant differences could be observed for all phenolic sub-class equivalents (p-value ≤ 0.05), except for flavanols, flavonols, and lignans. In general, small-molecular-weight phenolics and phenolic acids displayed the highest concentrations, ranging from 0.655 mg eq./g to 0.232 mg eq./g and from 0.247 mg eq./g to 0.079 mg eq./g. Overall, T6-treated samples reported the highest values for anthocyanins (0.038 ± 0.005 mg eq./g), phenolic acids (0.214 ± 0.033 mg eq./g), and stilbenes (0.024 ± 0.003 mg eq./g). On the contrary, T5-treated samples showed the highest concentrations of flavanones/flavones/isoflavonoids/dihydroflavonols (0.011 ± 0.002 mg eq./g), flavanols (0.033 ± 0.021 mg eq./g), flavonols (0.025 ± 0.009 mg eq./g), and small-molecular-weight phenolics (0.527 ± 0.001 mg eq./g). Finally, the highest content of lignans (0.073 ± 0.005 mg eq./g) was registered following the T1 application. Among phenolic profiles, it is worth emphasizing the remarkable variance between T6 and T4 for anthocyanins and phenolic acids and between T6 and T5 for flavanones/flavones/dihydrochalcones/isoflavonoids and small-molecular-weight phenolics.

Non-significant differences (p-value ≤ 0.05) could be observed for yield, pH, titratable acidity, malic acid, and °Brix among treatments ( Table S7 ), even if statistical analysis with a less restrictive threshold (p-value ≤ 0.10) on this last response profile would not reject the H0 hypothesis, having obtained a p-value of 0.09.

Contrarily, weight and tartaric acid values were differently affected after AMF-based biostimulants were applied (p-value ≤ 0.05). Concerning tartaric acid, Tukey’s test highlighted remarkable differences between T6 and T3 treatments, respectively revealing the highest (3.725 ± 0.365 g/L) and lowest (2.798 ± 0.329 g/L) average values. Significant differences in weight variables were assessed between T4 (98.608 ± 8.480 g) and T2 (68.823 ± 4.836 g).

Nevertheless, MANOVA showed a significant effect on the treatment factor by assessing the multiple dependent variables simultaneously (p-value = 9.269e⁻⁰⁵). Subsequently, CDA was performed to thoroughly examine the discriminant power of the independent variables in distinguishing between treatment groups. CDA graphical output provided a low-dimensional visualization of between-group variation and vectors reflecting the weights of the fruit qualitative and quantitative variables in a canonical space ( Figure 4 ). The output showed that a substantial amount of the variance (62%) of the between-treatment mean differences was accounted for by the first canonical variable (Can1). In comparison, the second canonical variable (Can2), orthogonal and uncorrelated to the first one, explained the remaining 18.5%. Eight colored circles, each representing the 95% confidence interval around the treatment mean, were projected on a plane that showed the largest variation between them. The clear separation between these circles indicated that the independent variables significantly contributed to the prediction of grape samples membership to the treatment group and revealed the effective power of the identified discriminant functions for classification and prediction purposes. In contrast, overlapping or closely clustered circles indicated lower discriminant power and potential similarities between treatments. Interestingly, Can1 largely separated Control from T6, T7, and T2, while Can2 mainly differentiated T1 from Control. Contrarily, T4 and Control samples were not well disjoined by canonical variables, thus indicating similar qualitative and quantitative profiles. The seven response variables—pH, °Brix, tartaric and malic acids, titratable acidity, weight, and yield—were represented by vectors, whose angles with Can1 and Can2 reflected their correlations. The relative lengths of these vectors pointed to their contribution to discrimination among treatments. Can1 was positively correlated with weight and yield, but °Brix also separated T7 and T6 from control in the opposite direction. Can2 largely reflected tartaric acid and pH, with T1 and T6 having higher sugar content and lower acidity than T2 and T3.