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.