Dissimilarity in community composition among forest sites
Dissimilarity in community composition among sampling sites was high and mainly determined by spatial turnover, with a very limited contribution of nestedness (ASVs: βSIM = 0.974, βSNE= 0.004; OTUs: βSIM = 0.961, βSNE = 0.006). This pattern was consistent when each taxonomic group was separately analyzed (all βSIM >0.952, all βSNE <0.015). The community dissimilarity matrices (βSOR or βSIM) of the three taxonomic groups were correlated among them at both ASV and OTU level (Mantel test, all r >0.240, all p-values <0.001). As the contribution of nestedness was minimal, we only report the results of the statistical analyses obtained using the community dissimilarity matrices calculated with the Simpson dissimilarity index (βSIM).
The NMDS analysis grouped the sampling points according to their respective forest habitat type, except those from Black Pine (Pn) and Stinking Juniper (Jn) forests which exhibited a greater overlap in the ordination (Figure 3). PERMANOVA analyses detected significant differences in community composition among forest habitat types, a factor explaining over 25-34% of the overall variation in community dissimilarity (Figure 3). These differences were significant for all habitat pairs (all p-values <0.020), as confirmed by pairwise comparisons (all p-values <0.019). The NMDS ordinations obtained from different community dissimilarity matrices (βSOR vs. βSIM and ASVs vs. OTUs) converged on highly similar solutions and were significantly concordant according to Procrustes tests (all r = 0.917, all p-values <0.001).
When dbRDA analyses were applied across all habitats, all three sets of explanatory variables (forest habitat type, spatial and topoclimatic variables) were retained as significant predictors of community dissimilarity (βSIM) (Table 2), although the largest fraction of the variation was clearly explained by habitat type (Figure 4). Pairwise comparisons testing for the effect of forest habitat type in community dissimilarity confirmed significant differences between all habitat pairs (all p-values <0.020). When dbRDA were applied at within-habitat scale, community dissimilarity (βSIM) of most forest types showed a significant relationship with topoclimatic predictors, except for the case of the Stinking Juniper (Jn) habitat where no predictor was significant (Table 2; Table S5). Once the Jn sampling sites were analysed together with the Pn sites according to the great overlap observed in the NMDS-based ordinations (Figure 3), a significant correlation between community dissimilarity (βSIM) and topoclimatic predictors was confirmed (Table 2). The predominant role of topoclimatic variation in explaining βSIM diversity patterns within each habitat was consistent across ASVs and OTUs levels (Table 2; Table S5).
Multivariate GLMs (mvGLMs) provided highly concordant inferences with those obtained using dbRDA (Table S6). Forest habitat type was the predictor that explained the largest fraction of the variation in community composition at across-habitat scale (R >0.121), with the spatial and topoclimatic predictors also included in the final model, although exhibiting much less explanatory power (R <0.044; Table S6). Significant differences in community composition between all habitat pairs were supported by pairwise comparisons (all p-values <0.013). Analyses at within-habitat scale showed that topoclimatic predictors significantly explained community composition within most habitat types. In those cases in which spatial predictors were also retained in the final model, their univariate contribution to the explained variance was usually lower than the one of the topoclimatic predictors (Table S6).
Matrix regressions with randomization (MRR) showed that community dissimilarity (βSIM) among Golden Oak (Qa) sampling sites at both ASV and OTU levels was positively correlated with the IBR matrix based on the spatial configuration of Qa patches (FRAIBR) and the topoclimatic distance matrix (ENVPC1-2). Both predictors were significantly retained in the best-fit models (Table 3), and there was no correlation between them (Mantel test, r <0.103, p-value >0.188). According to the sensitivity analyses, the FRAIBR scenarios that best explained the observed patterns of βSIM variation were those in which the non-Quercus matrix offered much higher resistance (10- or 100-fold) than the target habitat (Table S7). In concordance with the MRR analyses, distance decay models adjusted using GLMs showed a roughly linear decrease in community similarity (βSIM) with increasing FRAIBR or ENVPC1-2 distances, with both predictors being significant (Figure 5).