**2.4 Statistical data analysis**

Statistical analyses were executed using the Plymouth Routines in Multivariate Ecological Research (Primer v5.0) software package [42, 43]. For each microcosm treatment, all univariate indices were considered—abundance, species number (S), Shannon-Wiener diversity index (H'), Margalef's species richness (d), and Pielou's evenness (J'). Data were first tested to fulfill parametric requirements, Gaussian normality, and homogeneity of variances. Two tests were necessary: The Kolmogorov– Smirnov test to evaluate the first condition and the Bartlett to check the second condition. Once our data approximated normality, one-way ANOVA (1-ANOVA) was useful to determine the total significant difference between conditions. For multiple comparisons, the test of Tukey's HSD was applied (Statistica version 5.1).

For multivariate analyses: a non-metric Multi-Dimensional Scaling ordination (nMDS) by the Bray-Curtis matrix of similarity (square-root transformed abundance) were performed to detect a possible trend in the distribution of treatment, depending on the responses of nematode taxa or those of biological traits after ECs exposure. Hierarchical cluster analyses (CLUSTER) were used to confirm the results provided by the nMDS. The analysis of Similarity (ANOSIM) was considered to test the dissimilarities significance eventually noted between the nematode assemblages or each biological trait. Finally, The SIMPER analysis was useful in determining the contribution (cumulative contribution of 70%) of each species or functional group to the mean dissimilarity between treatments.

To detect the different relationships between nematode taxa and functional groups, PAST v3.26 software was used to perform correspondence analysis (CA) using two-dimensional plots. Also, the principal component analysis (PCA) (performed after transforming the data into log (x + 1)) is associated with the Pearson correlation (adopted via XLSTAT. 2019) to determine the targeted relationships.
