Statistical Analysis
All the experiments were carried out in a completely randomised design with 5 replications. The univariate analysis of morphological and physiological parameters was carried out using XLSTAT 2014.5.03. Data were analysed through one-way ANOVA using Duncan’s test as post hoc (P ≤ 0.05).
Concerning metabolomics, the post-acquisition data analysis was carried out using the software Mass Profiler Professional 12.6 (Agilent Technologies); the compounds were log2 transformed, normalised at the 75th percentile, and baselined against the median. Afterwards, both unsupervised and supervised multivariate statistics were applied for interpretations. According to Euclidean distance and Ward’s linkage, the unsupervised hierarchical cluster analysis was used to underline the relatedness across the different treatments. In addition, the supervised orthogonal projection to latent structures discriminant analysis (OPLS-DA) was carried out, and the model parameter (goodness-of-fit R2Y and goodness-of-prediction Q2Y) were recorded. Also, the OPLS-DA model was cross-validated (CV-ANOVA), inspected for outliers (Hotelling’s T2), and the overfitting was excluded through a permutation test (n = 200). Then, the Variable Importance in Projection (VIP) analysis was used to select the metabolites having the highest discriminant potential (VIP score > 1.2). Finally, the differential compounds obtained from the ANOVA and fold-change analysis (FC) ( p < 0.05, Bonferroni multiple testing correction and Fold-Change FC ≥ 2) were exported into the Omic Viewer Pathway Tool of PlantCyc (Plant Metabolic Network) (Stanford, CA, USA) software for interpretation (Karp et al.2009; Caspi, Dreher & Karp 2013)
Regarding metabarcoding data, downstream analysis was performed using RStudio with R version 4.1.1. Phyloseq v1.38.0 was used to handle ASV sequences and tables. Samples were split into compartments (soil, root) and separately analysed. After outlier removal (one sample from control treatment on root and in soil, one sample from drought treatment and two from Drought-Heat treatment due to small sampling size), ASVs were filtered for mitochondria, and unassigned sequences were removed. Only ASV that (a) occurred in at least 3 samples and (b) occurred >10 times in total were retained for downstream analysis, leading to 1455 ASV in soil samples and 915 ASV in root samples. Alpha diversity indices (number of observed ASVs, Inverse Simpson index) were calculated using rarefied and filtered samples from root and soil datasets and plotted by treatment. Bray-Curtis dissimilarity indices were calculated on rarefied relative abundances and used to perform principal coordinate analysis (PCoA) and permutational analysis of variance (PERMANOVA) to investigate treatments’ effect on the bacterial community structure. For PERMANOVA, 999 permutations were carried out per dataset. Linear discriminant analysis of effect size (LEfSe) was applied to the root and soil datasets aggregated to the Genus level to identify keystone taxa that drive the differences between treatments (Segata et al. 2011). LEfSe was run with a Wilcoxon and Kruskal-Wallis cut-off value of 0.01. An LDA cut-off value of 2 resulted in 182 and 145 marker genera for root and soil, respectively, and an LDA cut-off of 4 in 43 and 10 marker genera. In addition, analysis of compositions of microbiota with bias correction (ANCOMBC) was used to identify differentially abundant features (Lin & Peddada 2020). The Holm-Bonferroni method was applied to adjust p-values, and features with adjusted P values <0.01 were considered significant, resulting in 193 and 313 genus-level markers for root and soil datasets, respectively.