Dispersal limitations and long-term persistence drive differentiation from haplotypes to communities within a tropical sky-island: evidence from community metabarcoding

For Review Only Dispersal limitations and long-term persistence drive differentiation from haplotypes to communities within a tropical sky-island: evidence from community metabarcoding Journal: Molecular Ecology Manuscript ID MEC-20-1200.R2 Manuscript Type: Original Article Date Submitted by the Author: 05-Jul-2021 Complete List of Authors: Gálvez-Reyes, Nancy; Universidad Nacional Autonoma de Mexico Instituto de Ecologia, Departamento de Ecología Evolutiva; Universidad Nacional Autónoma de México, Programa de Doctorado en Ciencias Biomédicas Arribas, Paula; IPNA, Island Ecology and Evolution Andujar, Carmelo; IPNA, Island Ecology and Evolution Emerson, Brent; IPNA, Island Ecology and Evolution Piñero, Daniel; Universidad Nacional Autonoma de Mexico Instituto de Ecologia, Departamento de Ecología Evolutiva Mastretta-Yanes, Alicia; CONABIO; CONACYT Keywords: neutral theory, arthropods, Nevado de Toluca, community metabarcoding, tropical mountains Molecular Ecology For Review Only 1 1Article type: Original Article 2 3 4Dispersal limitations and long-term persistence drive differentiation from 5haplotypes to communities within a tropical sky-island: evidence from community 6metabarcoding 7 8Nancy Gálvez-Reyes1,2*, Paula Arribas3, Carmelo Andújar3, Brent C. Emerson3, Daniel 9Piñero1 and Alicia Mastretta-Yanes4,5* 10 11 1. Departamento de Ecología Evolutiva, Instituto de Ecología, Universidad Nacional 12 Autónoma de México, CDMX, México. 13 2. Programa de Doctorado en Ciencias Biomédicas, Universidad Nacional 14 Autónoma de México, CDMX, México. 15 3. Island Ecology and Evolution Research Group, Instituto de Productos Naturales 16 y Agrobiología (IPNA-CSIC), La Laguna, Santa Cruz de Tenerife, Spain. 17 4. Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO), 18 Liga Periférico Insurgentes Sur 4903, Col. Parques del Pedregal, Tlalpan, CDMX, 19 México. 20 5. Consejo Nacional de Ciencia y Tecnología, Benito Juárez (CONACYT), CDMX, 21 México, Avenida Insurgentes Sur 1582, Crédito Constructor, Benito Juárez, 22 CDMX, México. 23 24 25 *Corresponding author: [email protected], [email protected] Page 1 of 45 Molecular Ecology For Review Only 2 26 Abstract 27 Neutral theory proposes that dispersal stochasticity is one of the main drivers of local 28 diversity. Haplotypes-level genetic variation can now be efficiently sampled from across 29 whole communities, thus making it possible to test neutral predictions from the genetic 30 to species-level diversity, and higher. However, empirical data is still limited, with the 31 few studies to date coming from temperate latitudes. Here, we focus on a tropical 32 mountain within the Transmexican Volcanic Belt to evaluate spatially fine-scale patterns 33 of arthropod community assembly to understand the role of dispersal limitation and 34 landscape features as drivers of diversity. We sampled whole-communities of arthropods 35 for eight orders at a spatial scale ranging from 50 m to 19 km, using whole community 36 metabarcoding. We explored multiple hierarchical levels, from individual haplotypes to 37 lineages at 0.5, 1.5, 3, 5, 7.5% similarity thresholds, to evaluate patterns of richness, 38 turnover, and distance decay of similarity with isolation-by-distance and isolation-by39 resistance (costs to dispersal given by landscape features) approaches. Our results showed 40 that distance and altitude influence distance decay of similarity at all hierarchical levels. 41 This holds for arthropod groups of contrasting dispersal abilities, but with different 42 strength depending on the spatial scale. Our results support a model where local-scale 43 differentiation mediated by dispersal constraints, combined with long-term persistence of 44 lineages, is an important driver of diversity within tropical sky islands. 45 46 47 48 Keywords: neutral theory, arthropods, Nevado de Toluca, community metabarcoding, 49 tropical mountains 50 Page 2 of 45Molecular Ecology For Review Only 3 51 Introduction 52 Biodiversity is a hierarchical concept, embracing levels of diversity from genes to 53 ecosystems (Bailey et al., 2009). However, the spatial structuring of biodiversity has 54 largely been addressed considering each level independently (e.g., genes, species, 55 communities), and at disparate spatial and temporal scales. To understand the distribution 56 of diversity within an integrated spatio-temporal framework, neutral theories of diversity 57 propose that stochastic, but spatially limited dispersal, may be the primary driver of local 58 diversity (Bell, 2001; Vellend, 2010; Rosindell et al., 2011), but patterns have been 59 described for each of the levels of biodiversity in an isolated fashion. Similarly, in the 60 neutral theory of molecular evolution (Kimura 1968, 1969), changes in genomes may 61 accumulate passively through time via genetic drift in alleles that are selectively neutral. 62 Both of these neutral theories have been integrated by linking stochastic population and 63 species-level processes under common assumptions of migration, speciation/mutation 64 and extinction (Vellend & Gerber 2005; Papadopoulou et al., 2011). However, 65 incorporating empirical data to test neutrality predictions requires sequencing entire 66 communities, which has only recently become feasible. The few empirical studies to date 67 (Baselga et al., 2013; 2015; Gómez-Rodríguez & Baselga, 2018; Arribas et al., 2020) 68 have found support for the generality of neutral processes across different hierarchical 69 levels of diversity. The generality of these observed patterns across geographical scales, 70 taxa and regions of the world remains to be assessed (Baselga et al., 2013; Arribas et al., 71 2020). 72 To examine neutrality predictions at the different levels of biodiversity, it has been 73 proposed that neutral processes (including neutral mutation, dispersal limitation, and and 74 birth and death of lineages) should act uniformly across hierarchical levels (that is from 75 haplotypes to species and higher taxonomic levels), with non-neutral processes differing Page 3 of 45 Molecular Ecology For Review Only 4 76 among levels (Baselga et al., 2013). This can be tested by analyzing the typically neutral 77 (but see Gerber et al., 2001) haplotype variation of the mitochondrial COI gene, which, 78 when assessed for entire communities, results in a largely regular decay of similarity with 79 spatial distance (Baselga et al., 2013). Under a scenario where dispersal constraints are a 80 dominant driver of spatial variation in community structure, assemblage turnover at the 81 species level should mirror, in a fractal way, distance-decay of haplotype-level similarity 82 (Baselga et al., 2013; Baselga, Gómez-Rodríguez, & Vogler, 2015). 83 To cross-fertilize predictions from the neutral theory of molecular evolution and 84 the neutral theory of biodiversity as proposed above, genetic data sampled across entire 85 communities are required. Whole community metabarcoding (cMBC) with bulk 86 sequencing of the mitochondrial COI gene is particularly promising for accomplishing 87 such a task (Andújar et al., 2018). Recent improvements for denoising metabarcoding 88 datasets (e.g., Edgar, 2016; Callahan et al., 2016) and evaluating the prevalence of 89 sequencing errors and co-amplified pseudogenes (Andújar et al., 2021) raise the prospect 90 of read-based, haplotype-level analyses with mitochondrial COI cMBC data, representing 91 a step change for the study of diversity patterns through whole-community genetic 92 analyses (Andújar et al., 2018; Arribas et al., 2020). Haplotype data from communities 93 can be used directly for analyses of genetic diversity or aggregated into species-level 94 entities for analyses of species diversity. This allows for the joint analysis of turnover 95 (beta diversity) at multiple hierarchical levels, based on analyses of distance decay of 96 similarity across lineages (established at different clustering thresholds). This analytical 97 approach has been exploited to determine whether the composition across multiple beetle 98 taxa assemblages is predominantly driven by dispersal (Baselga et al., 2013; Baselga, 99 Gómez-Rodríguez, & Vogler, 2015) and has been proposed as a useful way to compare 100 relative dispersal constraints between lineages from different taxonomic groups (GómezPage 4 of 45Molecular Ecology For Review Only 5 101 Rodríguez et al., 2019, Múrria et al., 2017). These and other studies using a multi102 hierarchical approach (e.g., Baselga et al., 2013, 2015; Gómez-Rodríguez et al., 2019), 103 have focused from regional to continental-scale distances. However recent work has also 104 exploited this framework to analyze community assemblage at finer spatial scales (<15 105 km) using cMBC data (Arribas et al., 2020). Thus, the reliable haplotype information 106 from mitochondrial COI cMBC approach seems well suited to test for neutrality 107 predictions from genes to communities across geographic scales, taxa, and biomes around 108 the world, in particular in the tropics which have yet not been studied with this approach. 109 Here, we focus on the potential role of dispersal limitations within tropical 110 mountain regions, which are characterized as hyperdiverse (Fjeldså et al., 2012; Rahbek, 111 et al., 2019a; Rahbek, et al., 2019b). In addition to the patchy distribution of habitat 112 resulting from topoclimate heterogeneity, two additional factors may also be important 113 determinants of this diversity. First, because tropical taxa are expected to have narrow 114 physiological tolerances to temperature, tropical mountain passes are more effective 115 barriers to dispersal than those in temperate regions, thus promoting speciation through 116 physical disruption of gene flow (Janzen, 1967; Polato et al., 2018; Sheldon et al., 2018). 117 Second, coupling altitudinal gradients with tropical latitudes allows populations within 118 tropical mountains to persist relatively in situ despite global climate fluctuations 119 (Mastretta-Yanes, et al., 2018; Rahbek, et al., 2019a; Rahbek, et al., 2019b). These 120 processes of isolation and long-term persistence have been widely used to explain 121 diversification among tropical mountain peaks across the world (e.g., Fjeldså et al., 2012; 122 He et al., 2019; Knowles, 2001; Mastretta-Yanes et al., 2018; McCormack et al., 2009; 123 Uscanga et al., 2021). However, recent evidence suggests that dispersal limitation within 124 a single mountain could also play an important role in generating endemism (Bray & 125 Bocak, 2016), but there is a lack of empirical data to test if this is indeed the case. Page 5 of 45 Molecular Ecology For Review Only 6 126 Studying evolutionary processes in mountain ranges across fine spatial scales is 127 challenging due to their complexity with regard to topography and geological and climatic 128 history. However, the insular nature of sky-islands reduces the complexity of 129 mountainous regions, thus providing more simplified systems to analyze evolutionary 130 processes. Within true oceanic islands, habitat discontinuity has been demonstrated to 131 have an important role in driving geographical diversification (García‐Olivares et al., 132 2019; Goodman et al., 2012). Additionally, Salces‐Castellano et al. (2020) have 133 demonstrated that, when dispersal ability and climate tolerance are restricted, strong 134 geographic isolation over distances of only a few kilometers can be found for multiple 135 co-occurring arthropod species of an oceanic island. These results support that in addition 136 to allopatric diversification between different islands, island systems could also promote 137 high levels of intra-island geographical diversification, acting as a local-scale diversity 138 source. Arthropod communities are ideal systems to test evolutionary processes at fine 139 spatial scales within sky-islands because they are often locally abundant and diverse and 140 relatively easy to sample in bulk. They also harbour a broad diversity of groups, including 141 winged and non-winged species and varying body-sizes, thus allowing to compare taxa 142 with different dispersal abilities. 143 We hypothesise that dispersal limitations within a single mountain, coupled with 144 long-term (glacial/interglacial) population persistence, play an important role in 145 generating tropical mountain diversity, from haplotypes to communities. Montane 146 ecosystems within the sky islands of the Transmexican Volcanic Belt (TMVB), including 147 the Nevado de Toluca study site, have been shown to persist locally during climate 148 fluctuations (Mastretta-Yanes et al., 2018). Within this framework we evaluate the roles 149 of dispersal limitation and landscape features as drivers of differentiation within a single 150 sky-island, by examining fine-scale community patterns of arthropod fauna at different Page 6 of 45Molecular Ecology For Review Only 7 151 hierarchical levels. To do this, we performed a systematic sampling consisting of 840 152 pitfall traps across Abies religiosa forests within Nevado de Toluca, a sky-island from the 153 TMVB. We generated haplotype-level metabarcoding data for 42 arthropod communities 154 distributed in sampling blocks separated from 50 m to 19 km, and evaluate patterns of 155 richness, turnover, and distance decay in community similarity at multiple hierarchical 156 levels (haplotype, putative species and supra-specific levels). As we are interested in the 157 effect of ecological and topographic features on dispersal limitation, in addition to testing 158 for the effect of isolation-by-distance (IBD), we also performed isolation-by-resistance 159 (IBR) analyses, incorporating costs to dispersal for landscape features such as altitude or 160 habitat type (McRae, 2006). This yields biologically more informative distance decay 161 relationships than Euclidean distance alone (McRae, 2006). Our analytical framework 162 thus allows us to assess the role of dispersal constraints, within a local landscape setting, 163 on tropical mountain diversity. 164 165 Materials and methods 166 2.1. Study area and bulk sampling 167 The study area comprises Abies religiosa forests which grow at around 2800-3500 m.a.s.l. 168 in the Nevado de Toluca volcano, which is a sky-island in the TMBV (Mastretta-Yanes 169 et al., 2015; Rzendowski, 2006). We sampled arthropods from 14 sampling points 170 distributed across four sites: Tlacotepec (TLC; 3 sampling points), San Juan de las 171 Huertas (SJH; 3 sampling points), San Bartolo (ASB; 5 sampling points) and Agua 172 Bendita (AAB; 3 sampling points; Figure 1), all within the conservation zone of the 173 natural protected area of Nevado de Toluca (Figure 1). Sampling was performed during 174 the rainy season of 2015 during mid-August and September. Pitfall traps were active for 175 15 days, and both the placement and collection of all traps was undertaken within a twoPage 7 of 45 Molecular Ecology For Review Only 8 176 day period to reduce stochastic heterogeneity among samples due to temporal climatic 177 variation. At each of the 14 sampling points, three sampling blocks of 20 x 15 m were 178 established, resulting in a total of 42 community samples (Figure 1). Each sample 179 consisted of all the specimens collected by 20 pitfall traps distributed equidistantly inside 180 each sampling totalling 840 pitfall traps. Sampling blocks were separated by at least 50 181 m within each sampling point, and maximum distance among sampling points was 19 km 182 (Figure 1). Pitfall traps consisted of plastic cups of 13 cm height and 10 cm diameter, 183 with five 2 x 1.5 cm perforations, 2 cm apart from each other, at a height of 10 cm above 184 the base. Cups were buried to the height of the window and lids were added to prevent 185 the entry of rainwater. Cups were filled with a mixture of 185 ml of ethanol 70% and 15 186 ml of glycerin (Figure 1b; Figure S1a), which we previously tested was enough to prevent 187 evaporation and thus DNA degradation. Pitfall traps of each block site (20 traps) were 188 collected after 15 days, as suggested by Cardoso et al. (2008, 2009) and Emerson et al. 189 (2017) and pooled in a single bulk-sample in a bottle containing ethanol at 96%. Sampling 190 was performed with SEMARNAT permit No. SGPA/DGVS/02641/15. 191 192 2.2 Molecular laboratory processes 193 We cleaned up each sample (comprising 20 pooled pitfalls) following a Flotation194 Filtration-Stereoscope protocol (FFS) that allowed us to have ‘clean’ bulk specimens 195 ready for DNA extraction (see Figure S1b-g). The samples had considerably different 196 body sizes, so to prevent differences in biomass between specimens from causing biases 197 at the DNA extraction and PCR steps (Elbrecht & Leese, 2015), we divided specimens 198 by size class boundaries into: small (e.g., size Drosophila melanogaster <=3 mm), 199 medium (e.g., size adult Apis mellifera > 3mm and <=15 mm), or large specimens (e.g., 200 adult grasshopper >15 mm). We then divided each sample as follows: complete bodies Page 8 of 45Molecular Ecology For Review Only 15 350 This is equivalent to Euclidean distance but accounts for the finite size of the input 351 landscape being analyzed, and is thus more appropriate for comparison (Lee-Yaw et al., 352 2009). 353 For the vegetation type rasters, a high resolution digital land cover map of the 354 Nevado of Toluca was used, based on satellite images (SPOT 6/7) of 2015 defining the 355 Abies religiosa forest (sampled vegetation), Pinus hartwegii forest, alpine grassland, 356 agriculture, urbanization, water in the crater and Nevado de Toluca’s crater (González357 Fernández et al., 2019). After cross-validation with our sampling records, minor 358 modifications were performed to adjust Abies forest distribution in areas of high 359 topographic complexity (see Supporting Material S2 for details). Conductance grids were 360 then constructed for the vegetation heterogeneity analyses assigning different 361 conductance values to each vegetation type (maximum value of conductance is 1, 362 meaning no resistance to dispersion) as detailed in Table S2. 363 To build the conductance grids for altitude and slope, a similar approach was 364 followed, assigning different conductance values to altitudinal and slope ranges. In total 365 30 conductance grids were tested (Table S2; Figure S5). Additionally, distance decay of 366 similarity was tested for at smaller geographic distances, for which West (AAB, ASB) 367 and East (SJH, TLC) sites were analyzed separately, but using only the “flat” raster and 368 the raster from the previous analysis with the highest explanatory power. 369 A negative exponential function ‘decay.model’ in the R package betapart 370 (Baselga & Orme, 2012) was used to adjust a negative exponential function to a 371 generalized linear model (GML). We used Simpson similarity (1 – βsim; Baselga, 2010) 372 as a response variable, pairwise effective distances of each resistance surface as 373 predictors, log link and Gaussian error (Arribas et al., 2020; Gómez-Rodríguez & 374 Baselga, 2018). Finally, we evaluated the fractal pattern (i.e., self‐similar systems; Page 15 of 45 Molecular Ecology For Review Only 16 375 Baselga et al., 2015) by a log–log Pearson correlation of the haplotype level and, 376 independently: (1) number of lineages, (2) initial similarity (i.e., intercept), and (3) mean 377 similarity. This was undertaken for each of the six orders across the sampling area, and 378 for Collembola and Diptera at the East and West sections. High correlation values are 379 indicative of self-similarity in lineage branching (i.e., number of lineages) and/or spatial 380 geometry of lineage distributional ranges (i.e., initial and mean similarity; Baselga et al., 381 2015), supporting the presence of fractal patterns, which are predicted under a neutral 382 process of community evolution (Arribas et al., 2020). Analyses and graphical 383 representations of data were performed with the R packages vegan (Oksanen et al., 2019) 384 betapart (Baselga et al., 2018) and ecodist (Goslee & Urban, 2020). 385 386 Results 387 Phylogenetic groups and ASVs recovered by COI 388 Across the 42 libraries, MiSeq sequencing generated a total of 11,639,999 paired reads. 389 We obtained from 108,583 to 419,903 reads in each direction for each sample. Of these, 390 from 39,484 to 178,720 sequences remained after quality filtering (totalling 4,990,334). 391 After read merging and sequence filtering to a length of 418 bp, each sample comprised 392 between 33,609 and 155,634 sequences (totalling 4,305,390). Taxonomic assignments 393 with usearch showed high similarity to a broad range of arthropod species (Figure 2a). 394 The 4,305,390 sequences comprised 1,277 ASVs (unique variants: Diptera, 385; 395 Collembola, 270; Arachnida, 155; Coleoptera, 136; Hymenoptera, 133; Hemiptera, 116; 396 Myriapoda, 51; Lepidoptera, 31) (Figure 2b). The GMYC threshold values obtained were: 397 Diptera, 0.9%; Collembola, 2.9%; Arachnida, 1.3%; Coleoptera, 0.7%; Hymenoptera, 398 1%; Hemiptera, 1.8%; Myriapoda, 1.5%; Lepidoptera 1.5% (Table S1). 399 Page 16 of 45Molecular Ecology For Review Only 17 400 Arthropod ASVs richness at multi-hierarchical levels 401 Rarefaction curves indicate that sampling effort was sufficient to achieve from 60 to 91% 402 completeness at the levels of individual haplotypes and lineages at 3%, for four of the 403 eight taxa (Diptera, Collembola, Arachnida and Hymenoptera Figure S3). The groups that 404 yielded higher richness within communities (alpha diversity by sample) were Diptera 405 (mean = 87 SD ± 23.1 haplotypes and mean = 21 ± 5.15 lineages at 3%), Collembola 406 (mean = 38 ± 9.44 haplotypes and mean = 16 ± 2.91 lineages at 3%), and Arachnida 407 (mean = 12 ± 6.28 haplotypes and mean = 9 ± 3.64 lineages at 3%). Diptera, Collembola, 408 Arachnida, Hymenoptera and Coleoptera presented significant differences in richness 409 among sampling sites but patterns were not consistent across groups (Figure 3). 410 411 Community composition at multi-hierarchical levels 412 The overall turnover (βsim) among sampling points was high for all the arthropod groups 413 across all levels analysed (close to 1, βsim haplotypes from 0.904 to 0.957; βsim 3% from 414 0.896 to 0.945; Figure S4) and average pairwise βsne was close to zero, βsne haplotypes 415 from 0.010 to 0.026; βsne 3% from 0.014 to 0.033. The mean value for βsim for haplotypes 416 of Coleoptera (n=136) was 0.957, the average pairwise βsne was 0.014, and βtotal was 0.970 417 (βsor). For Diptera (n=385), the mean value for βsim was 0.904, the average pairwise βsne 418 was 0.024, and the average pairwise βsor was 0.928. 419 NMDS ordination plots revealed differences in community composition among 420 sampling sites, particularly among those located in the West (AAB and ASB) and East 421 (SJH and TLC) sampling points within Nevado de Toluca (Figure 4). ANOSIM revealed 422 that differences of dissimilarity were significant for the different groups and lineages at 3 423 and 5%, except on Myriapoda (Figure 4). The orders Diptera, Collembola, Hemiptera and 424 Coleoptera each formed two geographic groups: (i) AAB-ASB (West) and (ii) SJH-TLC Page 17 of 45 Molecular Ecology For Review Only 18 425 (East). In contrast, three groups were recovered within Arachnida: (i) AAB, (ii) ASB and 426 (iii) SJH-TLC (Figure 4). The variation in Collembola, Arachnida and Diptera among 427 sites largely contributed to the ordinations of community similarity among sites. The 428 highest values were observed in Collembola (Haplotype r2 = 0.914, p < 0.001; lineages at 429 3% r2 = 0.826, p < 0.001), followed by Diptera (Haplotype r2 = 0.595, p < 0.001; lineages 430 at 3% r2 = 0.372, p < 0.001) and Arachnida (Haplotype r2 = 0.569, p < 0.001; lineages at 431 3% r2 = 0.452, p < 0.001; Figure 4). 432 433 Landscape community connectivity 434 Across all sampling sites (max distance 19 km), pairwise similarity within communities 435 at all levels of clustering decreased with Euclidean (“flat”) distance (Figure 5ac, Figure 436 S7). The fit of IBD was higher in Collembola (from r2 = 0.704, b = -2.07, p < 0.001 at the 437 haplotype level, to r2 = 0.580, b = -0.92, p < 0.001 at 7.5% lineages; Table S3; Figure 5a) 438 than Diptera (from r2 = 0.293, b = -0.56, p < 0.001 at the haplotype level, to r2 = 0.036, b 439 = -0.14, p < 0.001 lineages at 7.5%; Table S4; Figure 5c). Similar IBD results were found 440 for Arachnida, Coleoptera, Hemiptera and Hymenoptera (Figure S7, Table S7). However, 441 variation was better explained by IBR (Figure 5bd; Figure S5). For Collembola, the 442 explanatory power of the resistance surface “Altitude 3,000” was slightly higher than the 443 “flat surface” at different lineages (from r2 = 0.723, b = -1.63, p < 0.001 at the haplotype 444 level, b = -0.71, p < 0.001 at 7.5% lineages; Table S3; Figure 5b). In Diptera the highest 445 explanatory power was given by the resistance surface “Mean-elevation peak 446 (symmetrical)” (from r2 = 0.286, b = -0.15, p < 0.001 at the haplotype level, b = -0.05, p 447 < 0.001 at 7.5% lineages; Table S4, Figure 5d). 448 Distance-decay relationships (DDR) at a finer geographic scale within the West 449 and East regions remained strong in Collembola, but were weaker or non-significant in Page 18 of 45Molecular Ecology For Review Only 19 450 Diptera (Figure 6). East sampling points (<5 km), also showed distance decay of 451 similarity in Collembola (Table S5, Figure 6a), and Diptera, but coefficients were lower 452 in the second (Table S5, Figure 6c). When considering IBR analysis, very similar results 453 were found (Table S5, Figure S6ac). In the West section of our sampling (<2 km), a 454 similar pattern to the East was found for Collembola (Table S5, Figure 6b) but Diptera 455 showed non-significant results both for the IBD and IBR analyses (Table S5, Figure 6d, 456 Figure S6d). Exponential decay curves yielded higher slopes for Collembola than Diptera. 457 All lineages of Diptera and Collembola show the same pattern of distance decay, but each 458 lineage shows equal or greater pairwise similarity than at the previous one (Figure 5,6; 459 Table S3, S4). Regarding the multi-hierarchical analysis, the distance decay at the 460 haplotype level showed a significant log–log correlation with the number of lineages, 461 initial similarity, and mean similarity of communities. This is true for both our entire 462 sampling area and focusing on the West and East areas (Table S6). The log-log linear 463 correlations suggest that the patterns of assemblage variation across hierarchical levels 464 can be described by a fractal geometry (Baselga et al., 2013, 2015). Thus, we found that 465 community variation across genetic similarity levels can be described by a fractal 466 geometry in each group in all the geographic scales sampled. 467 468 Discussion 469 We recovered 1,277 ASVs, 496 lineages at 3%, 441 at 5% and 570 using GMYC across 470 the eight arthropod orders. Composition of arthropod communities exhibited high 471 turnover among sampling blocks and ANOSIM tests also showed significant differences 472 among sites within blocks. Across all sampling points, spanning a maximum distance of 473 only 19 km, we found that distance (pure geographical distance and corrected by 474 elevation) is a key variable explaining community structure, from the level of haplotypes Page 19 of 45 Molecular Ecology For Review Only 20 475 through to lineages. This pattern holds at finer geographic distances (<2 km), but only on 476 it up with considerably low dispersal ability (non-winged Collembola). 477 478 Arthropod communities recovered by metabarcoding and sampling blocks 479 The estimation of species richness in any ecological setting, and especially in forested 480 environments, can be challenging due to the rarity of some species, differences in 481 detection probabilities, and the field effort necessary to collect enough samples or species 482 to ensure meaningful coverage (Andújar et al., 2017; Arribas et al., 2016; Creedy et al., 483 2019). In our study, we used pitfall traps in sampling blocks, maximizing the probability 484 of detecting arthropod species by sampling intensively at multiple sites in one mountain, 485 covering eight orders of arthropods, and characterising these from haplotypes to 486 communities. Our sampling method and size sorting step allowed the recovery of eight 487 major arthropod orders, congruently with other metabarcoding analyses (Elbrecht et al., 488 2017; Elbrecht et al., 2018; Creedy et al., 2019). According to our rarefaction curves, our 489 sampling detected different taxa per sample (Figure S3), demonstrating the utility of 490 cMBC and our sampling design to study a region of high biological diversity and 491 ecological complexity. It is difficult to compare our results against morphological studies 492 in Nevado de Toluca because there are no complete checklists of arthropods for Mexican 493 highlands. For instance, out of 29 dung beetle species found by a recent survey of four 494 sky-islands, more than 10% were new species (Arriaga-Jiménez et al., 2018). 495 The metabarcoding pipeline we implemented (Arribas et al., 2020) allowed us to 496 analyze diversity patterns for each order separately and to analyze community assembly 497 within a multi-hierarchical framework, from individual haplotypes to lineage divergence 498 at 7.5%. Our approach identified 1,277 ASVs from 42 bulk samples, allowing us to 499 estimate community composition and turnover without the bias introduced by traditional Page 20 of 45Molecular Ecology For Review Only 21 500 taxonomy (Creedy et al., 2019), including small taxa (< 0.5 mm), such as Collembola, 501 and specimens that break easily, such as Diptera (Figure 3). Although the cMBC method 502 allowed us to explore a large dataset of hundreds of pitfall traps, there are potential 503 limitations which may introduce bias such as the use of universal primers for PCR (Piñol 504 et al., 2014), biomass differences (Elbrecht, Peinert, & Leese, 2017), and/or tag jumping 505 (Schnell, Bohmann, & Gilbert, 2015). We have attempted to reduce potential bias by 506 conducting PCRs in triplicate, dividing our specimens by size classes, and by using HLD 507 desalted primers to decrease tag jumping. 508 509 Strong turnover of arthropods communities at multi-hierarchical levels 510 Our study provides haplotype-level and higher lineages data of entire communities 511 (Figure 4, Figure S4). High beta diversity was found across communities in all orders, 512 and was dominated by lineages turnover (βsim) instead of nestedness (βsne), from 513 haplotypes to higher hierarchical levels (Figure S4). We found significant differentiation 514 among sites and vegetation type (Abies forests) according to site haplotype composition, 515 and this differentiation persisted when molecular entities that conservatively represent 516 species were considered. Diptera and Collembola presented the highest differentiation 517 among sites, specifically dividing the East (SJH-TLC) and the West (AAB-ASB) 518 sampling points, as shown by the NMDS and ANOSIM analyses (Figure 4). Again, this 519 occurs from the level of haplotype to higher lineages (Figure 4). The West-East division 520 in community structure within Nevado de Toluca (Figure 4) corresponds to opposing 521 hillsides. In mountain landscapes, East-facing slopes with morning sun may provide 522 different conditions from cold and foggy West-facing slopes (Rahbek et al., 2019a), 523 which could influence community assembly. While we cannot discard the role of 524 environmental heterogeneity (environmental distance) driving community differences Page 21 of 45 Molecular Ecology For Review Only 22 525 between western-eastern sampling points, it should be emphasized that all sampling 526 points were located within A. religiosa forests with no apparent differentiation within the 527 Nevado de Toluca. In fact, A. religiosa forests are expected to grow under similar 528 conditions within a single mountain, having a quite restricted environmental niche 529 associated with moist and cold sites (Rzendowski, 2006). 530 531 Dispersal limitations drive community structure across multi-hierarchical levels at fine 532 geographical scales 533 To investigate whether dispersal limitations shape distance-decay within a mountain, we 534 examined DDR at multi-hierarchical levels, including haplotypes which are expected to 535 behave neutrally, even across environmental gradients. IBD analyses were conducted for 536 six orders, and IBR analyses were conducted for Diptera and Collembola (which have 537 contrasting dispersal abilities) using 30 conductance matrices. For Diptera and 538 Collembola, we considered our entire sampling in Nevado de Toluca (19 km max distance 539 among sampling points) and finer geographic scales within the East (<5 km) and West 540 (<2 km) subsets of our sampling. We found that decay of similarity of communities 541 decreases with spatial distance at the level of haplotypes, at all levels of lineage 542 divergence (from 0.5 to 7.5%), and when haplotypes were segregated to putative 543 molecular species using GMYC (Figure 5). Distance decay was high in all orders (Figure 544 6, Figure S7), but it is more marked in groups with more limited dispersal abilities (e.g., 545 Collembola and Arachnida) than in groups with higher dispersal abilities (e.g., Diptera, 546 Coleoptera, Hemiptera and Hymenoptera). Interestingly, for Collembola our results hold 547 both considering our entire sampling as well as finer (<2 km) geographic distances 548 (Figure 5, 6), which is consistent with genetic studies within Collembola showing genetic 549 differentiation over very short geographic distances (Cicconardi et al., 2013; Arribas et Page 22 of 45Molecular Ecology For Review Only 23 550 al., 2020). 551 High dispersal ability is expected to enhance community similarity (Baselga et al., 552 2012). Our results support this, as distance decay is higher in the wingless Collembola (r2 553 = 0.704 and r2 = 0.599 at the haplotype and GMYC levels, respectively; Table S3; Figure 554 5a) than in the winged Diptera (r2 = 0.293 and r2 = 0.195 at the haplotype and GMYC 555 levels, respectively; Table S4; Figure 5c). Similar patterns of higher distance decay 556 relationships at multi-hierarchical levels in poorly dispersing organisms than in better 557 dispersers have been found in European water beetles (Baselga et al., 2013), Iberian leaf 558 beetles (Baselga et al., 2015) and European beetles (Gómez-Rodríguez & Baselga, 2018) 559 at much larger (hundreds of km) geographical scales than here (but see Gómez-Rodríguez 560 et al., 2019 where the pattern was not clear for terrestrial molluscs), and only recently 561 even at the scale of a few kilometers, in Collembola, Acari and Coleoptera (Arribas et al., 562 2020). Communities of good dispersers are more homogeneous not only because they can 563 disperse larger distances, but also because they can more easily overcome geographical 564 barriers between suitable habitat (Thompson & Townsend, 2006; Vellend, 2010). Thus, 565 as dispersal abilities become limited, landscape features impeding dispersal are expected 566 to play a more important role in structuring diversity, which can be explicitly tested 567 including landscape features in analyses such as IBR. 568 IBR quantifies ‘effective distances’ between communities that may yield more 569 biologically informative DDR than Euclidean distance (McRae, 2006; McRae et al., 570 2008). We found that altitudinal differences better explain similarity decay than distance 571 alone (“flat” landscape), slope or vegetation type. The resistance surface “flat” (i.e., IBD) 572 has slightly less explanatory power for Collembola (r2 = 0.704, p < 0.001 at the haplotype 573 level; Table S3; Figure 5a) than “Altitude 3,000” (r2 = 0.723, p < 0.001 at the haplotype 574 level; Table S3; Figure 5b), the best fitting resistance surface. Altitude 3,000 resistance Page 23 of 45 Molecular Ecology For Review Only 24 575 surface corresponds to the elevation at which the Nevado de Toluca volcano massif begins 576 (Figure S5; Table S2), suggesting that Collembola followed a pattern of IBD and that 577 within their limited dispersal ability, landscape features do not represent an impediment. 578 For Diptera, the highest explanatory power was provided by the resistance surface “Mean579 elevation peak (symmetrical)” (r2 = 0.286, p < 0.001 at the haplotype level; Table S4; 580 Figure 5d). This resistance surface assumes maximum conductance at the mean altitude 581 of our sampling and a gradual decrease until reaching altitudes outside of our sampling 582 range, but still where Abies forest can be found (Table S2). This suggests that even though 583 they are good dispersers, the efficiency of Diptera movement through relatively 584 unsuitable conditions (different altitudes) is compromised. We thus find that for Diptera, 585 it is not distance alone that drives community structure, but also landscape features. 586 Although our sampling blocks are separated by short distances from 50 m to 19 km, 587 effective distances among sites depends upon the elevation model used to set the 588 conductance values (Table S2; Figure S5). This is consistent with Janzen’s prediction of 589 “mountain passes being higher in the tropics” (Janzen, 1967), and adds to the recent 590 empirical data corroborating it (Polato et al., 2018). However, although our results show 591 that landscape connectivity contributes to dispersal limitation, geographic distance seems 592 to play a more dominant role both for both orders. This is consistent with dispersal 593 limitation acting over evolutionary time, as has been suggested to explain the small spatial 594 scale diversification of Scarelus beetles within tropical mountains (Bray & Bocak, 2016). 595 Distance decay patterns at the species level could reflect spatially correlated 596 environmental heterogeneity (i.e., between western and eastern sides of the Nevado de 597 Toluca). While some degree of environmental distance could explain biodiversity 598 patterns (but see above on the homogeneity of the sampling study habitat), the following 599 findings point to dispersal limitation within this single single sky-island as a major driver Page 24 of 45Molecular Ecology For Review Only 31 754 R. (2018). betapart: Partitioning Beta Diversity into Turnover and Nestedness 755 Components (1.5.1) [Computer software]. https://CRAN.R756 project.org/package=betapart 757 Bell, G. (2001). Neutral macroecology. Science, 293(5539), 2413-2418. https://doi.org/ 758 10.1126/science.293.5539.2413 759 Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina 760 sequence data. Bioinformatics, 30(15), 2114–2120. 761 https://doi.org/10.1093/bioinformatics/btu170 762 Bray, T. C., & Bocak, L. (2016). Slowly dispersing neotenic beetles can speciate on a penny 763 coin and generate space-limited diversity in the tropical mountains. Scientific Reports, 764 6(1), 33579. https://doi.org/10.1038/srep33579 765 Callahan, B. J., McMurdie, P. J., & Holmes, S. P. (2017). Exact sequence variants should 766 replace operational taxonomic units in marker-gene data analysis. The ISME Journal, 767 11(12), 2639–2643. https://doi.org/10.1038/ismej.2017.119 768 Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., & Holmes, S. P. 769 (2016). DADA2: High-resolution sample inference from Illumina amplicon data. 770 Nature Methods, 13(7), 581–583. https://doi:10.1038/nmeth.3869 771 Cicconardi, F., Fanciulli, P. P., & Emerson, B. C. (2013). Collembola, the biological species 772 concept and the underestimation of global species richness. Molecular Ecology, 22(21), 773 5382–5396. https://doi.org/10.1111/mec.12472 774 Clarke, K. R. (1993). Non-parametric multivariate analyses of changes in community structure. 775 Australian Journal of Ecology, 18(1), 117–143. https://doi.org/10.1111/j.1442776 9993.1993.tb00438.x 777 Creedy, T. J., Ng, W. S., & Vogler, A. P. (2019). Toward accurate species-level metabarcoding 778 of arthropod communities from the tropical forest canopy. Ecology and Evolution, 9(6), 779 3105–3116. https://doi.org/10.1002/ece3.4839 780 Edgar, R. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. 781 Nature methods, 10(10), 996-998. https://doi. org/10.1038/nmeth.2604 Page 31 of 45 Molecular Ecology For Review Only 32 782 Edgar, R. C. (2016). UNOISE2: Improved error-correction for Illumina 16S and ITS amplicon 783 sequencing. BioRxiv, 81257. https://doi.org/10.1101/081257 784 Edgar, R. C., & Flyvbjerg, H. (2015). Error filtering, pair assembly and error correction for 785 next-generation sequencing reads. Bioinformatics (Oxford, England), 31(21), 3476– 786 3482. https://doi.org/10.1093/bioinformatics/btv401 787 Elbrecht, V., & Leese, F. (2015). Can DNA-based ecosystem assessments quantify species 788 abundance? Testing primer bias and biomass—sequence relationships with an 789 innovative metabarcoding protocol. PloS one, 10(7), e0130324. 790 https://doi.org/10.1371/journal.pone.0130324 791 Elbrecht, V., Peinert, B., & Leese, F. (2017). Sorting things out: Assessing effects of unequal 792 specimen biomass on DNA metabarcoding. Ecology and evolution, 7(17), 6918-6926. 793 https://doi.org/10.1002/ece3.3192 794 Elbrecht, V., Vamos, E. E., Steinke, D., & Leese, F. (2018). Estimating intraspecific genetic 795 diversity from community DNA metabarcoding data. PeerJ, 6, e4644. 796 https://doi.org/10.7717/peerj.4644 797 Faria, C. M. A., Shaw, P., & Emerson, B. C. (2019). Evidence for the Pleistocene persistence of 798 Collembola in Great Britain. Journal of Biogeography, jbi.13610. 799 https://doi.org/10.1111/jbi.13610 800 Favreau, J. M., Drew, C. A., Hess, G. R., Rubino, M. J., Koch, F. H., & Eschelbach, K. A. 801 (2006). Recommendations for assessing the effectiveness of surrogate species 802 approaches. Biodiversity & Conservation, 15(12), 3949-3969. https://10.1007/s10531803 005-2631-1 804 Fjeldså, J., Bowie, R. C. K., & Rahbek, C. (2012). The role of mountain ranges in the 805 diversification of birds. Annual Review of Ecology, Evolution, and Systematics, 43(1), 806 249–265. https://doi.org/10.1146/annurev-ecolsys-102710-145113 807 García‐Olivares, V., Patiño, J., Overcast, I., Salces‐Castellano, A., López de Heredia, U., 808 Mora‐Márquez, F., Machado, A., Hickerson, M. J., & Emerson, B. C. (2019). A 809 topoclimate model for Quaternary insular speciation. Journal of Biogeography, 46(12), Page 32 of 45 Molecular Ecology For Review Only 33 810 2769–2786. https://doi.org/10.1111/jbi.13689 811 Gerber, A. S., Loggins, R., Kumar, S., & Dowling, T. E. (2001). Does nonneutral evolution 812 shape observed patterns of DNA variation in animal mitochondrial genomes?. Annual 813 review of genetics, 35(1), 539-566. 814 https://doi.org/10.1146/annurev.genet.35.102401.091106 815 Gómez-Rodríguez, C., & Baselga, A. (2018). Variation among European beetle taxa in patterns 816 of distance decay of similarity suggests a major role of dispersal processes. Ecography, 817 1825–1834. https://doi.org/10.1111/ecog.03693 818 Gómez‐Rodríguez, C., Miller, K. E., Castillejo, J., Iglesias‐Piñeiro, J., & Baselga, A. (2019). 819 Understanding dispersal limitation through the assessment of diversity patterns across 820 phylogenetic scales below the species level. Global Ecology and Biogeography, 28(3), 821 353-364. https://doi.org/10.1111/geb.12857 822 González-Fernández, A., Manjarrez, J., García-Vázquez, U., D’Addario, M., & Sunny, A. 823 (2018). Present and future ecological niche modeling of garter snake species from the 824 Trans-Mexican Volcanic Belt. PeerJ, 6, e4618. https://doi.org/10.7717/peerj.4618 825 Goodman, K. R., Welter, S. C., & Roderick, G. K. (2012). Genetic divergence is decoupled 826 from ecological diversification in the Hawaiian Nesoydne planthoppers. Evolution, 827 66(9), 2798–2814. 828 Goslee, S., & Urban, D. (2020). ecodist: Dissimilarity-Based Functions for Ecological Analysis 829 (2.0.5) [Computer software]. https://CRAN.R-project.org/package=ecodist 830 He, K., Gutiérrez, E. E., Heming, N. M., Koepfli, K.-P., Wan, T., He, S., Jin, W., Liu, S.-Y., & 831 Jiang, X.-L. (2019). Cryptic phylogeographic history sheds light on the generation of 832 species diversity in sky-island mountains. Journal of Biogeography, 46(10), 2232– 833 2247. https://doi.org/10.1111/jbi.13664 834 Hubbell, S. (2001). The Unified Neutral Theory of Biodiversity and Biogeography. Princeton 835 University Press; JSTOR. https://doi.org/10.2307/j.ctt7rj8w 836 Huson, D. H., Beier, S., Flade, I., Górska, A., El-Hadidi, M., Mitra, S., Ruscheweyh, H.-J., & 837 Tappu, R. (2016). MEGAN Community Edition—Interactive Exploration and Analysis Page 33 of 45 Molecular Ecology For Review Only 34 838 of Large-Scale Microbiome Sequencing Data. PLOS Computational Biology, 12(6), 839 e1004957. https://doi.org/10.1371/journal.pcbi.1004957 840 Janzen, D. H. (1967). Why Mountain Passes are Higher in the Tropics Author ( s ): Daniel H. 841 Janzen Source: The American Naturalist , Vol. 101 , No. 919 ( May—Jun ., 1967 ), pp. 842 233-249 Published by: The University of Chicago Press for The American Society of 843 Naturalist. The American Naturalist, 101(919), 233–249. 844 Ji, Y., Ashton, L., Pedley, S. M., Edwards, D. P., Tang, Y., Nakamura, A., Kitching, R., 845 Dolman, P. M., Woodcock, P., Edwards, F. A., Larsen, T. H., Hsu, W. W., Benedick, 846 S., Hamer, K. C., Wilcove, D. S., Bruce, C., Wang, X., Levi, T., Lott, M., … Yu, D. W. 847 (2013). Reliable, verifiable and efficient monitoring of biodiversity via metabarcoding. 848 Ecology Letters, 16(10), 1245–1257. https://doi.org/10.1111/ele.12162 849 Kimura, M. (1968). Evolutionary rate at the molecular level. Nature, 217(5129), 624-626. 850 Kimura, M. (1969). The rate of molecular evolution considered from the standpoint of 851 population genetics. Proceedings of the National Academy of Sciences, 63(4), 1181852 1188. 853 Knowles, L. L. (2001). Did the Pleistocene glaciations promote divergence? Tests of explicit 854 refugial models in montane grasshopprers. Molecular Ecology, 10(3), 691–701. 855 https://doi.org/10.1046/j.1365-294x.2001.01206.x 856 Lee-Yaw, J. A., Davidson, A., McRae, B. H., & Green, D. M. (2009). Do landscape processes 857 predict phylogeographic patterns in the wood frog? Molecular Ecology, 18(9), 1863– 858 1874. https://doi.org/10.1111/j.1365-294X.2009.04152.x 859 Mastretta-Yanes, A., Moreno-Letelier, A., Piñero, D., Jorgensen, T. H., & Emerson, B. C. 860 (2015). Biodiversity in the Mexican highlands and the interaction of geology, 861 geography and climate within the Trans-Mexican Volcanic Belt.c, 42(9), 1586–1600. 862 https://doi.org/10.1111/jbi.12546 863 Mastretta-Yanes, A., Xue, A. T., Moreno-Letelier, A., Jorgensen, T. H., Alvarez, N., Pinero, D., 864 & Emerson, B. C. (2018). Long-term insitu persistence of biodiversity in tropical sky 865 islands revealed by landscape genomics. Molecular Ecology, 27(2), 432–448. Page 34 of 45Molecular Ecology For Review Only 35 866 https://doi.org/10.1111/mec.14461 867 McCormack, J. E., Huang, H., Knowles, L. L., Gillespie, R., & Clague, D. (2009). Sky islands. 868 In Encyclopedia of islands (pp. 841–843). 869 McGill, B. J. (2010). Towards a unification of unified theories of biodiversity. Ecology letters, 870 13(5), 627-642. https://doi.org/10.1111/j.1461-0248.2010.01449.x 871 McRae, B. H. (2006). Isolation By Resistance. Evolution, 60(8), 1551–1561. 872 https://doi.org/10.1554/05-321.1 873 McRae, B. H., Dickson, B. G., Keitt, T. H., & Shah, V. B. (2008). Using circuit theory to model 874 connectivity in ecology, evolution, and conservation. Ecology, 89(10), 2712–2724. 875 https://doi.org/10.1890/07-1861.1 876 McRae, B. H., Shah, V. B., & Mohapatra, T. K. (2013). Circuitscape 4 User Guides (4.0) 877 [Computer software]. The Nature Conservancy. http://www.circuitscape.org 878 Múrria, C., Bonada, N., Vellend, M., Zamora‐Muñoz, C., Alba‐Tercedor, J., Sainz‐Cantero, C. 879 E., ... & Derka, T. (2017). Local environment rather than past climate determines 880 community composition of mountain stream macroinvertebrates across Europe. 881 Molecular ecology, 26(21), 6085-6099. https://doi.org/10.1111/mec.14346 882 Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. 883 R., O’Hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. H., Szoecs, E., & 884 Wagner, H. (2019). vegan: Community Ecology Package (2.5-6) [Computer software]. 885 https://CRAN.R-project.org/package=vegan 886 Papadopoulou, A., Anastasiou, I., Spagopoulou, F., Stalimerou, M., Terzopoulou, S., Legakis, 887 A., & Vogler, A. P. (2011). Testing the species–genetic diversity correlation in the 888 Aegean archipelago: toward a haplotype-based macroecology?. The American 889 Naturalist, 178(2), 241-255. 890 Polato, N. R., Gill, B. A., Shah, A. A., Gray, M. M., Casner, K. L., Barthelet, A., Messer, P. W., 891 Simmons, M. P., Guayasamin, J. M., Encalada, A. C., Kondratieff, B. C., Flecker, A. S., 892 Thomas, S. A., Ghalambor, C. K., Poff, N. L., Funk, W. C., & Zamudio, K. R. (2018). 893 Narrow thermal tolerance and low dispersal drive higher speciation in tropical Page 35 of 45 Molecular Ecology For Review Only 36 894 mountains. Proceedings of the National Academy of Sciences, 115(49), 12471–12476. 895 https://doi.org/10.1073/pnas.1809326115 896 Pons, J., Barraclough, T. G., Gomez-Zurita, J., Cardoso, A., Duran, D. P., Hazell, S., Kamoun, 897 S., Sumlin, W. D., & Vogler, A. P. (2006). Sequence-Based Species Delimitation for 898 the DNA Taxonomy of Undescribed Insects. Systematic Biology, 55(4), 595–609. 899 https://doi.org/10.1080/10635150600852011 900 Rahbek, C., Borregaard, M. K., Colwell, R. K., Dalsgaard, B., Holt, B. G., Morueta-Holme, N., 901 Nogues-Bravo, D., Whittaker, R. J., & Fjeldså, J. (2019a). Humboldt’s enigma: What 902 causes global patterns of mountain biodiversity? Science, 365(6458), 1108–1113. 903 https://doi.org/10.1126/science.aax0149 904 Rahbek, C., Borregaard, M. K., Antonelli, A., Colwell, R. K., Holt, B. G., Nogues-Bravo, D., 905 Rasmussen, C. M. Ø., Richardson, K., Rosing, M. T., Whittaker, R. J., & Fjeldså, J. 906 (2019b). Building mountain biodiversity: Geological and evolutionary processes. 907 Science, 365(6458), 1114–1119. https://doi.org/10.1126/science.aax0151 908 Rambaut, A. (2012). FigTree v.1.4.2 Available http://tree.bio.ed.ac.uk/software/figtree/ 909 Rosindell, J., Hubbell, S. P., & Etienne, R. S. (2011). The unified neutral theory of biodiversity 910 and biogeography at age ten. Trends in ecology & evolution, 26(7), 340-348. 911 https://doi.org/10.1016/j.tree.2011.03.024 912 Rzendowski, J. (2006). Bosque de coníferas. Vegetación de México, 295–327. 913 Salces‐Castellano, A., Patiño, J., Alvarez, N., Andújar, C., Arribas, P., Braojos‐Ruiz, J. J., ... & 914 Manolopoulou, I. (2020). Climate drives community‐wide divergence within species 915 over a limited spatial scale: evidence from an oceanic island. Ecology Letters, 23(2), 916 305-315. https://doi.org/10.1111/ele.13433 917 Sheldon, K. S., Huey, R. B., Kaspari, M., & Sanders, N. J. (2018). Fifty Years of Mountain 918 Passes: A Perspective on Dan Janzen’s Classic Article. The American Naturalist, 919 191(5), 553–565. https://doi.org/10.1086/697046 920 Shokralla, S., Porter, T. M., Gibson, J. F., Dobosz, R., Janzen, D. H., Hallwachs, W., Golding, 921 G. B., & Hajibabaei, M. (2015). Massively parallel multiplex DNA sequencing for Page 36 of 45Molecular Ecology For Review Only 37 922 specimen identification using an Illumina MiSeq platform. Scientific Reports, 5(1), 923 9687. https://doi.org/10.1038/srep09687 924 Staton, E. (2019). Sestaton/Pairfq [Perl]. https://github.com/sestaton/Pairfq (Original work 925 published 2013) 926 Sunny, A., González-Fernández, A., & D’Addario, M. (2017). Potential distribution of the 927 endemic imbricate alligator lizard (Barisia imbricata imbricata) in highlands of central 928 Mexico. Amphibia-Reptilia, 38(2), 225–231. https://doi.org/10.1163/15685381929 00003092 930 Thompson, R., & Townsend, C. (2006). A truce with neutral theory: Local deterministic factors, 931 species traits and dispersal limitation together determine patterns of diversity in stream 932 invertebrates. Journal of Animal Ecology, 75(2), 476–484. 933 https://doi.org/10.1111/j.1365-2656.2006.01068.x 934 Uscanga A, López H, Piñero D, Emerson BC, Mastretta-Yanes A. (2021) Evaluating species 935 origins within tropical sky-islands arthropod communities. Journal of Biogeography, 936 00:1–12. https://doi.org/10.1111/jbi.14144 937 Vellend, M. (2010). 2010 Community Ecology. 85(2), 183–206. https://doi.org/10.1086/652373 938 Vellend, M., & Geber, M. A. (2005). Connections between species diversity and genetic 939 diversity. Ecology letters, 8(7), 767-781. 940 Wickham, H., Chang, W., Henry, L., Pedersen, T. L., Takahashi, K., Wilke, C., Woo, K., 941 Yutani, H., Dunnington, D., & RStudio. (2020). ggplot2: Create Elegant Data 942 Visualisations Using the Grammar of Graphics (3.3.2) [Computer software]. 943 https://CRAN.R-project.org/package=ggplot2 944 Yu, D. W., Ji, Y., Emerson, B. C., Wang, X., Ye, C., Yang, C., & Ding, Z. (2012). Biodiversity 945 soup: Metabarcoding of arthropods for rapid biodiversity assessment and 946 biomonitoring. Methods in Ecology and Evolution, 3(4), 613–623. 947 https://doi.org/10.1111/j.2041-210X.2012.00198.x 948 949 Page 37 of 45 Molecular Ecology For Review Only 38 950 Figure 1. Study site and sampling design. a) Location of study area within Mexico. b) 951 Nevado de Toluca volcano natural protected area over a topographic map. c) At each 952 sampling point (red circle) we sampled 3 sampling blocks of 20 x 15 m using 20 pitfall 953 traps equidistantly distributed. Sampling blocks were separated by 50 m (within sampling 954 points) to 19 km (between sampling points). Sampling points were distributed in four 955 sites to the West (ABB: Agua Bendita and ASB: San Bartolo) and East (SJH: San Juan 956 de las Huertas, and TCL: Tlacotepec) hillsides of Nevado de Toluca. Vegetation types 957 come from González-Fernández et al., (2019), see text for details. 958 959 Figure 2. ASVs diversity at Nevado de Toluca Abies forests. a) Phylogenetic tree 960 constructed by taxonomic composition of the arthropods estimated with the lowest 961 common ancestor (LCA) algorithm on MEGAN. b) Number of ASVs recovered at 962 haplotype level for each order. Numbers above the bar show putative species by GMYC 963 and values inside of brackets are lineages at 3%. 964 965 Figure 3. Arthropods richness at different clustering levels by sampling sites. AAB 966 = Agua Bendita, ASB = San Bartolo, SJH = San Juan de las Huertas, TLC = Santiago 967 Tlacotepec. Same letters indicate that there are no significative differences among those 968 sites p<0.05. Significance codes *:p<0.05, **:p<0.01, ***:p<0.001, ns: non-significant. 969 970 Figure 4. Non-Metric Multidimensional scaling (NMDS) ordinations of community 971 similarity (Simpson index, βsim) at the haplotype, 3% and 5% lineages for each of the 972 eight taxonomic orders studied. The asterisk represents significant differences of 973 community structure among sites estimated with ANOSIM at *:p<0.05, **:p<0.01 974 ***:p<0.001 and ns: non-significant values. Page 38 of 45Molecular Ecology For Review Only 39 975 976 Figure 5. Distance decay of community similarity at Nevado de Toluca for 977 Collembola (a,b) and Diptera (c,d) at multiple levels of genetic similarity. Decay of 978 similarity is shown against flat effective distance (‘flat’ surface; a,c) and against effective 979 distance using the resistance surface “Altitude 3,000” for Collembola (b) and ‘Mean980 elevation peak (symmetrical)’ for Diptera (c), which were the ones showing the highest 981 explanatory power (significance levels in Table S3). Dots represent similarity between 982 pairs of sampling sites and lines are the fitted model for each lineage. Altitude rasters as 983 in Figure S5. Significance levels and r2 are in Table S3 and S4. Notice that effective 984 distances are resistances and hence are not comparable to km. 985 986 Figure 6. Distance decay of community similarity at fine geographic distances within 987 Nevado de Toluca for Collembola (a,b) and Diptera (c,d). Analysis at <5 km (a,c) 988 using the Eastern subset of sampling sites (SJH and TLC) and at <2 km (b,d) using the 989 Western subset (AAB and ASB). Decay of similarity is shown against geographic 990 distance (‘flat’ surface). Dots represent pairs of sampling sites for each clustering level 991 and lines are the fitted model. Flat rasters are in Figure S5. Significance levels and r2 are 992 in Table S5. Page 39 of 45 Molecular Ecology For Review Only Figure 1. Study site and sampling design. a) Location of study area within Mexico. b) Nevado de Toluca volcano natural protected area over a topographic map. c) At each sampling point (red circle) we sampled 3 sampling blocks of 20 x 15 m using 20 pitfall traps equidistantly distributed. Sampling blocks were separated by 50 m (within sampling points) to 19 km (between sampling points). Sampling points were distributed in four sites to the West (ABB: Agua Bendita and ASB: San Bartolo) and East (SJH: San Juan de las Huertas, and TCL: Tlacotepec) hillsides of Nevado de Toluca. Vegetation types come from GonzálezFernández et al. (2019), see text for details. 300x270mm (96 x 96 DPI) Page 40 of 45Molecular Ecology