Seeds of tomato (Solanum lycopersicum L.–Hybrid F1 Chicco Rosso) were sown in trays with size of pots of 100 ml each containing a commercial peat-based substrate (Substrate 2, Klasmann-Deilmann GmbH, Germany) having the following characteristics: density, 160 kg m³; total pore space, 85% v/v; total carbon, 55%; pH 5.5; N, 210 mg L⁻¹; P, 105 mg L⁻¹; K, 224 mg L⁻¹; and Mg, 100 mg L⁻¹; trace elements in chelated forms. Substrate was watered to water holding capacity. Trays with seeds were kept for 2 days at 4°C in the dark. Trays with seeds were placed in the controlled growth chamber (FS-WI, PSI, Czechia) at a 16-h day/8-h night regime, 22°C day/20°C night, 60% relative humidity, and with cool-white LED (250 μmol photons m⁻² s⁻¹) and far-red LED (5.5 μmol photons m⁻² s⁻¹) lighting.

Prior to plant transplanting into 3-L pots, trays were uniformly watered at 6, 7, 12, and 14 days after placement of trays in a controlled growth chamber. On day 7 and day 14, plants were fertigated with a solution containing 1.04 g L⁻¹ calcium nitrate (15.5% N; 28% CaO), 0.04 g L⁻¹ ammonium nitrate (34% N), 0.14 g L⁻¹ monopotassium phosphate (52% P₂O₅, 34% K₂O), 0.18 g L⁻¹ potassium sulfate (50% K₂O, 45% SO₃), 0.5 g L⁻¹ magnesium sulfate (10% N, 16% MgO), and 0.5 ml L⁻¹ FloraMicro (5% N, 1% K₂O, 5% Ca, 0.01% B, 0.001% Cu, 0.1% Fe, 0.05% Mn, 0.0008% Mo, and 0.015% Zn).

Twenty-day-old plants were selected with uniform growth characteristics and transplanted into 3-L pots (mixture of Substrate 2 Klasmann soil and river sand in 3:1 ratio was used). The pots were labeled with unique identification codes for each plant replicate and treatment. For determining the water content at container capacity, one set of substrate pots was dried for 3 days at 80°C and another set was saturated with water and left to drain for 1 day before weighing 100% water holding capacity (Awlia et al., 2016). Water content at container capacity was calculated as the difference between substrate weight at water holding capacity and dried substrate. On the day before transplantation, soil was prepared, and moisture content was adjusted to 60% of container capacity. Twenty-one-day-old plants were transplanted into the prepared substrate mixture with 60% of container capacity. Following the transplantation, plants were regularly watered to reference weight (40% of container capacity) defined as low water availability condition by using the automated watering and weighing unit of the PlantScreen^TM Modular System (Photon Systems Instruments (PSI), Czechia).

Plant-derived PH biostimulant Trainer^® was provided by Italpollina Company (Rivoli Veronese, Italy). The plant-derived PH Trainer^® is a commercial PH obtained through enzymatic hydrolysis of proteins derived from legume seeds. Briefly, the seeds are ground, and the flour was dispersed in acidified water to extract the soluble compounds. Filtration and centrifugation are then used to separate the protein concentrate from the other organic compounds. Enzymatic hydrolysis is used to release the amino acids and peptides from protein concentrate. Insoluble residual compounds are separated from amino acids and peptides by centrifugation. The resulting PH is concentrated through water evaporation (Colantoni et al., 2017). The final product contains mostly peptides and amino acids and, with a less extent, soluble carbohydrates, mineral elements and phenolic compounds. Trainer^® has a density of 1.21 kg L⁻¹, a dry matter of 46%, and a pH of 4.0. It contains 310 g kg⁻¹ of free amino acids and soluble peptides (Rouphael et al., 2018a). The aminogram of the product (in g kg⁻¹) was as follows: Ala (12), Arg (19), Asp (33), Cys (4), Glu (54), Gly (13), His (8), Ile (12), Leu (24), Lys (19), Met (4), Phe (16), Pro (15), Ser (17), Thr (11), Trp (4), Tyr (13), and Val (16). The antioxidant activity of Trainer^®, as measured by ferric-reducing antioxidant power (FRAP), was 41.9 mmol Fe²⁺ g⁻¹ f.w., while the total phenolics and flavonoids, determined following the methods reported by Borgognone et al. (2014), were 8.93 mg of gallic acid equivalent per gram of f.w. product and 0.95 mg of quercetin equivalent per gram of f.w. product, respectively. The Trainer^® content of soluble sugars was 90 g kg⁻¹ f.w., and its elemental composition was as follows (g kg⁻¹ f.w): N (50.0), P (0.9), K (41.1), Ca (10.9), Mg (0.5), Fe (0.024), Zn (0.010), Mn (0.001), B (0.005), and Cu (0.001) (Colla et al., 2017a). The Trainer^® content of N–NO₃ and N–NH₄ was 3.13 and 6.00 μg g⁻¹ f.w., respectively (Ceccarelli, 2018). No detectable phytohormones in Trainer^® have been reported (Luziatelli et al., 2019).

Plants were randomly distributed into three groups with six biological replicates per group. Three groups each containing six plants were identified as follows: no application, foliar application, and drench application of PH. Each plant was labeled with a unique barcode identifier used for registration of the plants in the PlantScreen^TM Modular System.

The PH was applied either as foliar spray or as substrate drench (Supplementary Figure S1B) as water solution containing a non-ionic surfactant Triton X-100 at 0.1%. A control group (no application) was sprayed with distilled water containing 0.1% Triton X-100. PH application was performed twice: 5 days after transplanting (DAT) referred to as Treatment 1 (T1) and 12 DAT referred to as Treatment 2 (T2). For 24 h prior to and following spraying, humidity in the cultivation chamber was kept at 85% relative humidity. For foliar spray treatments, 2 ml of PH was diluted in 500 ml of distilled water with 0.1% Triton X-100, and 60 ml of solution was applied by homogeneous foliar spray over the entire plant surface per plant replica. Substrate of each pot was covered with aluminum foil during and upon spraying and was removed prior to the next phenotypical analysis in the PlantScreen^TM Modular System. For drenching treatment, 4 ml of biostimulant was diluted in 1,000 ml of 0.1% Triton, and 60 ml per plant replicate was applied by drenching. At both PH application times (T1 and T2), plants in control treatment and those foliarly sprayed with PH were irrigated with 60 ml of water each to avoid changes of substrate water status in comparison with plants treated by drench application of PH. Right after PH treatment, plants were taken back to fytoscope FS-WI.

Plant phenotypic measurements were performed using the PlantScreen^TM Modular System installed in semi-controlled greenhouse environment conditions in the PSI Research Center (PSI, Drásov, Czechia). The platform was operated in closed imaging loop located in a climatized environment with temperature ranging between 21°C and 24°C. The platform is equipped with four robotic-assisted imaging units, an automatic height measuring light curtain unit, an acclimation tunnel, and a weighing and watering unit. Plants placed in individual transportation disks were transported by moving belt toward individual imaging units and watering and weighing stations.

Twenty-two-day-old plants were randomly distributed into three batches, each batch containing 12 plants. Plant imaging started with 22-day-old plants (1 DAT, day 1 of phenotyping) and continued for 15 days (15 DAT, day 15 of phenotyping). Plants were imaged using the following protocol. Briefly, plants were manually transferred from the climate-controlled growth chamber to the manual loading station of the PlantScreen^TM Modular System and were transported through the acclimation tunnel with automatic height measuring unit. Prior to the imaging, plants were dark-adapted in acclimation tunnel for 15 min. Each batch of plants was automatically phenotyped for around 30 min by using kinetic chlorophyll fluorescence imaging measurement for photosynthetic performance analysis; top view and multiple-angle side view Red Green Blue (RGB) imaging for morphological, growth, and color analysis; and finally a thermal imaging unit for plant surface temperature quantification (Supplementary Figure S1A). Following the imaging, plants were automatically transported to the watering and weighting unit for maintaining precise soil water holding capacity. After completion of the phenotyping protocol, plants were manually moved back to the climate-controlled growth chamber until the subsequent phenotyping day. We used the automatic timing function of the PlantScreen^TM Scheduler (PSI, Czechia) to schedule the initiation of the phenotyping protocol at the same time of the diurnal cycle (after 3 h of illumination in the climate-controlled growth chamber). The phenotyping data were acquired twice prior to biostimulant application in days 1 and 3 (pre-T measurements), three times post-first biostimulant application in days 6, 8, and 10 (post-T1 application), and twice post-second biostimulant application in days 13 and 15 (post-T2 application). The acquired images were automatically processed using Plant Data Analyzer (PSI, Czechia), and the raw data exported into CSV files were provided as input for further analysis.

Kinetic chlorophyll fluorescence (ChlF) measurements were acquired using an enhanced version of the FluorCam FC-800MF pulse amplitude modulated (PAM) chlorophyll fluorometer (PSI, Czechia) with an imaging area in top view position of 800 × 800 mm, as described in Tschiersch et al. (2017). We assessed the photosynthetic performance in the plants by quantifying the rate of photosynthesis at different photon irradiances using the light curve protocol (Henley, 1993; Rascher et al., 2000). The measuring protocol described previously (Awlia et al., 2016) was optimized for the tomato plants from early to later developmental stage. For the light curve characterization, three actinic light irradiances (Lss1–170 μmol photons m⁻² s⁻¹, Lss2–620 μmol photons m⁻² s⁻¹, and Lss3–1,070 μmol photons m⁻² s⁻¹) were used with a duration of 30 s in order to quantify the rate of photosynthesis.

From the fluorescence data, a range of parameters were extracted as described in detail by Awlia et al. (2016). Additionally, 1 - q_P was calculated, which reflects the proportion of PSII reaction centers that are closed (Maxwell and Johnson, 2000; Na et al., 2014).

To assess digital biomass of the plants, RGB imaging was done from top view (RGB2) and side view from multiple angles (RGB1). The RGB imaging unit is a light-isolated box equipped with turning table with precise angle positioning and two RGB cameras (top and side) mounted on robotic arms, each supplemented with LED-based lighting source to ensure homogeneous illumination of the imaged object.

Projected shoot area (PSA) parameter, together with regularly determined weight of the plants, was used to estimate TUE. TUE was defined by the ratio of aboveground biomass produced per unit of water transpired and depends on the characteristics of the plants and on the environment where the plants grow (Al-Tamimi et al., 2016). TUE was estimated from transpiration defined by measures of water loss and growth from PSA by plant-specific pixel counts quantification.

To assess leaf surface temperature of the plants, a thermal imaging unit based on side view imaging was used. The thermal imaging unit incorporated in the PlantScreen^TM System consists of a light-isolated box with one side view camera mounted on a robotic arm, precise plant positioning, and a background heated wall with an integrated temperature sensor to increase contrast for the image processing step. The imaged area is 1,205 × 1,005 mm (height × width). To assess spatio-temporal variations in temperature over plant surface, we used FLIR A615 thermal camera with 45° lens and resolution 640 × 710 pixels, with high-speed infrared windowing option and <50 mK thermal sensitivity (FLIR Systems Inc., Boston, MA, United States). The thermal images were acquired in line scan mode with each image consisting of 710 pixels with a scanning speed of 50 Hz (lines per second). Thermal images were acquired in darkness. Image acquisition conditions, plant positioning, and camera settings were fixed throughout the experiment. Leaf surface temperature of each plant was automatically extracted with Plant Data Analyzer software (PSI, Czechia) by mask application, background subtraction, and pixel-by-pixel integration of values across the entire plant surface area. To minimize the influence of the environmental variability and the difference in the image acquisition timing among individual plants, the raw temperature of each plant (°C) was normalized by the actual background temperature and expressed as ΔT (°C).

Plant material was harvested 19 DAT for metabolomic analysis by harvesting and combining the third and fourth fully expanded leaves from the top of each plant. Additionally, the final biomass of each plant was determined by measuring fresh weight and dry weight of the remaining shoot.

Plant samples were homogenized in pestle and mortar using liquid nitrogen, and then an aliquot (1.0 g) was extracted in 10 ml of 0.1% HCOOH in 80% aqueous methanol using an Ultra-Turrax (Ika T-25, Staufen, Germany) (Borgognone et al., 2016). The extracts were centrifuged (12,000 × g) and filtered into amber vials through a 0.22-μm cellulose membrane for analysis. Thereafter, metabolomic analysis was carried out through a ultra-high performance liquid chromatograph (UHPLC) coupled to a quadrupole-time-of-flight mass spectrometer (UHPLC/QTOF-MS). The metabolomic facility included a 1290 ultra-high-performance liquid chromatograph, a G6550 iFunnel Q-TOF mass spectrometer, and a JetStream Dual Electrospray ionization source (all from Agilent Technologies, Santa Clara, CA, United States). The untargeted analysis was carried out as previously described (Rouphael et al., 2016). Briefly, reverse-phase chromatography was carried out on an Agilent Zorbax Eclipse-plus C18 column (100 × 2.1 mm, 1.8 μm) and using a 34-min linear elution gradient (5% to 95% methanol in water, with a flow of 220 μL min⁻¹ at 35°C). The mass spectrometric acquisition was done in SCAN (100–1,000 m/z) and positive polarity (Pretali et al., 2016).

Features deconvolution and post-acquisition processing were done in Agilent Profinder B.06. Mass and retention time alignment followed by a filter-by-frequency postprocessing filter were done to retain only those compounds that were present in >75% of replications within at least one treatment. Compound annotation was done using the “find-by-formula” algorithm, i.e., using monoisotopic accurate mass, isotope spacing, and isotope ratio, with a mass accuracy tolerance of <5 ppm. The database PlantCyc 12.5 (Plant Metabolic Network¹) was used for annotation purposes. Based on the strategy adopted, identification was carried out according to Level 2 (putatively annotated compounds) of the COSMOS Metabolomics Standards Initiative². The classification of differential compounds into biochemical classes was carried out following PubChem (NCBI³) and PlantCyc information.

For automatic image data processing, we used the data processing pipeline Plant Data Analyzer, which includes preprocessing, segmentation, feature extraction, and postprocessing of acquired images. Values for projected shoot area were calculated from images taken in the visible light spectrum and correspond to plant volume estimation. The plant volume was used as a proxy for the estimated biomass of the plants. Data were processed using MVApp application. Statistical differences between treatments and time points were determined by one-way analysis of variance (ANOVA) with post hoc Tukey’s Honest Significant Difference (HSD) test (p-value < 0.05) performed using appropriate scripts in MVApp tool. Data are displayed as mean ± standard error of the six independent plants per treatment.

Elaboration of metabolomic data was carried out using Mass Profiler Professional B.12.06 as previously described (Salehi et al., 2018). Briefly, compounds’ abundance was Log2 transformed and normalized at the 75th percentile and then baselined against the median. Unsupervised hierarchical cluster analysis was carried out using the fold-change-based heatmap, setting similarity measure as “Euclidean” and “Wards” linkage rule. Thereafter, the dataset was exported into SIMCA 13 (Umetrics, Malmö, Sweden), Pareto-scaled, and elaborated for Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA). This latter supervised statistic allowed the separation of variance into predictive and orthogonal (i.e., ascribable to technical and biological variation) components. Outliers were excluded using Hotelling’s T2 and adopting 95 and 99% confidence limits, for suspect and strong outliers, respectively. Model cross-validation was done through CV-ANOVA (p < 0.01), and permutation testing (N = 300) was used to exclude overfitting. Model parameters (goodness-of-fit R²Y and goodness-of-prediction Q²Y) were also produced. Finally, Variable Importance in Projection (VIP) analysis was used to select the metabolites having the highest discrimination potential. A subsequent fold-change analysis and two-way ANOVA were finally performed from VIPs to identify extent and direction of the changes in accumulation related to the use of the biostimulants.

Chemical Similarity Enrichment Analysis (Barupal and Fiehn, 2017) was finally performed on VIP metabolites to critically highlight the chemical nature of the discriminant compounds, as previously described (Showalter et al., 2018). Such enrichment analysis is based on chemical similarities and used Tanimoto substructure chemical similarity coefficients to cluster metabolites into non-overlapping chemical groups. In our elaborations, OPLS-DA VIP scores were used instead of individual p-values, and the regulation (up- or down-accumulation) of discriminant metabolites was compared across treatments following chemical enrichment. The online web-app tool⁴ was used for this analysis.

Integrative phenotyping facilities provide an opportunity to combine various methods of automated, simultaneous, non-destructive analyses for assessment of plant growth, morphology, and physiology. Here, we used the PlantScreen^TM Modular System (PSI, Czechia) available in the PSI Research Center (Drásov, Czechia) for simultaneous analysis of multiple morpho-physiological traits in tomato plants treated with plant-derived PH biostimulant substances (Supplementary Figure S1A). Tomato plants were cultivated under control conditions and were phenotyped by using RGB imaging to capture plant growth dynamics, morphology and color, by chlorophyll fluorescence (ChlF) imaging to quantify photosynthetic performance and by thermal imaging to analyze leaf surface temperature prior to and following the PH treatment (Figure 1). Finally, an automated watering and weighing unit was used to maintain constant low water availability conditions in the tomato plants treated with PH by both drenching and spraying applications (Supplementary Figure S1B).

Visible RGB digital color imaging was used for the assessment of range of visual traits in control plants (no application) and plants treated with PH by either drenching (drench application) or spraying application (foliar application) (Figure 1A,B). RGB imaging was used to quantify the effect of the PH on growth status, biomass accumulation, and color of tomato plants cultivated under limited water availability conditions (Figure 2A). Simple image stacks acquired from top view and two side view images were used to extract and calculate shoot volume as a proxy of shoot digital biomass and quantify shoot color throughout the cultivation period. The morphological traits were assed dynamically and were used to calculate growth rates (Figure 2B).

The analysis of the growth-related above-mentioned traits revealed that tomato plants cultivated under low water availability conditions and treated with PH by either spraying or drenching grew better than control plants. The best-performing plants treated with PH were those where PH was applied as substrate drench. At the end of the phenotyping period, the digital shoot biomass was significantly increased (Figure 2A and Supplementary Tables S1–S3) as well as the height and width of the plants (Supplementary Tables S4, S5). In addition, the growth rate calculated over the entire phenotyping period was also strongly enhanced in drench treated plants compared to foliarly sprayed ones under limited water availability (Figure 2B), suggesting that overall growth performance of the plants was improved following the drenching application of PH. The image-based data could be further confirmed by destructive plant biomass assessment as both fresh and dry weights of the PH-treated plants harvested at the end of the experiment were increased (Supplementary Figure S2A). Measurements of projected shoot area obtained using HTP imaging approach were strongly correlated with fresh and dry weights of the plants, and there was no indication of any deviation from a linear relationship even at the highest biomasses measured in this experiment (Supplementary Figures S2B,C).

The variation in shoot color of the tomato plants over the phenotyping period was assessed by quantification of greenness hue abundance from the color-segmented RGB images (Supplementary Figure S3). The analysis algorithms were calibrated by using RGB images from all treatments and all measurements as described previously (Awlia et al., 2016). Some minor changes were observed in the analyzed green hues, but no clear trend could be observed except for the slight increase in darker green hues at the end of the phenotyping period for the drench application variant (Supplementary Table S6).

To assess the effect of PH application on photosynthetic performance of tomato plants under water-limiting conditions, chlorophyll fluorescence measurements were acquired using automated chlorophyll fluorescence imaging setup (Figure 1C and Supplementary Figure S1). The rate of photosynthesis at different photon irradiances was quantified using the light curve protocol reported by Henley (1993) and Rascher et al. (2000). From the measured fluorescence transient states, the basic ChlF parameters were derived (i.e., F_o, F_m, F_t, and F_v), which were used to calculate a range of parameters characterizing plant photosynthetic performance (i.e., F_v/F_m, NPQ, q_P, and ΦPSII) [for an overview, refer to Paul et al. (2011); Awlia et al. (2016); Tschiersch et al. (2017)]. In addition, photochemical quenching (1 - q_P) and photosynthetic electron transport rate (ETR) parameters were calculated, which refer to proportion of closed PSII reaction centers (Maxwell and Johnson, 2000) and ETR of photosystem II and indicate the efficiency of linear electron flow route in the photosynthetic machinery for producing energy-rich molecules adenosine triphosphate (ATP) and the reduced form of nicotinamide adenine dinucleotide phosphate (NADPH), respectively.

A few of the parameters were selected to dynamically characterize the photosynthetic function of PSII in the tomato plants prior to and after the biostimulant treatment under limited water availability (Figure 3): the maximum quantum yield of PSII photochemistry in the dark-adapted state (F_v/F_m), the photochemical quenching coefficient that estimates the fraction of open PSII reaction centers (q_P), steady-state non-photochemical quenching (NPQ), and ETR correlating to the quantum yield of the CO₂ assimilation mechanisms and to the overall photosynthetic capacity of the plants (Genty et al., 1989). No significant changes of those parameters between the control and PH-treated plants (Figure 3 and Supplementary Table S7) were recorded during the phenotyping period. However, minor dynamic changes in lower actinic irradiance of the 1 - q_P parameter were observed at the end of the phenotyping period on day 15 (Supplementary Figure S4). 1 - q_P was used as an indicator of the closed PSII reaction center and as an estimate of the relative PSII excitation pressure to which an organism is exposed (Maxwell and Johnson, 2000), suggesting that PH application induced a higher redox status than control treatment, resulting in slightly lowered ETRs (Supplementary Figure S4).

Plant water status is determined by the equilibrium between root water uptake and shoot transpiration (Berger et al., 2010). Under limited water availability in tomato plants, triggering of shoot transpiration and root respiration has been carried out by commercial PH provided to the plant by foliar and drenching application, respectively. Imaging thermography approach was used to measure the whole plant temperature in an automated manner, and the image data were utilized to assess the leaf transpiration of plants (Figure 1D).

To minimize the influence of the environmental variability and the difference in the image acquisition timing among individual plants, the raw temperature of each plant (°C) was normalized by the actual background temperature and expressed as ΔT (°C) (Paul et al., 2016). Experimental data showed that leaf surface temperature of the tomato plants was not influenced by PH treatment, and no difference compared to control plants was observed throughout the entire phenotyping period (Supplementary Figure S5A). In addition to leaf surface temperature we assessed TUE that increased in drenching PH-treated plants in comparison with foliar and control treatments (Supplementary Figure S5B).

A strong correlation was reported between plant transpiration rate and stomatal conductance (Berger et al., 2010). As stomatal conductance is the measure of the CO₂ entering or leaving the stomata of a leaf, higher TUE observed in PH-drench application might suggest that more CO₂ might get fixed and generate more organic matter, thereby increasing in biomass compared to other treatment methods.

An untargeted UHPLC/QTOF-MS metabolomic analysis was carried out to elucidate the molecular mechanisms underlying the effect of PH application on leaves of tomato plants grown under limited water availability. Multivariate statistics from the metabolomic dataset pointed out similarities/dissimilarities among phytochemical profiles. The use of an untargeted profiling followed by annotation on the basis of a comprehensive database (namely, PlantCyc) produced over 1,900 compounds annotated, overall. These compounds exhibited a large chemical diversity and included metabolites from a wide range of biochemical classes and metabolic processes.

The first step of interpretation was a hierarchical clustering, produced from the fold-change-based heatmap according to Euclidean distances. This unsupervised clustering approach allowed describing similarities/dissimilarities among treatments, as shown in Figure 4. As provided, two main clusters were generated–one comprising drench application and the other including foliar application and control. In this latter cluster, two distinct subclusters could be identified, thus indicating different metabolic profiles between foliar application of the biostimulant and control plants. Even though the application of PHs resulted in distinctive profiles in tomato under limited water availability, the naive (unsupervised) hierarchical clustering of metabolomic signatures suggested that the application method of the PH was an additional and relevant factor determining the actual difference in such phytochemical profiles.

A consistent outcome could be produced through the supervised OPLS-DA multivariate modeling. This analysis allowed separating predictive and orthogonal components (i.e., those components ascribable to technical and biological variation) of variance. Therefore, OPLS-DA effectively discriminated among the three groups into the score plot hyperspace. The OPLS-DA score plot (Figure 5) indicated a complete separation among control, foliar, and drench applications. The model parameters of the OPLS-DA regression were excellent, being R²Y = 0.99 and Q²Y = 0.94, respectively. The model was validated (CV-ANOVA p = 2.47 × 10⁻¹⁰) and overfitting could be excluded through permutation testing (N = 100). Validation through a misclassification table indicated a 100% model accuracy (Fisher’s probability 3.5 × 10⁻⁷). Furthermore, Hotelling’s T2 allowed us to exclude suspect and strong outliers. Given the validated model outcomes, the variable selection method called VIP (Variable Importance in Projection) was used to identify compounds explaining the differences observed. The discriminating compounds having a VIP score >1.25 were exported and subjected to fold-change analysis to identify the trends of regulation altered by the treatments. Thereafter, one-way ANOVA (Tukey post hoc) was used to describe significance of the differences. The discriminant compounds, together with their VIP score, P, and fold-change values, were grouped into chemical classes to facilitate the discussion of results (Table 1).

Notably, relatively few biochemical classes included most of the discriminant metabolites. In more detail, lipids (including membrane lipids, sterols, carotenoids, and other terpenes) were the most represented class of compounds among VIP discriminants, followed by phytohormones, polyamine conjugates, prenyl quinones, and chlorophyll-related compounds. Among hormones, brassinosteroids, indole conjugates, salicylate, cytokinins, and two gibberellins were identified among discriminant compounds of treatments (Table 1). Furthermore, abietane diterpene resin acids, as well as pteridins and few other compounds, could be outlined by VIP analysis. Interestingly, two osmolytes (trehalose and glycine betaine) were identified among VIP discriminants (Table 1).

The following chemical enrichment analysis carried out in chemRICH highlighted sterols (cholestanes, cholestadienols, and hydroxycholesterols), carotenoids, unsaturated fatty acids and phosphatidic acids, terpenes, and coproporphyrins as the most represented chemical groups (Supplementary Figure S6). The analysis, carried out separately for each application method (foliar or drench as compared to control), represented differences in accumulation for the selected metabolites. Most of the classes reported exhibited a down-accumulation following biostimulants treatment, as compared to control, except for terpenes (foliar application treatment) and unsaturated fatty acids (drench application treatment).