Seeds from the tomato cultivar, “Star 9001” was obtained from a tomato breeding program for resistance against R. solanacearum (Stark Ayres, Pty. Ltd., Bredell, South Africa)¹ and cultivated in germination mixture (Culterra, Muldersdrift, South Africa). Plants were grown under greenhouse conditions: a light/dark cycle of 12 h/12 h, with the light intensity set at 80 μmol/m²/s and the temperature regulated to between 22 and 24^°C. Plants were rotated on a daily basis to prevent any positional effects.

The csp22 peptide elicitor with sequence ATGTVKWFNETKGFGFITPDGG (Wei et al., 2018), was synthesized at ≥90% purity (GL Biochem, Shanghai, China). A stock solution of the peptide elicitor was made to 1 mg/mL in 8 mM MgSO₄ and used as a diluted sample during the inoculation procedures.

The leaves of 6 w old mature tomato plants were treated with csp22 and subsequently stained with a 3,3′-diaminobenzidine (DAB, Sigma, St. Louis, United States) solution to visualize and detect the presence of hydrogen peroxide (H₂O₂). The protocol was performed with slight modifications as previously described (Bach-Pages and Preston, 2018). Briefly, the abaxial side of the leaves were treated with 500 nM csp22 by means of pressure infiltration and incubated for 30 min. The corresponding abaxial side of the leaves were treated with 8 mM MgSO₄, which served as a negative control. Care was taken to avoid excess wounding or mechanical damage during pressure infiltration. The DAB solution (1 mg/mL in water, pH 3.8) was prepared 1 h before use. The leaves were excised from the plant and placed within a DAB solution under light at 23^°C for 8 h with constant agitation. After the incubation period, the leaves were removed and immersed in boiling 70% ethanol for 10 min. After cooling, the leaves were transferred into absolute ethanol at room temperature and left overnight. The leaves were then imaged. The visible brown polymerized precipitate in the host tissue was produced as a result of the oxidation of DAB by H₂O₂.

Leaf disks (0.4 cm²) were punched out from fully expanded leaves using a cork borer. The leaf disks were floated adaxial side up on 200 μL MilliQ water in a white 96-well microtiter plate (Nunc, Roskilde, Denmark) which was placed under light at room temperature for 24 h. After the incubation period the water from each well was completely removed and replaced with a 100 μL of a master mix solution composed of: 34 μg/mL luminol (Sigma, St. Louis, United States) and 20 μg/mL horseradish peroxidase (Sigma, St. Louis, United States) and 1 μM csp22 in 8 mM MgSO₄. Special consideration was taken to limit mechanical damage of the leaf disks during the floating disk and water removal steps. A negative control composed of the above-mentioned master mix, excluding the csp22 elicitor, was added. The negative control was supplemented with 8 mM MgSO₄ to maintain sample volumes. Luminescence was measured every 2 min for 60 min using a Synergy HT Biotek microplate reader (Biotek Instruments, Vermont, United States). The luminescence data was exported to an excel file for further analysis. To account for natural variability three leaf disks per plant were taken. In total of 24 leaf disks were used per treatment condition.

Six-week-old tomato plants were watered generously 5 h prior to elicitor inoculation to open leaf stomata and facilitate inoculation. Three plants, selected for uniform size and appearance, were reserved for each treatment/control condition prior to the treatment process. To ensure sample consistency, leaves from the fourth node branching point of the plants were selected for elicitor/control inoculation. The respective plants were treated with (500 nM) elicitor solution by pressure infiltration into the leaves using a blunt-ended syringe. Separate plants were treated with 8 mM MgSO₄ functioning as a negative control. In each instance, the entire leaf surface was supplied with elicitor/control treatment solution to minimize biological variation. It should though be noted that the fragility and complex reticulate venation inherent with tomato leaves complicate the inoculation process and that care should be taken to avoid/limit wounding or mechanical damage. After inoculation, the plants were incubated for 16, 24, and 32 h, respectively. After each incubation time the inoculated leaves were harvested from the three selected plants, snap-frozen in liquid nitrogen to quench metabolic activity and stored at −80^°C until further use. The experimental design consisted of three biological replicates (n = 3) that were created for each elicitor/control treatment at each time point (16, 24, and 32 h), thus constituting 18 biological sample groups (3 × 2 × 3) in total, covering all conditions.

The leaf tissues frozen with liquid nitrogen were pulverized with a mortar and pestle. Two grams of leaf powder were extracted with 80% methanol in a 1:10 (w/v) ratio. The samples were sonicated twice in a sonicator bath for 30 min, with the temperature controlled at 20^°C. Cell debris was pelleted with a bench top swinging-bucket centrifuge set at 5,525 × g and 5^°C for 20 min. The supernatants were evaporated to 1 mL using a rotary evaporator at 55^°C, carefully transferred into 2 mL microcentrifuge tubes and dried in a heating block overnight at 55^°C. The samples were then reconstituted in 500 μL of 50% HPLC-grade Methanol: MilliQ water solvent (1:1, v/v). The 18 samples were filtered through 0.22 μm nylon syringe filters into vials fitted with 500 μL inserts and stored at 4^°C until analyzed. Three pooled quality control (QC) samples consisting of aliquots of all samples, as well as 50% methanol blanks were included in the sample list. The QC samples were added to check sample stability, feature legitimacy, assess intensity drifts that occur during data the acquisition process, and monitor instrumental efficiency and robustness. Each of the 18 biological samples (that were prepared in triplicate as biological repeats) were analyzed in triplicate (technical repeats) on the UHPLC-MS instrument to gain precision and accuracy. Although the method described above has been widely used in scientific metabolomics literature, it should be noted that the extraction of highly polar compounds e.g., phosphates and sugars, highly non-polar compounds e.g., several lipid and sterol species, as well as volatile compounds may only partially extract or not at all. There are currently no extraction method that can recover the entire metabolome with a high level of robustness and reproducibility (Tugizimana et al., 2013). Previous literature has shown that the described method is able to recover many of the secondary metabolites of interest present within the tomato metabolome (Gómez-Romero et al., 2010; Roldan et al., 2014).

Two microliter of each sample extract was analyzed on an UHPLC-quadrupole time-of-flight (qTOF) high-definition MS system equipped with an electrospray ionization (ESI) source (Synapt G1, Waters Corporation, Manchester, United Kingdom). The analytes were separated on an Acquity HSS T3 reverse-phase column (2.1 × 150 mm × 1.7 μm; Waters Corporation, Milford, MA, United States) using a binary solvent system consisting of acetonitrile (Romil Chemistry, Cambridge, United Kingdom): MilliQ water, with both solvents containing 0.1% formic acid (FA, Sigma, Munich, Germany) and 2.5% isopropanol (Sigma, Munich, Germany). A gradient elution method was used over a 30 min run with a flow rate set to 0.40 mL/min. The elution was started at 2% (v/v) acetonitrile from 0 to 1 min, raised to 70% acetonitrile from 1 to 22 min, taken up to 95% from 22 to 23 min then kept constant at 95% acetonitrile from 23 to 26 min. The composition of the mobile phase was then reverted to 2% acetonitrile from 26 to 27 min, for column cleaning and equilibration from 27 to 30 min. The chromatographically separated metabolites were detected with the aid of a SYNAPT G1 high definition mass spectrometer (Waters Corporation, Manchester, United Kingdom) set to acquire data in both positive and negative ionization modes. The MS conditions were as follows: capillary voltage of 2.5 kV, sample cone voltage of 30 V, microchannel plate detector voltage of 1,600 V, desolvation temperature of 450^°C, source temperature of 120°C, cone gas flow of 50 L/h, desolvation gas flow of 550 L/h, m/z range of 50–1,500, scan time of 0.2 s, interscan delay of 0.02 s, mode set as centroid. The lockmass flow rate was 0.1 mL/min, with leucine encephalin as reference calibrant (50 pg/mL, [M + H]+ = 556.2771 and [M – H]^– = 554.2615), continuously sampled every 15 s, producing an average intensity of 350 counts/scan in centroid mode, with typical mass accuracies between 3 and 5 mDa and a mass accuracy window of 0.5 Da. High purity Helium was used as desolvation-, cone-, and collision gas. The MS analyses performed in an unfragmented as well as four fragmenting experiments (MS^E) simultaneously using ramping of the collision energy from 10 to 50 eV. The data acquisition at these collision energies was used to facilitate metabolite fragmentation and ease downstream structure elucidation and compound annotation. Each of the three biological replicates was analyzed in triplicate on the UHPLC-MS system, creating the technical replicates (n = 3). Thus, data for nine samples were obtained (n = 9) that was further processed by multivariate data analyses (MVDA).

The UHPLC-ESI-MS data sets were analyzed with Masslynx XS™ software (Waters Corporation, Manchester, United Kingdom), with the addition of other statistical programs for MVDA. The raw UHPLC-ESI-MS data was processed with MarkerLynx XS™ 4.2 software, with the following parameters: 0.60–21 min retention time (Rt) range of the chromatograms and m/z domain of mass range 50–1,500 Da. The Rts were allowed to differ by ±0.20 min and the m/z values by ±0.05 Da. The mass tolerance was 0.01 Da and the intensity threshold was 10 counts. Only the data matrices with noise level less than 50% (MarkerLynx cut off) were retained for downstream data analyses. The MarkerLynx application uses the patented ApexPeakTrack algorithm to perform accurate peak detection and alignment. Furthermore, MarkerLynx performs sample normalization, based on total ion intensities of each defined peak. Prior to calculating intensities, the software performs a modified Savitzky-Golay smoothing and integration (Zhou et al., 2012; Chen et al., 2013; Tugizimana et al., 2016).

The generated data matrices were imported into SIMCA (soft independent modeling of class analogy) software, version 14.0 software with the “omics” skin (Sartorius Stedim Data Analytics AB, Umeå, Sweden) for statistical analyses. To put all variables on equal footing, and adjust for measurement errors, the Pareto-scaling method was applied to data prior to chemometric modeling. A non-linear iterative partial least squares algorithm (built-in the SIMCA software) was used to handle missing values, with a correction factor of 3.0 and a default threshold of 50%. Unsupervised and supervised learning methods were applied. As part of the unsupervised methods, both principal component analysis (PCA) and hierarchical cluster analysis (HiCA) were applied. PCA is an unsupervised projection-based statistical method, used for reducing the multi-dimensionality of the data matrix. This is done by projecting the data matrix into lower dimensional space (2D) where global and qualitative visual representation of the observations can be observed. This results in the discovery and summarization of underlying clusters, trends, or sample outliers, as well as displaying the systematic variation present within the data matrix. Orthogonal projection to latent structures discriminant analysis (OPLS-DA modeling), as the chosen supervised method and binary classifier, was applied as the variable selection method to compare the samples of the elicitor treatment to that of the MgSO₄ controls—leading to identification of metabolite features positively correlated to the discrimination between the two groups (Tugizimana et al., 2013). The OPLS-DA separates multivariate relationships into predictive (positively correlated to the csp22-treatment) and orthogonal (positively correlated to the control) variation. As a tuning procedure in computing the models, a sevenfold cross-validation (CV) method was applied (Tugizimana et al., 2016). Thorough model validation steps were consistently applied; and only statistically valid models were examined and used in data mining for metabolite annotation. For variable selection, the OPLS-DA loading S-plots were used. The loading plot displays an S-shape when the data is centered and Pareto-scaled. The loadings plot facilitates the identification of features with variability between groups (discriminating variables), i.e., variables situated at the upper right or lower left sections in the S-plot. Discriminant features with a |p(corr)| of ≥0.5 and a co-variance value of |(p1)| ≥0.05 were selected for further analysis (Trygg et al., 2007). It should be noted that these selection parameters are largely data dependent. These selection parameters were chosen based on the application of descriptive statistics (ANOVA) to the downstream metabolite features, where a statistical cut-off of p-value <0.05 was used. To avoid selection bias, the statistical significance of each potential marker was investigated with the application of univariate descriptive statistics. These analyses included the generation of average intensity values, standard deviations, p-values, fold changes, and each feature’s coefficient of variation in control and treatment samples. The overall selectivity of the OPLS-DA models was assessed by constructing receiver operating characteristic (ROC) curves. The ROC curves illustrated the supervised model’s ability to discriminate between features correlated to the different sample conditions. The predictive capacity of the supervised models was assessed with a n = 100 response permutation test. The permutation test consisted of comparing the Q² obtained for the original data with the distribution of the Q²-values calculated during the randomly assigned permutations (Triba et al., 2014).

The chemical—and structural identities of the metabolites were elucidated using their respective mass spectral properties and patterns obtained during the MS analyses. MS spectral-based metabolite identification was performed based on sufficient and accurate mass fragment information, accurate calculation of each feature’s elemental composition and database searches for possible metabolite annotation. MassFragment, a built in Markerlynx software tool, was utilized for assigning possible structures to observed fragment ions of the precursor metabolite features using novel algorithms. The putative empirical formula of each statistically significant extracted ion peak (XIC) in the mass spectra was obtained and searched in databases (ChemSpider,² Dictionary of Natural Products,³ PubChem,⁴ the KEGG Compound database,⁵ and Metacyc⁶ database for the identification of possible compound matches (Gómez-Romero et al., 2010; Mhlongo et al., 2016; Zeiss et al., 2018). The analysis and identification of lipids was facilitated using the Lipidmaps database.⁷ Metabolites were annotated and tentatively identified to level 2 of the Metabolomics Standards Initiative (MSI) (Sumner et al., 2007; Spicer et al., 2017).

The DAB staining protocol was performed, as a qualitative validation method, to verify that the elicitor was perceived by the plant and that a defense-related physiological response was triggered. The left abaxial side of the leaves was treated with the elicitor solution, while the right abaxial half was supplied with the MgSO₄ control. The formation of a color product (Figure 1A) revealed the presence of H₂O₂ and associated ROS. The associated production of ROS in planta hints to the trigger of the oxidative burst and suggests the activation of initial plant immune responses. The time-dependent kinetics of the elicitor-linked oxidative burst was investigated using the chemiluminescence assay (Figures 1B,C). Univariate statistics were applied, in the form of a Student’s t-test, to determine the statistical significance (P ≤ 0.0001) between the two conditions. The data generated from the DAB stain and chemiluminescence assay confirmed the plant’s ability to perceive of the csp22 peptide elicitor. This is in contrast to the non-perception of the polymorphic flg22 from R. solanacearum, regarded as part of the pathogen’s immune evasion strategy (Wei et al., 2018, 2020).

As an analytical platform, reverse-phase UHPLC-MS is especially suitable to reflect the overall phytochemical abundance of plants, including secondary/specialized metabolites which have a defensive or protective function (Tugizimana et al., 2013). Following leaf inoculation, harvesting and metabolite extraction with 80% methanol (able to extract a wide range of semi-polar and non-polar metabolites), the sample extracts were analyzed on such a high definition/accurate mass UHPLC-MS system. The base peak intensity (BPI) chromatograms highlighted the relative peak intensities and adequate resolution of individual peaks as well as the increased ionization of S. lycopersicum metabolites in ESI negative mode (Figure 2). Qualitative variation is reflected by the peak size—where the y-axis represents the relative intensity of metabolites on the x-axis at their respective Rts (min). Changes in peak intensities and/or the presence/absence of peaks reflect differential variation in csp22-induced leaf metabolism. A comparison of the chromatograms between the control and treated conditions highlights the induced metabolic perturbations over the described incubation points in the 20 min Rt window (Figures 2A,B).

Multivariate data analysis was performed as the chosen method of data exploration by: (1) revealing trends within the metabolome across the incubation time points, and (2) detecting similarities/differences in the metabolite profiles of the treatment conditions. Processing of the chromatographic data yielded a data matrix containing the relative intensities of detected metabolite features (variables) present in each of the samples (observations). The data matrix was subjected to PCA (Figure 3A—an unsupervised projection-based statistical method, used for reducing the multi-dimensionality of the data matrix. By projecting the data embedded in the matrix into lower dimensions (2D), a global and qualitative visual representation of the observations can be observed—leading to the discovery and summarization of underlying clusters, trends, or sample outliers, thus displaying the systematic variation present within the data matrix. Additionally, the data matrices were subjected to HiCA (Figure 3B), an unsupervised hierarchical-based statistical method that was used complimentary to PCA. HiCA was performed: (1) to build a hierarchy of the sample observations and, (2) to subsequently reveal trends within the matrix that may be overlooked during PCA analysis. The HiCA model was constructed based on Ward’s linkage method, considering distance clusters between- and within-samples. The PCA scores plot (Figure 3A) revealed the partial overlap between some of the observational groups e.g., the 32 h control and treated groups. The partial overlapping pattern observation could be attributed to a single immune-inducing elicitor treatment rather than a cocktail of elicitors frequently associated with live pathogen infection. The immune response produced by csp22-treatment may leave the flux of some metabolic pathways unaffected resulting in metabolites with similar cellular levels in both the control and treated groups. This pattern of similarity in some metabolite classes may contribute to the partial overlapping pattern observed in the PCA scores plot.

In contrast to unsupervised learning methods that evaluate the patterns within data, supervised methods are designed for the prediction, classification and discovery of biomarkers (Ren et al., 2015). OPLS-DA models were constructed, as the selected supervised statistical method, to inform on class separation (Figure 4A). The statistical significance for the observed class separation in the OPLS-DA models were measured by calculating the cross-validated analysis of variance (CV-ANOVA) p-values, as a tuning method, applying a cut-off of p <0.05 (Trygg et al., 2007; Eriksson et al., 2008). The p-values of each computed supervised model were tabulated (Supplementary Table 1). The corresponding OPLS-DA loadings S-plot (Figure 4B) showed variables that were positively correlated to class separation, i.e., csp22 treatment at the selected incubation time points. These multidimensional analyses were then applied toward the identification of significant features to consider for annotation.

Discriminant features with a | p(corr)| of ≥ 0.5 and a co-variance value of |(p1)| ≥0.05 were selected for further analysis (Trygg et al., 2007). These selection parameters were chosen based on the application of descriptive statistics (ANOVA) to the downstream metabolite features (Figure 4B), where a statistical cut-off of p-value <0.05 was used. This improved the selection of statistically relevant features, while simultaneously excluding the selection of false positives. The overall performance of the OPLS-DA models, in terms of selectivity, was assessed with the construction of a ROC curve (Figure 4C), showing that the supervised models, as binary classifiers, had perfect discrimination with regards to sensitivity and specificity (Supplementary Table 1). The predictive capacity of the supervised models was assessed with a n = 100 response permutation test (Figure 4D; Eriksson et al., 2008). The permutation test revealed that the presented supervised models had higher calculated R² and Q²-values (Supplementary Table 1) in comparison to the model permutations, concluding that the OPLS-DA model obtained was statistically superior to the generated permutations (Eriksson et al., 2008; Saccenti et al., 2014). The corresponding sets of figures for the 24 and 32 h time points are presented as Supplementary Figures 1, 2.

The OPLS-DA S-plot (Figure 4B) was used for the selection of metabolite features positively correlated to the elicitor treatment. The cut-off values of |p(corr)| of ≥ 0.5 and a co-variance value of |(p1)| ≥0.05 were determined based on the application of univariate descriptive statistics (the calculation of control and treatment averaged peak intensities, the standard deviation, the coefficient of variation, fold change between the two groups, and the p-value) on the selected features. Descriptive statistics were used to provide basic information and summarize the characteristics of the variables between the samples from control and treated tissues. Metabolite features that were found to contain a p-value >0.05 and a coefficient of variation >30 were removed from the generated list of features. Using this process, the optimal values threshold values for |p(corr)| and |(p1)| were determined.

From 1,575 features (a combination of individual mass signals, including parental ions, isotopologs and possible MS fragments) initially acquired through ultra-high performance liquid chromatography coupled the mass spectrometry analyses, a total of 36 metabolites that were positively correlated to the csp22 treatment were putatively identified in the leaf tissue of S. lycopersicum (Table 1). The metabolites have been previously described in literature and databases related to tomato or plant species within the Solanaceae family. The metabolites were tabulated according to increasing Rt with corresponding m/z values. Among the identified metabolites, several of the phytochemicals belong to the phenylpropanoid class (including cinnamates, benzoates, flavonoids and coumarins), conjugates such as amides of hydroxycinnamic acids (HCAs) and chlorogenic acids (CGAs), organic acids and amino acid derivatives and lipid classes. The metabolites were annotated based on: (1) accurate mass, utilized for calculating an empirical formula, (2) analysis of each compound’s mass fragmentation, and (3) comparative analysis with existing literature.

Below, the structural elucidation of acetyl tryptophan rhamnoside (a deoxyhexoside) is shown, demonstrating the steps involved in metabolite annotation. Each metabolite’s mass fragmentation pattern was investigated in both ionization modes, in conjunction with MS fragmentation at different energies (MS^E), to illustrate how the MassFragment plugin of the MassLynx software facilitated the annotation of the fragment ions and verify the overall mass fingerprints. The elemental compositions of fragments were also calculated as a secondary method of validating each compound’s structural identity (Figure 5).

The relative peak intensities of ferulic acid derivatives (feruloyl dopamine, feruloyl quinic acid and feruloyl hexoside) were compared between the MgSO₄ control and csp22 elicitor treatments over the described incubation times (Figure 6). The data revealed that feruloyl dopamine biosynthesis (Figures 6A,B) was increased by a factor of ≥ 4 throughout the 16, 24, and 32 h incubation intervals, with the overall production of the compounds peaking at the 24 h time point. This observed result overlaps with literature relating the tomato-Pseudomonas syringae pathosystem (Zacarés et al., 2007). Feruloyl quinic acid (Figures 6C,D) exhibited its highest cellular abundance (a fold increase of ≥ 8) during the early 16 h time point which slowly decreased to a new cellular homeostasis level during the 24 and 32 h intervals. A similar trend was observed with the feruloyl hexoside derivative (Figures 6E,F). The hexoside conjugate showed the highest cellular abundance in the leaf tissue at the 16 h time point with a fold increase of ≥ 2. Following the 16 h time point, the levels of feruloyl hexoside decreased, returning to a newly established cellular homeostasis. The increased production of the feruloyl derivatives over the incubation period points to an elicitor-induced reprogramming at that occurs at the metabolome level. The observed trend indicates that the tomato leaf tissue produced most of its defense-associated compounds during the 16 h time point that decrease, either by metabolite degradation or conversion, during the 24 and 32 h incubation points.

The relative abundance of the coumarin, methylumbelliferone was compared in the elicitor treatment and the MgSO₄ control over the three incubation times (Figure 7A). The data indicates that methylumbelliferone production was increased up until the 16 h interval followed by a decline to a new cellular homeostasis. The metabolite’s synthesis increased by a factor of ≥ 8 during the 16 h incubation interval indicating an early involvement in the host immune response (Figure 7B).

Similarly, the cellular levels of the phytohormone derivative, hydroxyjasmonic acid hexoside, were also monitored between the two conditions over the incubation intervals (Figures 7C,D). The phytochemical displayed two cellular increases during the 16 h and 32 h intervals, with a drop in the overall abundance during the 24 h point (Figure 7C). It should be noted that the metabolite levels at the 16 h time points, although statistically validated, demonstrated high levels of intra-sample variability. The jasmonic acid derivative showed a ≥ 2-fold change during the 16 h interval followed by a delayed but more intense ≥ 4-fold change in the 32 h interval (Figure 7D). A similar pattern was observed in the cellular levels of coumaroyl tyramine hexoside (Figures 7E,F). This HCA amide (HCAA) produced two statistically significant increases during the 16 and 32 h time intervals, while the 24 h point showed levels overlapping with the MgSO₄ control. The compounds produced a ≥ 2-fold change during the 32 h time interval (Figure 7E).