*2.6. Statistical Analyses*

R and Statistica 12.0 (StatSoft Inc., Poland, Krakow) software packages were used for all statistical analyses. The effects of the experimental factor (biostimulant/fertilizer), and the development phase term (based on the BBCH scale) on enzymatic activity in the soil were tested with three-way ANOVA (Tables S1–S5). Nitrogenase activity and agronomic parameters were tested using two-way ANOVA (Tables S6–S16). Homogeneous subsets of mean were identified by means of Duncan's test, at a significance level of α = 0.05.

As *year* was a variable, we used soil biochemical activity parameters (model 1) and agronomic parameters, together with biological nitrogen fixation (BNF) (model 2) ANOVA mixed models. The impact of two or three explanatory variables on the response variable, respectively, was assessed. In both cases, the following models were used:

$$y\_{i\bar{i}\bar{k}} = \mu + a\_{\bar{i}} + \beta\_{\bar{j}} + \gamma\_{\bar{k}} + (\beta\gamma)\_{\bar{j}\bar{k}} + (a\beta)\_{i\bar{j}} + (a\gamma)\_{\bar{k}} + (a\beta\gamma)\_{i\bar{j}\bar{k}} + \varepsilon\_{i\bar{j}\bar{k}} \pmod{1} \tag{1}$$

$$y\_{i\bar{k}l} = \mu + a\_{\bar{i}} + \gamma\_{\bar{k}} + (a\gamma)\_{\bar{i}k} + \varepsilon\_{i\bar{k}l} \pmod{2} \tag{2}$$

where: μ—is the overall average value, <sup>α</sup>*i*—is the effect of the operation of the random factor of *year* at level *i* (*i* = 1, 2, 3), β*j*—is the effect of the action of the fixed factor *term* at level *j* (= 1, 2 ... , 4), γ*k*—is the effect of the fixed fertilization factor at level *j* (= 1, 2 ... , 9 *j* (= 1, 2 ... , 9), with appropriate interactions of these factors, and *eijkl*—is the residual error.

In cases where the interaction of *year* with the other factors was significant, an analysis was carried out for each year separately. To estimate the cause-and-effect relationship between the studied soil biochemical activity parameters and agronomic parameters, principal component analysis (PCA) was used for each year separately, as well as for combined years. PCA was performed with the use of an appropriately scaled correlation matrix. PCA analysis was used to demonstrate the similarities between independent variables and determines the components that are a linear combination of the

variables considered. Accurate analysis of the principal components allows the identification of the initial variables that are the reference system for the remaining variables.

A heat map (using the heatmaply function in R), was proposed as a graphical presentation of appropriately transformed data of soil biochemical activity parameters, agronomic parameters, and biological nitrogen fixation (BNF). Data transformation using 'normalise' was used to compare and group different data.

Fertilization data were represented by colours. Cluster analysis allowed for the grouping of both soil biochemical activity and agronomic parameters after the application of biostimulators, and the effect of the fertilizers/biostimulators so that the degree of connection between the applied fertilization treatments *within* one group was the largest, while the degree of connection *between* groups was the smallest.

Grouping of tree diagrams was obtained by using the Ward Hierarchical Clustering method and the Euclidean distance measurement.
