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.