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 2 >0.121), with the spatial and
topoclimatic predictors also included in the final model, although
exhibiting much less explanatory power (R 2 <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).