Canoco version 5

Patterns of spatiotemporal variability of physical and chemical variables were assessed via PCA, using the software CANOCO version 5 (Microcomputer Power, Ithaca, NY, USA). Variables were log transformed (except pH due to its logarithmic nature), centered, and standardized prior to the analysis. Comparisons among the three sampling periods for each physical and chemical variable were performed using repeated-measures ANOVA. Data were Box-Cox transformed when needed to fit homogeneity of variance assumption.

Environmental conditions were characterized by mean Landolt indicator values (Landolt 2010) of accompanying vascular plant species in herb and shrub layers calculated for circular plots centered at each target individual. Landolt indicator values describe ecological requirements of species in terms of climate (temperature, T; continentality, K; light, L) and soil parameters (moisture, F; reaction, R; nutrients, N; humus content, H; aeration, D; moisture variability, W) ranging from 1 (low) to 5 (high). Niche differences among the two parental species and their hybrids (identified with NewHybrids, see Results) were tested by comparing averaged indicator values of accompanying species in a multivariate analysis of variance (MANOVA) in R version 3.0.2. Visualization of particular ecological niches was performed using a principal component analysis (PCA) with Canoco version 5 (Microcomputer Power, Ithaca, US), based on standardized indicator values. Differences in additional environmental factors recorded at sites were inferred in R, either using an analysis of variance (ANOVA) followed by Tukey's HSD post hoc test (coverage of vegetation layers, species diversity) or, if the assumptions of ANOVA were not met, using a Kruskal–Wallis rank_ sum test followed by a multiple comparison test available in the package Pgirmess version 1.5.9 (Giraudoux 2014; bare rock cover).

Microsoft Excel Software was used to analyze the abundance and distribution of the key rhizosphere bacterial communities at the class and order level. Also, Shinyheatmap online tool (www1.heatmapper.ca/expression/) was used for the plotting of heatmaps using the relative abundance values. Equitability_J, Fisher_alpha, Berger-Parker and Chao-1 indices for diversity assessment were used for samples across rhizosphere and bulk soil samples, and the comparison of these indices was performed using Kruskal–Wallis test. The analyses were carried out using PAST version 3.20 (Hammer et al. 2001 ). Beta (β) diversity of the bacterial communities of the samples was determined using the Euclidean-based principal coordinate analysis (PCoA) and ANOSIM through 999 permutations. p < 0.05 was considered statistically significant for all datasets. Both PCoA and Principal Component Analysis (PCA) graphs were used to evaluate the relationship between bacterial communities and the measured physiochemical parameters, which were plotted using CANOCO version 5.0 (Microcomputer Power, Ithaca, NY). The effects of pH, total N and OM on the bacterial community distribution were determined using Canonical Correspondence analysis (CCA). Graphs and one-way ANOVA for functional measurements were achieved using GraphPad Prism 7 (GraphPad Software, California).