INTRODUCTION
Understanding the drivers of metacommunity structure across heterogeneous landscapes remains a fundamental question in ecology (Meynard et al., 2013; Viana & Chase, 2019). Under a niche selection perspective (Chase & Leibold, 2003; Hutchinson, 1959) metacommunity structure results from species sorting via environmental filtering, while under the neutral paradigm (Chase, 2005; Hubell, 2001) metacommunities are structured by stochastic dispersal and ecological drift. These processes are not mutually exclusive, as their relative importance in community assembly varies along a niche-dispersal continuum (Brown, Sokol, Skelton, & Tornwall, 2017; Cottenie, 2005; Gravel, Canham, Beaudet, & Messier, 2006). Scale dependency has been suggested to affect the perception of which process predominates, with environmental filtering prevailing at larger spatial scales and the role of stochastic processes increasing at finer scales (Viana & Chase, 2019). Yet, the apparent prevalence of environmental over spatial processes at regional scales (Chase, 2014) could be biased by the tendency of most metacommunity studies to quantify dispersal as a simple function of geographic distance (Biswas & Wagner, 2012), without considering landscape features such as topography or matrix heterogeneity. Accounting for matrix resistance when quantifying connectivity (McRae, Dickson, Keitt, & Shah, 2008) may enforce the role of dispersal as the primary driver of metacommunity structure (e.g., Resasco & Fletcher, 2021). This is particularly relevant when focusing on montane regions, where spatial connectivity among habitat patches can be strongly influenced by high topographic complexity and environmental heterogeneity due to steep elevational gradients (Graf, Kramer-Schadt, Fernández, & Grimm, 2007; Liu, Dudley, Xu, & Economo, 2018). At the same time, elevational gradients in topoclimatic parameters such as temperature or precipitation can impose strong environmental filtering on montane metacommunities (Hoiss, Krauss, Potts, Roberts, & Steffan-Dewenter, 2012; Leingärtner, Krauss, & Steffan-Dewenter, 2014) and are often considered to drive patterns of species richness (Peters et al., 2016; Rahbek, 1995) and community uniqueness (Wang et al., 2020) at regional scales. As montane communities are rapidly changing due to rising temperatures and anthropogenic disturbance (Rahbek et al., 2019; Steinbauer et al., 2018), it is crucial to gain a better understanding of how elevation and landscape mediate the interplay between environmental filtering and dispersal limitation as drivers of metacommunity structure across mountainous regions (e.g., Gálvez-Reyes et al., 2020).
Disentangling the drivers of metacommunity structure requires comprehensive empirical datasets from different functional groups, geographic regions and landforms, as the relative importance of spatial and environmental constraints is expected to be taxon- and context-dependent (He et al., 2020; Tonkin et al., 2018). It is thus difficult to extrapolate conclusions from the limited number of well-characterised montane metacommunities that have been studied to date (e.g., Benito, Fritz, Steinitz-Kannan, Vélez, & McGlue, 2018; Brodie & Newmark, 2019; Tonkin et al., 2017). More empirical studies from a diversity of montane biota across the globe are required to understand what general principles may be at play. Expanding metacommunity research to understudied geographic areas, and/or poorly known taxonomic groups, is complicated by the “taxonomic impediment” (Cicconardi, Fanciulli, & Emerson, 2013; Young, Proctor, deWaard, & Hebert, 2019), and has been often compromised by coarse taxonomic resolution (He et al., 2020; Verleyen et al., 2009). The development of metabarcoding provides new opportunities to accelerate studies on metacommunity structure across underexplored fractions of biodiversity and greatly increases their taxonomic resolution (Arribas, Andújar, Salces-Castellano, Emerson, & Vogler, 2021b; Bush et al., 2020; Martin et al., 2021; Zinger et al., 2019). Recent advances in field, laboratory and bioinformatic protocols for whole organism community DNA (wocDNA, Creedy et al., 2021) metabarcoding have led to improvements in both the efficiency for the generation of such high-resolution taxonomic inventories (Arribas, Andújar, Hopkins, Shepherd, & Vogler, 2016; Elbrecht, Vamos, Steinke, & Leese, 2018), and the reliability of α and β diversity estimates (Andújar et al., 2018; Creedy, Ng, & Vogler et al., 2019). New bioinformatic tools for the removal of noise generated by amplification and sequencing errors (e.g., Callahan et al., 2016; Edgar, 2016) and the filtering of spurious sequences resulting from co-amplification of nuclear mitochondrial pseudogenes (Andújar et al., 2021) allow to move beyond classical OTU (Operational Taxonomic Unit) clustering and define haplotype-level entities, or Amplicon Sequence Variants (ASVs; Callahan, McMurdie, & Holmes, 2017). In contrast to OTUs, ASVs have intrinsic biological significance and offer the possibility for direct comparisons among studies that use the same marker (Callahan et al., 2017; Porter & Hajibabaei, 2020), while they potentially improve community diversity estimates (Joos et al., 2020). The availability of reliable whole-community ASV and OTU datasets allows comparisons of spatial structure at haplotype and species levels, which can provide insights into the prevalence of stochastic vs. niche-based processes in community dynamics, as similar diversity patterns at both levels are predicted under neutral scenarios (Papadopoulou et al., 2011; Baselga, Gómez-Rodríguez, & Vogler, 2015). All the latest advances in the field of wocDNA metabarcoding are expected to lead to a better understanding of species diversity and community processes, particularly in historically intractable habitats, such as the soil (Arribas et al., 2021b).
Soil biodiversity is among the most complex and poorly known terrestrial biotas on Earth (Decaëns, 2010). The high structural complexity and heterogeneity of the edaphic environment are thought to facilitate species coexistence and drive patchy distributions of soil organisms at multiple spatial scales (Berg, 2012; Thakur et al., 2020), providing a particularly interesting template for metacommunity studies. Soil microarthropods, including the highly abundant and diverse groups of Acari, Collembola and Coleoptera, represent a major component of below-ground communities, with a broad range of functional roles in soil ecosystem services (Nielsen, 2019) and high levels of cryptic diversity (e.g., Cicconardi et al., 2013; Young et al., 2019). Recent metabarcoding studies (Arribas et al., 2021b; Zinger et al., 2019) have revealed an important role of stochastic processes and dispersal limitation in community assembly of soil microarthropods at within-habitat scale, in contrast to previous research that had emphasized environmental filtering in response to soil attributes as a major driver of community composition (e.g., Caruso, Schaefer, Monson, & Keith, 2019; Caruso, Taormina, & Migliorini, 2012; Gao, Liu, Lin, & Wu, 2016). Metabarcoding data also supports relatively low dispersal rates (Zinger et al., 2019) and high turnover of soil microarthropod assemblages even within short geographic distances (Arribas et al., 2021b), despite suggestions that long-distance passive dispersal might be prevalent in these small-bodied groups (Schuppenhauer, Lehmitz, & Xylander, 2019; Türke, Lange, & Eisenhauer, 2018). It remains to be assessed whether those differences among studies are due to the higher taxonomic resolution offered by metabarcoding, or due to context- or scale-dependent variation among systems (Berg, 2012; Ferrenberg, Martinez, & Faist, 2016). What both metabarcoding and morphology-based studies agree on is that habitat type can impose a strong filtering effect overriding other environmental and spatial processes, with forest vs. grassland habitats harbouring largely distinct metacommunities (Arribas et al., 2021b; Caruso et al., 2012; Rota et al., 2020). Yet equivalent comparisons among more structurally similar habitats (e.g., different forest types) remain limited.
Here we use both OTU- and ASV-level metabarcoding of entire communities to characterize soil microarthropod assemblages (Acari, Collembola and Coleoptera) across an isolated montane forest mosaic and evaluate the relative importance of forest type, topoclimatic variation and spatial/landscape factors in driving metacommunity structure. The Troodos mountain range harbours unique and understudied Mediterranean forest habitats within the island of Cyprus, one of Europe’s most vulnerable islands to climate change (Vogiatzakis, Mannion, & Sarris, 2016), and a major component of the Mediterranean biodiversity hotspot due to high levels of endemicity (Medail & Quezel, 1997). Troodos is characterised by complex topography and steep environmental gradients which, in combination with anthropogenic disturbance since ancient times (Delipetrou, Makhzoumi, Dimopoulos, & Georghiou, 2008), have created a mosaic consisting of five main forest habitat types that differ in area, altitudinal range and level of fragmentation (Figure 1). These comprise forests of: (i) the narrow endemic Cyprus Cedar - Cedrus brevifolia (Cb, hereafter), with a highly restricted distribution (~300 ha, between 900-1,400 m) in Western Troodos; (ii) the endemic Golden Oak - Quercus alnifolia (Qa), with a broad (~20,000 ha, 700-1,700 m) but highly fragmented distribution across Troodos; (iii) Black Pine -Pinus nigra pallasiana (Pn, ~3,500 ha, 1,450-1,950 m) and (iv) Stinking Juniper - Juniperus foetidissima (Jn, ~250 ha, 1,450-1,950 m), both narrowly distributed at the top of the highest peak (Chionistra, 1952 m) in Central Troodos; and finally (v) Calabrian Pine forest - Pinus brutia (Pb) the dominant habitat type, forming continuous and extensive forests across Troodos (~90,000 ha, 400-1,400 m). We generate metabarcode data for soil microarthropods across this habitat mosaic matrix to address the following questions: (i) Is forest habitat type the primary factor shaping metacommunity structure, in a similar way to that seen between grassland and forest? (Arribas et al., 2021b; Caruso et al., 2012), (ii) What is the relative contribution of spatial vs. environmental processes as drivers of within-habitat metacommunity structure? (Zinger et al., 2019), (iii) Focusing on the endemic Quercus alnifolia habitat, which is highly fragmented across Troodos, does habitat connectivity across the heterogeneous landscape play an important role in metacommunity structure? (Resasco & Fletcher, 2021). Finally, (iv) are α and β diversity patterns obtained using ASVs and OTUs equivalent and explained by similar spatial or topoclimatic factors? Apart from elucidating soil biodiversity dynamics in those poorly studied but highly vulnerable and precious Mediterranean forests, this system can provide insights into the utility of high-resolution community metabarcoding, in combination with fine-scale topoclimatic and landscape matrix information, for disentangling the drivers of metacommunity structure across mountainous regions and complex landscapes.