Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests

Disentangling the relative role of environmental filtering and spatial processes in driving metacommunity


| INTRODUC TI ON
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 archetype (Chase, 2005;Hubbell, 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 nichedispersal continuum (Brown et al., 2017;Cottenie, 2005;Gravel et al., 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 spatial isolation as a simple function of geographical distance (Biswas & Wagner, 2012), without considering landscape features such as topography or matrix heterogeneity.Accounting for matrix resistance when quantifying connectivity (McRae et al., 2008) may enforce the role of spatial processes as the primary driver of metacommunity structure (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 et al., 2007;Liu et al., 2018).At the same time, elevational gradients in topoclimatic parameters such as temperature or precipitation can impose strong environmental filtering on montane metacommunities (Hoiss et al., 2012;Leingärtner et al., 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 spatial processes as drivers of metacommunity structure across mountainous regions (Gálvez-Reyes et al., 2021).
Disentangling the drivers of metacommunity structure requires comprehensive empirical data sets from different functional groups, geographical 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 wellcharacterized montane metacommunities that have been studied to date (Benito et al., 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 geographical areas, and/or poorly known taxonomic groups, is complicated by the "taxonomic impediment" (Cicconardi et al., 2013;Young et al., 2019).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, et al., 2021;Bush et al., 2020;Martin et al., 2021;Zinger et al., 2019).Recent advances in field, laboratory and bioinformatic protocols for community metabarcoding have led to improvements in both the efficiency for the generation of such high-resolution taxonomic inventories (Arribas et al., 2016;Elbrecht et al., 2018), and the reliability of α and β diversity estimates (Andújar et al., 2018;Creedy et al., 2019).New bioinformatic tools for the removal of noise generated by amplification and sequencing errors (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 us to move beyond classical OTU (operational taxonomic unit) clustering and define haplotype-level entities, or amplicon sequence variants (ASVs; Callahan et al., 2017).While OTUs are traditionally considered as proxies of species-level entities, ASVs represent haplotypes with intrinsic biological significance and offer the possibility for direct comparisons among studies that use the same marker (Callahan et al., 2017;Porter & Hajibabaei, 2020), potentially improving community diversity estimates (Joos et al., 2020).The availability of reliable whole-community ASV and OTU data sets 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 (Baselga et al., 2015;Papadopoulou et al., 2011).All the latest advances in the field of community 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, Andújar, Salces-Castellano, et al., 2021).
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 amplicon sequence variants, circuit theory, Cyprus, habitat connectivity, habitat filtering, operational taxonomic units services (Nielsen, 2019) and high levels of cryptic diversity (e.g., Cicconardi et al., 2013;Young et al., 2019).Recent metabarcoding studies (Arribas, Andújar, Salces-Castellano, et al., 2021;Zinger et al., 2019) have revealed an important role of stochastic processes and dispersal limitation in community assembly of soil microarthropods at the 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 (Caruso et al., 2012(Caruso et al., , 2019;;Gao et al., 2016).Metabarcoding data also support relatively low dispersal rates (Zinger et al., 2019) and high turnover of soil microarthropod assemblages even within short geographical distances (Arribas, Andújar, Salces-Castellano, et al., 2021), despite suggestions that long-distance passive dispersal might be prevalent in these small-bodied groups (Schuppenhauer et al., 2019;Türke et al., 2018).It remains to be assessed whether those differences among studies are due to the higher taxonomic resolution offered by metabarcoding, or are due to context-or scale-dependent variation among systems (Berg, 2012;Ferrenberg et al., 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, Andújar, Salces-Castellano, et al., 2021;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 whole-organism community DNA (wocDNA, Creedy et al., 2021) metabarcoding at both OTU and ASV levels 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 (hereafter Troodos) harbours unique and understudied Mediterranean forest habitats within the island of Cyprus, one of Europe's islands most vulnerable to climate change (Vogiatzakis et al., 2016), and a major component of the Mediterranean biodiversity hotspot due to high levels of endemicity (Medail & Quezel, 1997).Troodos is characterized by complex topography and steep environmental gradients which, in combination with anthropogenic disturbance since ancient times (Delipetrou et al., 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), with a highly restricted distribution (~300 ha, between 900 and 1400 m) in Western Troodos; (ii) the endemic golden oak-Quercus alnifolia (Qa), with a broad (~20,000 ha, 700-1700 m) but highly fragmented distribution across Troodos; (iii) black pine-Pinus nigra pallasiana (Pn, ~3500 ha, 1450(Pn, ~3500 ha, -1950 m) m) and (iv) stinking juniper-Juniperus foetidissima (Jn, ~250 ha, 1450(Jn, ~250 ha, -1950 m) 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-1400 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, Andújar, Salces-Castellano, et al., 2021;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 Q. 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.

| Soil sampling and sample processing
During 2019, from mid-April to mid-June, we collected soil samples from 44 sites representing the main five forest habitat types of the Troodos mountain range (described above, Figure 1; Table S1).Our sampling scheme covered the full extent of the distribution and altitudinal range that the five tree species exhibit on Troodos, spanning over 65 km along an east-west axis and 1500 m of elevation range (Figure 1; Table S1).We collected two soil samples per sampling site corresponding to the superficial (1 m 2 of leaf litter and humus, 5 cm depth) and the deep layer (30 cm diameter, 30 cm depth, comprising ~20 L of soil) as described in Arribas, Andújar, Salces-Castellano, et al. (2021).The 88 soil samples were subsequently processed using a standardized flotation-Berleseflotation protocol to extract the soil mesofauna as detailed in Arribas et al. (2016), Arribas, Andújar, Salces-Castellano, et al. (2021).Once retrieved, each mesofauna sample was submitted to a final cleaning procedure, which involved flotation and posterior filtering using two sieves of 1-mm and 41μm mesh.This protocol allows discarding unwanted soil components (sediments and vegetal remains) while retrieving two subsamples of bulk arthropod specimens, retained in each sieve according to their body size (typically Acari and Collembola vs. Coleoptera), which are suited for "clean" extraction of wocDNA (Arribas, Andújar, Salces-Castellano, et al., 2021).During bulk-sample processing, we additionally selected "voucher" specimens of Acari, Collembola and Coleoptera representing the whole range of morphological variability of these groups observed across samples and habitats.A total of 176 bulk subsamples and 344 "voucher" specimens were preserved at −20°C in 100% ethanol for molecular analyses.

| DNA extraction, PCR amplification and sequencing
We extracted DNA from each bulk subsample using the Biosprint 96 DNA Blood Kit (Qiagen) on a Thermo KingFisher Flex automated extraction instrument.We quantified the DNA concentration using a NanoDrop spectrophotometer and combined the extracts of each pair of subsamples at a 1:10 ratio (Arribas, Andújar, Salces-Castellano, et al., 2021).Then, we PCR (polymerase chain reaction) amplified the 3′ end of the cytochrome c oxidase subunit I (COI) barcode region corresponding to the 418-bp bc3′ fragment using primers Ill_B_F (Shokralla et al., 2015) and Fol-degen-rev (Yu et al., 2012).Amplifications were performed following the PCR conditions described in Arribas et al. (2016).We carried out independent PCR replicates using as template two different DNA dilutions (1:10 and 1:100).PCR products were visualized on agarose gel and three amplicons (including two PCR replicates of the 1:10 DNA template dilution and one PCR replicate of the 1:100 DNA template dilution) were pooled per sample, which were purified using a magnetic beadbased protocol (Agencourt AMPure XP).We included three negative controls corresponding to different wet-laboratory steps (lysis, DNA extraction and PCR amplification).The 88 metabarcoding samples and the three negative controls were used for a dual-indexed library preparation following the Nextera XT DNA workflow (llumina) and were sequenced on a paired-end 2 × 300-bp lane of an Illumina MiSeq platform at the Earlham Institute (Norwich, UK).
We individually extracted DNA from each of the "voucher" specimens as described above.We amplified Folmer's COI barcode region (658-bp, overlapping with the 418-bp metabarcoding fragment), using the primers Fol-degen-for and Fol-degen-rev (Yu et al., .Green hatched areas represent zones where P. brutia (light green) forests are present according to our own surveys but whose extent has not been mapped in the publicly available cartography, probably because they are the result of plantations, have been affected by fires and/or present very low tree densities 2012) and following the PCR conditions described in Arribas et al. (2016).The PCR products were purified and bidirectionally Sanger sequenced (Macrogen).These "voucher" sequences were used for constructing a local reference database, which allowed us (i) to improve the higher-rank taxonomic assignment of the metabarcode reads and thus maximize the number of ASVs retained for each of the three taxonomic groups and (ii) to apply a recently developed method for read filtering (metamate, Andújar et al., 2021, see below).

| Illumina read processing and filtering
Detailed information on the read processing and filtering pipeline is summarized in Table S2.Briefly, we demultiplexed raw reads allowing no mismatch in the dual-index pair.Then, we used fastqc version 0.11.7 (Andrews, 2010) to quality check raw reads and cutadapt version 2.10 (Martin, 2011) to trim primers and filter out raw reads exhibiting any variation from expected primer length and composition.Subsequently, we used pear version 0.9.11 (Zhang et al., 2014) to merge forward and reverse reads.Each metabarcoding sample was then separately quality filtered, dereplicated discarding singletons, length filtered retaining only reads 416-420 bp, denoised using unoise3 and de novo chimera filtered using uchime3 as implemented in vsearch version 2.9.1 (Rognes et al., 2016).Once denoising was performed, reads from all metabarcoding samples were pooled and again dereplicated (discarding no sequences) to generate a catalogue of unique putative haplotypes (ASVs).Subsequently, we ran blast to compare all ASVs against a combined database composed of the NCBI nt collection (accessed November 2020) and a curated reference catalogue including the 344 Sanger sequences of the "voucher" specimens plus 561 previously available sequences corresponding to soil lineages of Acari, Collembola and Coleoptera (Arribas et al., 2016;Arribas, Andújar, Salces-Castellano, et al., 2021).Based on the blast output we assigned the ASVs to high-rank taxonomic levels, by applying the weighted lowest common ancestor algorithm in megan6 (Huson et al., 2016; see also Hleap et al., 2021).Only ASVs assigned to Acari, Collembola or Coleoptera were retained and used for downstream analyses.We further filtered the ASVs using metamate version 0.1b18 (Andújar et al., 2021), a novel approach aiming at removing putative nuclear copies of mitochondrial DNA (NUMTs; Lopez et al., 1994) and other types of low-frequency erroneous sequences from denoised metabarcoding data sets.This software allows the application of multiple read-abundance filtering strategies and posterior evaluation of their effects on the prevalence of known authentic mitochondrial haplotypes and presumed nonmitochondrial copies (e.g., those violating the reading frame or expected length, as expected for NUMTs and erroneous sequences) in the final filtered data set (Andújar et al., 2021).The designation of ASVs as either known mitochondrial haplotypes or presumed nonauthentic sequences was performed using both global (GenBank, NCBI nt collection) and local ("voucher" specimens) reference databases.We selected the most stringent filtering solution to ensure the removal of most erroneous sequences (see Appendix S1 for details on ASV designation and filtering in metamate).Subsequently, we used vsearch to generate a read-count community table of the metamate-filtered ASVs by matching them with a 100% identity value against the raw read data set before dereplicating, length filtering and denoising.We further filtered these community tables by removing ASVs showing abundances of two or fewer reads and also those whose contribution to the total number of reads per taxonomic group and library was <1%.Finally, filtered read-count community tables were converted to presence/absence tables (see Jurburg et al., 2021).Negative controls were processed alongside actual samples throughout the filtering workflow.
Using the ASV and OTU data sets, we calculated α diversity (richness) for each sampling site and habitat type.We also calculated community dissimilarity among sampling sites and between soil layers using the Sørensen index (β SOR ), which we decomposed into its additive components (Simpson dissimilarity index or spatial turnover without the effect of variation in richness, β SIM , and nestedness, Baselga, 2010).Then, we used these community dissimilarity matrices to calculate the local contribution to β diversity (LCBD), a comparative indicator of the uniqueness of each sample in terms of community composition (Legendre & De Caceres, 2013).We calculated the above diversity metrics using the betapart (Baselga & Orme, 2012) and adespatial (Dray et al., 2021) r packages.

| Characterization of sampling sites using spatial, environmental and landscape variables
A number of spatial, environmental and landscape variables were calculated to characterize the sampling sites and quantify the distances among localities, taking into account the high topographic complexity and environmental/landscape heterogeneity of the study area.We calculated pairwise weighted topographic distances (SPA TWD ) based on a digital elevation model (DEM) at 90-m resolution using the topoDistance r package (Wang, 2020) as described in the Appendix S1.We also generated a set of high-resolution environmental variables (at 90-m resolution) for Cyprus, by spatial interpolation of temperature and precipitation layers at lower resolution These variables are known to affect the water-energy dynamics and to explain patterns of diversity in several arthropod groups, including Coleoptera (Hawkins et al., 2003).We extracted values of each interpolated variable along with the elevation for all sampling sites and for 500 randomly distributed points throughout Cyprus, to avoid potential biases resulting from only considering conditions at focal sites.We then applied a principal component analysis (PCA; Figure S1) to reduce the dimensionality of the data set and eliminate covariance among variables and we retained the two first principal components (PCs).PC1 (accounting for 81.3% of the variation) was positively correlated with altitude and precipitation variables and negatively with temperature variables, and PC2 (8.0%) was positively correlated with the topographic wetness index (Table S3).The PC1 and PC2 scores for each sampling site were considered as topoclimatic predictors (ENV PC1 and ENV PC2 ) for downstream analyses.We also calculated a Euclidean distance matrix among sampling points based on the obtained scores of the two retained PCs, which was used as a topoclimatic predictor (ENV PC1-2 ) for matrix regression analyses.
Finally, we applied a circuit theory approach (McRae, 2006) to quantify habitat connectivity among the Quercus alnifolia sampling sites.We focused specifically on the Q. alnifolia habitat because it is broadly distributed but highly fragmented across Troodos, in contrast to the Cedrus brevifolia, Pinus nigra and Juniperus foetidissima forests which are very narrowly distributed and the Pinus brutia forest which is very extensive with largely continuous distribution (Figure 1).Based on the assumption that dispersal among isolated forest fragments can be impeded by the presumed lower habitat suitability of the surrounding landscape (Brodie & Newmark, 2019), we built an isolation-by-resistance (IBR) scenario of connectivity (FRA IBR ) defined by the distribution of Q. alnifolia forest patches according to the existing cartography (see Appendix S1).Two alternative IBR scenarios were constructed for comparison: the TRI IBR representing the topographic complexity of the study area as estimated by the terrain roughness index (Title & Bemmels, 2018) and the NULL IBR representing a completely "flat landscape" with a fixed resistance (=1) value assigned to all cells.Resistance distances among all Q. alnifolia sampling points (n = 11) were calculated under each alternative IBR scenario (FRA IBR , TRI IBR and NULL IBR ) in circuitscape version 4.0.5 (McRae & Beier, 2007).

| Statistical analyses
All statistical analyses were performed at sampling site level, after integrating the community tables of the two soil layer samples (leaf litter and deep soil) for all three taxonomic groups (Acari, Collembola and Coleoptera) into a single matrix.Each of the analyses was performed at both ASV and OTU levels.

| Richness and community uniqueness
We used ANOVAs to test for significant differences in richness (RICH) and community uniqueness (LCBD, local contribution to β diversity) between forest habitats.We used generalized linear mixed models (GLMMs) to analyse the relationship between RICH or LCBD per site and the topoclimatic variables (ENV PC1 and ENV PC2 ) as predictors, with latitude and longitude as covariates.We built GLMMs fitting forest habitat type as a random effect in order to account for nonindependence among samples from the same forest habitat (see Appendix S1).

| β diversity
The clustering of sampling sites according to their community composition was visualized using nonmetric multidimensional scaling (NMDS) assuming two dimensions (k = 2).Differences among forest habitat types were tested for significance using permutational multivariate analysis of variance (PERMANOVA).We used symmetric Procrustes analyses to statistically assess nonrandomness among NMDS ordinations calculated from different community dissimilarity matrices (β SOR vs. β SIM and ASVs vs. OTUs; see Peres-Neto & Jackson, 2001).These analyses were performed using the vegan (Oksanen et al., 2020) and pairwiseAdonis (Martinez Arbizu, 2020) r packages.
The effect of forest habitat type and of the spatial and topoclimatic predictors on β diversity patterns was tested using distancebased redundancy analysis (dbRDA; Legendre & Anderson, 1999), ).The best-fit model was selected after ensuring there were no issues of multicollinearity (variance inflation factors, VIF <10; e.g., Tonkin et al., 2016).Finally, the best-fit models for the across-habitat analyses were used to partition the variance explained exclusively by each variable group (forest habitat type, space and topoclimate) and their intersections using the adjusted coefficient of determination (R 2 ADJ ) (Zinger et al., 2019).These analyses were performed using the BiodiversityR (Kindt & Coe, 2005) and vegan r packages.
To validate the dbRDA inferences (Jupke & Schäfer, 2020), we also applied multivariate generalized linear models (mvGLMs) as implemented in the mvabund r package (Wang et al., 2012).Unlike dbRDA, which is a distance-based approach, this model-based method fits a separate GLM per species and performs resamplingbased hypothesis testing for community-level effects of predictors.Here, incidence (presence/absence) data sets were used as input and tested against the same sets of predictors as in dbRDA.
Models were fitted using a binomial error distribution and a log link function.We first assessed the predictor significance using single-term models and only those predictors showing a significant effect were used to assemble a full model.The best-fit model was built following a backward stepwise selection approach (How et al., 2020).Term significance was assessed using likelihood ratio tests, PIT-trap resampling and 999 bootstrapping iterations (Wang et al., 2012).
Finally, we assessed the effect of habitat connectivity on the community composition of the Quercus alnifolia (Qa) sampling sites by applying matrix regressions with randomization (MRR; Wang, 2013).Specifically, the community dissimilarity matrices based on the Simpson dissimilarity index (β SIM ) were used as response variables, while the explanatory variables were selected among: (i) resistance due to habitat fragmentation (FRA IBR ), (ii) resistance due to topographic complexity (TRI IBR ), (iii) resistance due to a "flat landscape" (NULL IBR ), (iv) weighted topographic (SPA TWD ) distances and (v) topoclimatic (ENV PC1-2 ) distances.We assembled a full model including all explanatory matrices and built the best-fit model following a backward stepwise selection approach, using 999 permutations for the significance tests (e.g., Ortego et al., 2015).The unique contribution of each predictor to the total variance explained by the best-fit model was quantified using hierarchical variance partitioning analysis in the hier.partr package (Walsh & Mac Nally, 2003).Once the best-fit model was selected, we visualized the relationship between community similarity (1 − β SIM ) and the explanatory distance/ resistance matrices and fitted a distance-decay of community similarity curve (Gómez-Rodríguez & Baselga, 2018;Nekola & White, 1999).Specifically, we fitted a negative exponential function to univariate GLMs assuming a Gaussian error distribution and a log link function as implemented in the betapart r package.

| Richness and uniqueness of soil microarthropod communities
After denoising and chimera filtering, the removal of putative spurious sequences by metamate, followed by additional filtering of com-

| Effect of topoclimatic and spatial variables on patterns of richness and uniqueness
Regression analyses (GLMMs/GLMs; see Appendix S1) showed a significantly negative relationship between average richness of ASVs per site and longitude (Lon), indicating that community richness decreased towards the east of the Troodos mountain range (Table 1; Table S4).To ensure that this relationship was not biased by the Cedrus brevifolia (Cb) sites, which are geographically restricted to the westernmost part of the study area and have the highest ASV richness (Figures 1 and 2), additional analyses excluded these sites.
These analyses consistently supported the significant effect of longitude on ASVs richness across forest habitats (95% confidence interval [CI]: [−34.968] to [−4.032]).Conversely, the richness of OTUs per site was only explained by the topoclimatic predictor ENV PC2 (Table 1; Table S4).Similarly, we found that LCBD estimates at both ASV and OTU levels were significantly correlated with the topoclimatic variables, as summarized with the ENV PC1 and ENV PC2 predictors (Table 1; Table S4).

| 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 analysed (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 > .240,all p <.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 Pinus nigra (Pn) and Juniperus foetidissima (Jn) forests which exhibited a greater overlap in the ordination (Figure 3).PERMANOVAs 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 < .020),as confirmed by pairwise comparisons (all p < .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 = .917,all p < .001).

| Effect of topoclimatic and spatial variables on metacommunity structure
When dbRDAs 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 < .020).When dbRDAs 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 Juniperus foetidissima (Jn) habitat where no predictor was significant (Table 2; Table S5).Once the J. foetidissima sampling sites were analysed together with Pinus nigra 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 ASV and OTU levels (Table 2; Table S5).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 > .121),with the spatial and topoclimatic predictors also included in the final model, although exhibiting much less explanatory power (R 2 < .044;Table S6).Significant differences in community composition between all habitat pairs were supported by pairwise comparisons (all p < .013).
Analyses at within-habitat scale showed that topoclimatic predictors F I G U R E 2 Average richness (top panels) and community uniqueness (local contribution to β diversity-LCBD, bottom panels) across sampling sites, as estimated at ASV (left panels) and OTU (right panels) levels, per forest habitat type.Colours and habitat codes are as in Figure 1.Inset graphs show the contribution of each taxonomic group (Acari, dark grey; Collembola, grey; Coleoptera, light grey) to the cumulative richness (γ diversity) per forest habitat type.Shared letters below the box-plots indicate that differences between the respective habitats are not statistically significant (p > .05)after post hoc Tukey's tests 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 that of the topoclimatic predictors (Table S6).

| Effect of habitat connectivity on community composition of the Quercus alnifolia habitat
MRR showed that community dissimilarity (β SIM ) among Quercus alnifolia (Qa) sampling sites at both ASV and OTU levels was positively TA B L E 1 Results of model selection and averaging testing for the relationship between average richness (RICH) or local contribution to β diversity (LCBD) at ASV and OTU levels as response variables and topoclimatic variation (ENV PC1 and ENV PC2 ; see Table S3 Notes: Latitude (Lat) and longitude (Lon) were included as covariates.Predictors excluding the value 0 in their 95% confidence intervals (CI) are indicated in bold type and their effects were considered significant.For each final model including only predictors considered significant (in bold type), marginal (R 2 m , variance explained by fixed effects) and conditional (R 2 c , variance explained by both fixed and random effects) coefficients of determination are reported.
F I G U R E 3 Nonmetric multidimensional scaling (NMDS) ordination of sampling sites according to community dissimilarity (Simpson dissimilarity index, β SIM ), at ASV (left panel) and OTU (right panel) levels.Circles correspond to sampling sites, with circle size representing sample richness.The percentage of explained variation (R 2 ) and the significance of forest habitat type as a grouping factor based on PERMANOVA are reported on the top of each plot.Colours and habitat codes are as in Figure 1 correlated with the IBR matrix based on the spatial configuration of Q. alnifolia patches (FRA IBR ) and the topoclimatic distance matrix (ENV PC1-2 ).Both predictors were significantly retained in the bestfit models (Table 3), and there was no correlation between them (Mantel test, r < .103,p > .188).According to the sensitivity analyses, the FRA IBR 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 FRA IBR or ENV PC1-2 distances, with both predictors being significant (Figure 5).

TA B L E 2
Best-fit models from distance-based redundancy analyses (dbRDA) on community composition (Simpson dissimilarity index, β SIM ) of sampling sites at ASV and OTU levels, performed either with all sites (across habitats) or separately for each of the forest habitat types

ASVs (haplotypes)
OTUs (3% lineages) Note: The explanatory variables were forward selected from full models containing three sets of predictors: forest habitat type (HAB), spatial (SPA PCNMi ) and topoclimatic (ENV PCi ) variables.The adjusted coefficient of determination (R 2 ADJ ) for each model is provided.The Pinus nigra (Pn) and Juniperus foetidissima (Jn) sampling sites were both separately and jointly analysed according to results of NMDS-based ordinations (Figure 2).Predictor information for models with no significant variables (null) is replaced by dashes.

F I G U R E 4
Venn diagram illustrating the partitioning of explained variance in community composition among the three sets of explanatory variables (forest habitat type, spatial and topoclimatic predictors) and their intersections.Variance partitioning was conducted on the best-fit dbRDA models, which only included significant variables after a forward selection procedure (Table 2)  Notes: The explanatory variables were backward selected from full models containing the following distance matrices: "flat" scenario (NULL IBR ), topography (SPA TWD ), topoclimate (ENV PC1-2 ), topographic complexity (TRI IBR ) and forest fragmentation (FRA IBR ).The significance of both explanatory terms retained in the best-fit model and those rejected during the backward selection procedure are reported.The independent (R I 2 ) and joint coefficient of determination (R J 2 ) of each predictor retained in the best-fit final model are provided.For display purposes, the decay curves were fitted using generalized linear models (GLMs) with a negative exponential function.Regression slopes, coefficients of determination (pseudo-R 2 ) and p-values are provided on the top of each panel.See Table 3 for multivariate matrix regression analyses assessing the independent contribution of each predictor to community dissimilarity diversity eastwards due to anthropogenic disturbance.Therefore, species richness presumably responded to niche-based processes, while haplotype richness responded to historical contingencies.
Our results demonstrate the utility of combining OTUs with ASVs obtained by stringent filtering for characterizing diversity patterns of complex assemblages across heterogeneous landscapes.

| Forest type shapes regional metacommunity structure
The community composition of soil microarthropods across Troodos was largely explained by forest type, irrespective of spatial and topoclimatic variation (Table 2; Table S6; Figure 4), suggesting a primary role of habitat filtering as a driver of regional metacommunity structure.The importance of habitat type in shaping the composition of soil microarthropod assemblages has been demonstrated by previous research comparing forest vs. grassland habitats (Arribas, Andújar, Salces-Castellano, et al., 2021;Caruso et al., 2012) or different grass and shrub species (Coulson et al., 2003; see also Doblas-Miranda et al., 2009), but to our knowledge, this is the first study to compare different forest types in a systematic way.Given the mosaic nature of Troodos (Figure 1), our sampling scheme revealed high turnover among nearby sampling sites of different woodland habitats (Figure 3), highlighting a major role of underlying forestspecific edaphic features for the community assembly of soil microarthropods (Eissfeller et al., 2013).Leaf litter and vegetation type are documented to influence the soil environment through changes in microhabitat availability and physicochemical edaphic properties (Berg & McClaugherty, 2008), and such local abiotic features might be driving a scenario of species sorting (Leibold et al., 2004) or even the existence of largely separate metacommunities inhabiting each forest type.In contrast to the general pattern, the communities of the two highland habitats of Pinus nigra (Pn) and Juniperus foetidissima (Jn) are more similar to each other (Figure 3; although still significantly different based on pairwise comparisons), suggesting that harsh climatic conditions at those high altitudes might be imposing a stronger environmental filter than forest-associated soil attributes.Alternatively, the spatial configuration of this highland woodland (composed of isolated small stands of J. foetidissima embedded in a much larger matrix of P. nigra) might facilitate dispersal among habitat patches and partly counteract the effects of habitat filtering (Leibold et al., 2004;Logue et al., 2011).These two hypotheses are not mutually exclusive and may jointly contribute to the higher similarity in community composition between these highland assemblages.

| Topoclimate as a driver of within-habitat metacommunity structure
Analyses performed at the 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, Andújar, Salces-Castellano, et al., 2021;Bahram et al., 2016;Ingimarsdóttir et al., 2012;Lindo & Winchester, 2009).Disentangling the relative contribution of environmental vs. spatial factors 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 et al., 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 nonsignificant 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 et al., 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;Gan et al., 2019;Gao et al., 2016;Grear & Schmitz, 2005).However, they appear to contrast with the recent wocDNA metabarcoding study of Arribas, Andújar, Salces-Castellano, et al. (2021), where dispersal limitation was identified as the main driver of community assembly at the 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., Cedrus brevifolia, Pinus nigra) with observational scales slightly smaller (<7 km) than those of Arribas, Andújar, Salces-Castellano, et al. (2021).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, Andújar, Salces-Castellano, et al. (2021) 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 2 ADJ < 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).In particular, 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 pre-empting (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; Gan et al., 2019;Gao et al., 2016), although some of their variation is probably 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).

| Role of habitat connectivity in the assembly of the Quercus alnifolia metacommunity
Habitat fragmentation has been shown to alter the relative importance of spatial vs. environmental processes as drivers of metacommunity structure (Jamoneau et al., 2012), but the way we measure spatial distances among habitat patches and account for the effects of the surrounding matrix may affect our interpretation of the predominant processes (Resasco & Fletcher, 2021;Watling et al., 2011).In the case of the highly fragmented Quercus alnifolia (Qa) habitat, circuit theory-based connectivity modelling demonstrated that isolation estimates accounting for fragmentation and matrix resistance consistently performed better at explaining turnover than those based on topography, with the latter performing better than null models assuming a homogeneous matrix (Table 3; Table S7).This finding aligns with recent studies that have documented the ecological importance of dispersal corridors for community assembly across different taxonomic groups (Firmiano et al., 2021;Marrec et al., 2021;reviewed in Fletcher et al., 2016).However, this approach relies on the assumption that all species of the metacommunity respond similarly to landscape heterogeneity.Future studies integrating metabarcoding with morphological information derived from local "voucher" reference collections could facilitate the implementation of species-specific analyses accounting for ecological and trait variation (Brodie & Newmark, 2019;Hartfelder et al., 2020).Despite these limitations, our results provide empirical evidence of an important effect of habitat fragmentation on soil microarthropod metacommunity structure across the Q. alnifolia (Qa) forest patches, with the role of environmental filtering remaining equally significant (Table 3).Interestingly, we also obtained equivalent distancedecay of community similarity curves based on connectivity or on topoclimatic distances (Figure 5), although the topoclimatic variables were not spatially autocorrelated (Mantel test, r < .103,p > .188)and their shared variance with connectivity-based predictors was very low (<2%; Table 3).The similar slopes of the decay curves for ASVs and OTUs based on spatial distances are compatible with a neutral dispersal-constrained model (Baselga et al., 2015), while the equivalent pattern based on topoclimatic distances could be generated under certain scenarios of high dispersal and narrow ecological niches (Baselga et al., 2013).Taken together, our results suggest that spatial and niche-based processes have jointly shaped turnover patterns across the Q. alnifolia (Qa) habitat, although each process may affect different fractions of the metacommunity (e.g., some species and/or intraspecific entities may have low dispersal propensity and wide topoclimatic niches, while others might be good dispersers with narrow niches).Future species-specific analyses may help to refine these conclusions.

| Utility of ASVs for community ecology
All the above conclusions regarding the predominant processes shaping community composition of soil microarthropods across the Troodos forests were very similar when based on OTUs or ASVs, contributing to the broader discussion about whether OTUs should be replaced by ASVs in metabarcoding studies (Callahan et al., 2017;Porter & Hajibabaei, 2020).As recent read-filtering methods have overcome the need of clustering to account for amplification and sequencing errors, the user-defined OTUs (traditionally used as proxies of species-level entities) could be redundant.However, the exclusive use of ASVs could affect the biological conclusions drawn from biodiversity analyses, as patterns of haplotypic diversity can reflect demographic attributes of populations, and do not always coincide with diversity patterns at the species level (Martin et al., 2021).In our system, while β diversity estimates at ASV and OTU levels were correlated and probably shaped by the same ecological processes as described above, we observed distinct richness patterns between them, which were explained by statistically significant differences in geographical or topoclimatic predictors.OTU richness per site was primarily explained by topoclimatic conditions, with assemblages hosting fewer OTUs as elevation and precipitation increased and temperature decreased (Table 1; Figure S2), thus following the general rule of declining species richness with increasing elevation (Rahbek, 1995), commonly interpreted as an outcome of environmental filtering driven by temperature or productivity gradients (Graham et al., 2014;Peters et al., 2016).In contrast to OTUs, ASV richness varied significantly along a longitudinal axis with local communities harbouring more haplotypes westwards (Table 1; Figure S2), which might be interpreted as a signature of higher on average intraspecific genetic diversity in the western part of the mountain range, which has historically been less affected by anthropogenic disturbance (Delipetrou et al., 2008).This finding is in accordance with population genetic studies of forest trees that observed high genetic diversity in the western populations of Pinus brutia (Eliades et al., 2018) and Cedrus brevifolia (Eliades et al., 2011).This is potentially a consequence of the local topography facilitating the maintenance of higher effective population sizes in this region during the Pleistocene climatic oscillations, and/or of less intensive historical human impact (livestock grazing and logging) than in Eastern Troodos (Eliades et al., 2018).The decline of haplotypic diversity in soil microarthropod assemblages eastwards may therefore indicate incipient biodiversity loss, as genetic variation tends to be eroded more quickly than species diversity under scenarios of global change (Balint et al., 2011).However, such differences between species and haplotype diversity were not reflected in patterns of community uniqueness, as our ASV-and OTU-based estimates of LCBD (local contribution to β diversity) were correlated and similarly explained by topoclimatic variation (Table 1; Figure S2), without any clear signature of historical contingencies, as those probably affecting longitudinal haplotypic richness patterns across Troodos.Our results therefore highlight the complementarity of OTUs and ASVs for community metabarcoding, as such side-by-side comparisons can help to detect processes that produce uncoupled patterns between the two levels of diversity (Reisch & Schmid, 2019).
Our ASV-level analyses were facilitated by the application of the metamate tool that utilized local and public reference sequence databases to discard nonauthentic ASVs and retain only true biological sequence variants (Andújar et al., 2021).Although applying the most stringent filtering in metamate might have caused the removal of valuable rare biological haplotypes, appreciable intraspecific genetic variation (on average 2.35 ASVs per OTU) was still retrieved and produced reasonable haplotype diversity patterns as explained above.
Based on our results, we advocate stringent ASV filtering, as it can provide informative data sets without compromising the required reliability for haplotype-level metabarcoding.Even if we cannot be fully confident that all erroneous haplotypes were filtered out, as the performance of the approach depends on the completeness of the reference sequence catalogue (Andújar et al., 2021), the future incorporation of intraspecific genetic data in local reference databases will provide further confidence and the opportunities for intra-OTU analyses of ASV variation (Elbrecht et al., 2018;Zizka et al., 2020).

| CON CLUS IONS
This study highlights the power of combining ASV-and OTU-level community metabarcoding, enabled by the application of stringent filtering strategies, for the description of spatial biodiversity patterns of complex communities in understudied regions (Cyprus) and environments (soil), overcoming previous limitations of the taxonomic impediment, low-resolution data and noise due to the presence of spurious sequences.The wide implementation of harmonized field, laboratory and bioinformatic protocols for community metabarcoding of unexplored assemblages will increase the comparability of data sets from across the globe (Arribas, Andújar, Bidartondo, et al., 2021), providing the basis for broad-scale analyses of metacommunity patterns that would enable drawing more general conclusions on the consistency or context-dependency of ecological processes across spatial scales and fractions of biodiversity.Additionally, the ease with which all species in local communities can be characterized at the population genetic level using metabarcoding with stringent filtering raises the prospect for modelling demographic processes for each of the component species (Overcast et al., 2019(Overcast et al., , 2021)).Such an approach has the potential to elucidate historical and contemporaneous community responses to environmental heterogeneity and dispersal limitation at a much finer resolution than the summary statistics currently applied in whole-community metabarcoding.

F
Geographical location of sampling sites throughout the Troodos mountain range in Cyprus (top panel) and distribution of the main five forest habitat types (bottom panel).Sampling sites and forest habitat distribution are coloured as follows: Pb, Pinus brutia (light green); Qa, Quercus alnifolia (orange); Cb, Cedrus brevifolia (blue); Pn, Pinus nigra (black); Jn, Juniperus foetidissima (purple).Top panel: background map displays elevation at 90-m resolution (SRTM Digital Elevation Data, http://srtm.csi.cgiar.org/).Bottom panel: data from the Department of Forests (Ministry of Agriculture, Rural Development and Environment, Republic of Cyprus; https://www.data.gov.cy/) using the aforementioned DEM (see Appendix S1).Specifically, we interpolated six WorldClim (annual mean temperature, maximum temperature of warmest month, minimum temperature of coldest month, annual precipitation, precipitation of wettest quarter and precipitation of driest quarter; Fick & Hijmans, 2017) and three ENVIREM (climatic moisture index, Thornthwaite aridity index and topographic wetness index; Title & Bemmels, 2018) variables.
both across-habitats and within-habitat scales.Specifically, the response variables were the community dissimilarity matrices based on the Simpson dissimilarity index (β SIM ), and the explanatory variables were forward selected from full models containing the following sets of predictors: (i) forest habitat type (HAB); (ii) spatial variables (SPA PCNMi ) derived from the transformation of the topographic weighted distance matrix using Principal Coordinates of Neighbour Matrices (PCNM, see Appendix S1); and (iii) topoclimatic variables (ENV PC1 and ENV PC2 Figure2).
. Analyses were performed on the Simpson dissimilarity index (β SIM ) community dissimilarity matrices at ASV and OTU levels.Percentages refer to adjusted coefficients of determination (R 2 ADJ ).Percentages below 0.25% are not shown wocDNA metabarcoding enabled us to characterize the soil microarthropod assemblages of the understudied montane forests of Cyprus and provide insights into the drivers of metacommunity structure across the topographically complex region of Troodos.By integrating community metabarcoding with high-resolution topoclimatic and landscape information, we revealed that environmental filtering induced by forest habitat type and topoclimatic heterogeneity controls different facets of Troodos soil biodiversity (α diversity, β diversity and community uniqueness), while habitat connectivity mediates the spatial turnover across the highly fragmented Quercus alnifolia (Qa) habitat.Patterns of β diversity were very similar at OTU and ASV levels, but α diversity varied, with OTU richness following an altitudinal gradient and ASV richness a longitudinal one, probably indicating a decline of genetic TA B L E 3 Multiple matrix regression with randomization (MRR) analyses on community composition (Simpson dissimilarity index, β SIM ) of Quercus alnifolia (Qa) sampling sites at ASV and OTU levels

F I G U R E 5
Distance decay of community similarity across the Quercus alnifolia (Qa) habitat.Distances among Q.alnifolia sampling sites were calculated either based on the isolation-by-resistance (IBR) scenario reflecting the spatial distribution of Qa forest patches (FRA IBR , left panel) or on the topoclimatic distances among sampling sites (ENV PC1-2 , right panel).Community similarity was estimated using the Simpson similarity index (1 − β SIM ) at ASV (filled circles, solid regression line) and OTU (open circles, dashed regression line) levels.
ACK N OWLED G EM ENTS This work is a product of the iBioGen project, which has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 810729.We thank two anonymous referees for their constructive and valuable comments on an earlier version of the manuscript.We are grateful to the Department of Forests-Ministry of Agriculture, Rural Development and Environment, Republic of Cyprus-for providing sampling permission and cartography of the forest habitat distributions.We would like to express our gratitude to Savvas Savva and Heriberto López for assistance in designing and building soil processing facilities at the University of Cyprus, Zoe Makridou, Loudmila J. Lagou, Stylianos Mavrianos and Andreas Dimitriou for help during sampling and sample processing, Konstantinos Ntatsopoulos for support with Coleoptera taxonomy, and Angelina Ceballos-Escalera and Elena Lugli (Natural History Museum, London) for advice during metabarcoding library preparation.We also wish to thank "Centro de Supercomputación de Galicia" (CESGA) for access to computer resources.V.N. was supported by a postdoctoral contract under the iBioGen project (grant agreement No. 810729) and a "Juan de la Cierva-Formación" postdoctoral fellowship (grant FJC2018-035611-I) funded by MCIN/AEI/10.13039/501100011033.AUTH O R CO NTR I B UTI O N S V.N. and A.P. conceived and led the study.V.N. and E.M. performed fieldwork.V.N. performed laboratory work with help from E.M.A.C-I.generated climatic and topographic information and wrote the respective methodological section.V.N. and A.P. conceived the methodological approach and V.N. performed the analyses.V.N. wrote the manuscript with help from A.P. All authors contributed critically to the draft and gave final approval for publication.
) as explanatory variables