## 2 Materials and Methods
### VOCs Data  
The VOCs dataset was elaborated by using the R statistical software 3.5 [38].  
Differences among treatments, time of exposure and TrxTi interaction were inferred through PERMANOVA multivariate analysis (999 permutations) using the package _vegan_. Pairwise comparisons were calculated using a custom script and correcting \(p\) values using the False Discovery Rate (FDR) method.  
In order to identify VOCs key predictors that could constitute a molecular signature identifier among the treatments within each time of exposure, we used a preliminary unsupervised (Principal Component Analysis, PCA) and then supervised analysis (Sparse Projection to Latent Structure-Discriminant Analysis, sPLS-DA) by using the _mixOmics_ package [39]. Statistical algorithms are detailed in [PERSON] et al. [39] and they account for multiple comparisons inherent in biomarker datasets, where multiple classification features are considered for a relatively small number of specimens (\(p>>\) n). In particular, the sPLS-DA procedure constructs artificial latent components of the predicted dataset (VOCs Table denoted X (N \(\times\) P)) and the response variable (denoted Y with categorical information of samples, e.g., CTR, ABIO, BIO and COMB). To predict the number of latent components (associated loading vectors) and the number of discriminants, for sPLS-DA, we used the _perf.plsda()_ and _tune.splsda()_ functions, respectively. We finetuned the model using five-fold cross-validation repeated 10 times to estimate the classification error rates employing two metrics, overall error rates and balanced error rates (BER), between the predicted latent variables with the centroid of the class labels (categories considered in this study) and specifying the _max.dist_ (which gave the minimal classification rate in this study).  
Calculation of Additive, Synergistic or Antagonistic Effects in Combined Stress  
To determine whether stress combination yielded additive, synergistic or antagonistic effects respect to each of the single stress alone, we used the method of [PERSON] et al. [12]. To such aim, we compared the observed effects (Ob) with expected additive effects (Ex) for the plants exposed to the abiotic stress and herbivory combination (COMB) at 3 and 8 days of treatments, only. The Ob effect sizes were calculated as the absolute value of:  
\[\mathrm{Ob}=(\mathrm{ob}-\overline{\mathrm{x}CTR})/\overline{\mathrm{x}CTR} \tag{2}\]  
where Ob is the value of each measured trait in each plants and treatment and \(\overline{\mathrm{x}CTR}\) is the mean value of the same trait measured in CTR plants.  
For each of the traits considered, the Ex additive effect sizes for the COMB treatment were defined in two steps by first determining and then summing the independent effects (In) of each treatment. The In effect sizes were calculated as the absolute value of:  
\[\mathrm{Ind}=(\overline{\mathrm{x}stress}-\overline{\mathrm{x}CTR})/\overline {\mathrm{x}CTR} \tag{3}\]  
where \(\overline{\mathrm{x}stress}\) is the mean values of a given trait in the presence of a single stress, and \(\overline{\mathrm{x}CTR}\) is the corresponding mean value in CTR plants. Then, the Ex additive effect size for the COMB treatment was calculated by using a multiplicative risk model [40], that is the sum of the two In effects minus their product. Finally, the Ex additive values for COMB plants were compared to the actual Ob additive effects. In particular, we calculated a mean difference (\(\pm 95\%\) confidence interval) between the effect sizes of Ob and Ex was for COMB plants. When Ob-Ex > 0 and the lower 95% confidence limit was greater than zero, then the impact from the combination of both stressors was classified as synergistic. Antagonistic effects were defined when the Ob-Ex < 0 and the upper 95% confidence limit was less than zero. Finally, we classified additive effects when the 95% confidence interval crossed the zero line.