Topoclimate as a driver of within-habitat structure
Analyses performed at within-habitat scale revealed the importance of both environmental and spatial variables on community composition in most cases (Table 2; Table S6), which is in line with previous research on soil microarthropod communities (Arribas et al., 2021b; Bahram et al., 2016; Ingimarsdóttir et al., 2012; Lindo & Winchester, 2009). Disentangling the relative contribution of environmental vs. spatial processes has proven to be challenging as environmental variation is often spatially autocorrelated, which can lead to a spurious inflation of the inferred environmental contribution (Clappe, Dray, & Peres-Neto, 2018; Vellend et al., 2014). In our study, spatial predictors generally explained less variance than topoclimatic factors in mvGLMs (Table S6), and their effect became non-significant after applying a forward selection approach in dbRDA (Table 2; Table S5). These results along with the relatively low degree of collinearity between spatial and topoclimatic axes (VIF <7; Vittinghoff, Glidden, Shiboski, & McCulloch, 2012), emphasize the role of environmental filtering as a key driver of metacommunity structure (Brown et al., 2017). Our results would complement several morphology-based studies suggesting that community composition of soil microarthropods is driven by environmental filtering, primarily in response to gradients of edaphic parameters (Caruso et al., 2019; Gao et al., 2016; Gan, Zak, & Hunter, 2019; Grear & Schmitz, 2005). However, they appear to contrast with the recent wocDNA metabarcoding study of Arribas et al. (2021b), where dispersal limitation was identified as the main driver of community assembly at within-habitat scale. This discrepancy cannot be attributed to taxonomic resolution, as both studies used very similar protocols to retrieve ASVs and OTUs, but it could be partly explained by differences in sampling scale, as the generally broader sampling extent of our study could enhance the role of environmental filtering as a consequence of encompassing higher environmental heterogeneity (Chase, 2014). Yet we also found environmental filtering to prevail in our narrowly distributed habitats (e.g., Cb, Pn) with observational scales slightly smaller (<7 km) than those of Arribas et al. (2021b). Additionally, the overall stronger effect of environmental filtering in our study system may reflect context-dependency (Soininen, 2014), with environmental processes playing a more important role in systems characterized by high topoclimatic heterogeneity. While in Arribas et al. (2021b) there were only moderate altitudinal gradients (~200-670 m elevation difference), our sampling spanned a steep elevational (1470 m elevation difference) and environmental gradient, with topoclimatic conditions varying greatly even across short distances, both within and across habitats (Figure S1). This phenomenon may be common in topographically complex regions, where dispersal limitation may actually be imposed by environmental heterogeneity rather than by geography per se (Liu et al., 2018), and points out the relevance of detailed topoclimatic characterization for understanding metacommunity structure within mountainous landscapes. However, it is noteworthy that the total variance explained by some models was relatively low (R 2ADJ<5-10%, Table 2). This was not unexpected, as it is a common finding among metacommunity studies (Cottenie, 2005), and has been traditionally attributed to other ecological processes that are not frequently measured (Vellend, 2010). Particularly, stochastic demographic processes including ecological drift in the absence of dispersal limitation (Bahram et al., 2016; Zinger et al., 2019) or priority effects via niche preempting (Fukami, 2015) have been hypothesized as relevant forces potentially interfering with community assembly in the soil environment. Additionally, we have not considered explicitly the effect of edaphic variables (e.g., organic matter, nutrient content or pH; Gao et al., 2016; Gan et al., 2019), although some of their variation is likely captured by forest habitat type and by certain topoclimatic variables, which are thought to influence specific soil attributes (horizon depth, moisture; Florinsky, 2012; Hillel, 2008).