Neutral‐based processes overrule niche‐based processes in shaping tropical montane orchid communities across spatial scales

Tropical montane forests (TMF) are characterised by high endemism, species richness and turnover across elevations. A key question is how niche‐based processes, via adaptation of species to local environmental conditions, and neutral‐based processes from dispersal limitation shape community composition at different spatial scales across human‐modified landscapes. We expect that (1) communities are highly distinct even within the same habitat type and (2) niche‐based processes play the main role in compositional turnover. To address these expectations, we investigated the compositional turnover of orchids, one of the most prominent floristic elements of TMFs. We sampled orchids in 332 plots spanning over 270 km in the eastern Colombian Andes. Plots ranged between elevations of 1160–3750 m. We used two different spatial extents (whole gradient and two elevational bands), two grains of analysis for the first expectation (regional and local) and two spatial grains for the second expectation (broad and fine scales based on Moran's Eigenvector Maps [MEMs]). We found 331 orchid species in 171 (51%) plots. We found a strong pattern of high compositional turnover across scales (>72% of total beta diversity is given by species turnover), with 87.5% of the total species pool occurring in fewer than five plots, supporting our first expectation. Contrary to our second expectation, we found that community composition is mainly driven by geographical distance, while the relative influence of elevation, environmental variables and their combined fractions were negligible across habitats and spatial scales, rejecting niche‐based processes. Synthesis. High compositional turnover, even across habitats with the highest degree of human intervention, suggests that both forest‐dwelling and open‐habitat species do not easily disperse across habitats. Species dispersal is the major force of orchid community turnover and might be strongly dependent upon macroevolutionary processes and species life‐history traits over multiple scales. Dispersal limitation also draws attention to the importance of preserving habitat connectivity to halt biodiversity losses.


| INTRODUC TI ON
Mountains harbour a global hyperdiversity of species and endemism (Enquist et al., 2019;Küper et al., 2004;Peters et al., 2019). The complex geological history of mountains boosted the evolution of many plant groups by promoting speciation through, for instance, the provision of new habitats, or ecological and dispersal barriers to other species and other populations of the same species (Luebert & Weigend, 2014;Pérez-Escobar, Gottschling, et al., 2017). The replacement of species between sites-compositional turnover-is the prevailing pattern across African (Peters et al., 2019), European (Fontana et al., 2020) and Andean mountainous regions (Herrera-Pérez et al., 2019;Muñoz Mazón et al., 2021;Novillo & Ojeda, 2021).
However, the ecological processes governing high compositional turnover typical of mountainous ecosystems remain understudied at multiple spatial scales.
Patterns in compositional turnover can be explained by nichebased processes, that is the adaptation of species to local environmental conditions, and neutral-based processes such as dispersal limitation (Hubbell, 2001;Leibold et al., 2004;Tuomisto et al., 2003).
Niche-based processes act when communities respond to environmental heterogeneity due to the availability of niches and coexistence with other species (Chase & Leibold, 2003). In mountainous regions, environmental variability due to habitat patchiness, drastic changes in precipitation, and topography increase compositional turnover (Arellano et al., 2016;Du et al., 2021;Trujillo et al., 2019).
In contrast, neutral-based processes become important as species and their populations fluctuate at random due to demographic stochasticity and spatially limited dispersal (Hubbell, 2001;Rosindell et al., 2011). Species with limited dispersion can sometimes propagate and successfully colonise within narrow geographical ranges regardless of environmental heterogeneity. Strong dispersal limitation in mountainous regions is expected due to abiotic factors, such as high topographic variability and landscape fragmentation, especially in communities dominated by rare species (Arellano et al., 2017;Myers et al., 2013;Willig et al., 2011).
The pattern of high compositional turnover can emerge from niche-based and dispersal limitation processes alike. If ecological specialists respond to very specific environmental conditions, to dispersal limitation because species cannot disperse long distances, or to both, then it is reasonable to expect communities with a large pool of these specialists to exhibit high compositional turnover. Studies of woody plants across elevation gradients have shown that nichebased processes are important in temperate forests (e.g. Myers et al., 2013) and that both niche-and neutral-based processes drive community structure in tropical forests, probably due to recruitment limitation or strong environmental gradients (Arellano et al., 2017;Myers et al., 2013;Tello et al., 2015;Trujillo et al., 2019), while vascular epiphytes seem to respond to near-neutral processes (Janzen et al., 2020).
These processes can be further influenced by the spatial scale of study and anthropogenic activities, especially habitat loss (Arroyo-Rodríguez et al., 2013, 2017. The spatial scale of study plays an important role in how communities respond to environmental changes. Geographic distance shapes regional composition patterns across communities (neutral-based processes). Equally, environmental differences operate mostly at local scales (Karp et al., 2012) because longer and steeper environmental gradients allow a more marked sorting of species by their different preferences (Trujillo et al., 2019), even across human-modified landscapes (Du et al., 2021;Sreekar et al., 2020).
Human modification of natural landscapes also alters the processes governing community composition (Haddad et al., 2015;Venter et al., 2016). The novel conditions generated by human activities can trigger homogenization in species pools, the process by which a small set of organisms can cope and even flourish after the alteration of natural habitats, that is, biotic homogenization (McKinney & Lockwood, 1999). For instance, human transformation of the landscape reduces structural connectivity causing harsh conditions that inhibit dispersal of forest-specialist or dispersal-limited species , while habitat-specialists increase spatial aggregation due to the reduction of environmental heterogeneity (Audino et al., 2017;Karp et al., 2012). A better understanding of these patterns is required to quantify the influence of neutral-and niche-based processes across spatial scales and human-modified landscapes. Empirical research is now dedicated to disentangling the contributions of niche-based and neutral-based processes using pattern-to-process approaches across spatial scales, while controlling for spurious environment-space correlations (Clappe et al., 2018).
We quantified compositional turnover and the roles of nichebased and neutral-based (i.e. dispersal limitation) processes across different spatial scales in a large tropical montane region.
In particular, we quantified: (i) patterns of compositional turnover across regional and local scales, and elevations; and (ii) the independent and combined roles of environmental predictors (natural and anthropogenic) and distance (i.e. geographical and elevational distance) in shaping compositional turnover. We focus on the Colombian Andes, an area of highly diverse topographical and hydrological conditions (Bürgl, 1967) that cause large shifts limitation also draws attention to the importance of preserving habitat connectivity to halt biodiversity losses.

K E Y W O R D S
determinants of plant community diversity and structure, ecological process, neotropics, northern Andes, Orchidaceae, spatial scales in biological communities over short distances (Hazzi et al., 2018;Malizia et al., 2020;Rahbek et al., 2019). We use orchids as a model group because the tropical Andes is the centre of global orchid diversity (Pérez-Escobar et al., 2022) with one of the highest species richness worldwide, with a wide variation in species' range sizes from widespread to highly restricted (Calderón-Saenz, 2007;Jost, 2004;Moreno et al., 2020). The majority species use wind as their main dispersal strategy (anemochory syndrome), while interspecific competition, and the effects of herbivory and pathogens are expected to be negligible (Dressler, 2005;Zotz, 2016). We hypothesised that local orchid communities are highly distinct from one another, even within the same habitat type (high turnover). We expect that compositional turnover would be driven primarily by environmental variables, regardless of the spatial scale of study, because many orchid species have been documented to respond to environmental conditions (Calderón-Saenz, 2007;Mondragón et al., 2015;Zuleta et al., 2016).

| Study area
The study was located in the departments of Cundinamarca, Boyacá, Meta and Santander, spanning both the east and west flanks of the eastern Andes cordillera in Colombia (Figure 1). Between January 2018 and November 2019, we sampled 332 plots across 270 km of distance, covering a wide elevational range, 1163-3763 m (4-28°C and 879-3817 mm per year). The study area included three habitat types: (i) forests, (ii) pasture for cattle and (iii) paramo shrublands and grasslands. All habitats have undergone historical humanmediated disturbance given the easy accessibility by roads or footpaths, and information from local field assistants at each sample location (pers. comm. landowners). All sampled forests belonged to the network of protected areas in Colombia (Sistema Nacional de areas Protegidas-SINAP).
F I G U R E 1 Distribution of sampling plots and characteristics of the study area in the eastern cordillera of the Colombian Andes. The map shows the slopes across sites (Hillshade at 45% slope and 270% azimuth), elevation (c and d) jointly with forest cover (green) and the absence of forest cover (grey). Inset maps and photographs as follows:

| Sampling design
The sampling design followed the natural observed tree line along the elevational gradient ( Figure S2). In the Andes, the treeline is considered the upper-limit boundary of the establishment of arboreal vegetation between upper Andean montane forest and low-stature vegetation (Bader et al., 2007). Below the observed treeline (at elevations of 1163-2900 m), we paired forests with nearby pastures, and above the treeline, natural paramo shrublands, grasslands and high elevation elfin forests were paired with high elevational pasturelands ( Figure 1; henceforth paramo >2900 m).
Below the tree line, we randomly placed between 1 and 18 sam- Each sampling plot consisted of a 10 m × 30 m plot. We recorded orchids up to 2 m above-ground within each plot on the ground, standing trees, tree trunks, fallen tree trunks and branches, vines, lianas, leaves on standing trees or herbaceous plants, palm trees, tree ferns and cycads across the understorey. Canopy orchid species might have been sampled from fallen branches, but since we used the same sampling across sites, all plots would have the same probability of capturing canopy orchids. Only adult individuals were recorded. Adult plants were all large individuals, that is, within the observed variation in the local population, with evidence of developed floral or reproductive structure, large ramets, various stems per stand and prominent root systems. Species with no flowers were taken to nurseries in the vicinity of the sampled plots, and later visited for identification. Identification to species or morphospecies was conducted following specialised literature and consultancy with local experts at the Herbarium VALLE (Palmira, Colombia).

| Environmental predictors
Habitat structure was composed of tree diameter at breast height (average DBH per plot) and tree density (trees per plot) reduced to one single predictor using principal component analysis (PCA), hereafter "Forest structure" (first PCA axis accounted for 61% of the variation, SD = 1.1043, supplemental material). Cloud cover We quantified the percent of forest cover in the local landscape as the effects of habitat loss (Fahrig, 2013). The percentage of forest cover is a proxy of habitat loss that calculates the amount of habitat surrounding each of our sampling plots. Forest cover is a metric of landscape composition that quantifies the habitat available in the landscape as a proxy of resource availability (Fahrig, 2013) but cannot inform how the habitat in the landscape is connected (structural connectivity). Thus, the operational concept we use in our study was habitat availability in the landscape as an environmental driver of community composition (e.g. Banks- Leite et al., 2014;Püttker et al., 2015). Forest cover was measured from the 30-m resolution global change forest map (Hansen et al., 2013), converted into a binary forest/non-forest cover map (non-forest <= 50, and >50% "forest" in 2018; Supporting Information). We selected a 1000 m radius to collate forest cover to match the climatic predictors, as we did not evaluate the scaleof-effects (Jackson & Fahrig, 2015).

| Distance predictors
We defined two independent distance predictors based on (i) geographical coordinates (XY coordinates) and (ii) elevation (m) measured in the field and later standardised with the 12.5-m ALosPalsar's Radiometric terrain model (Tadono et al., 2014). All spatial predictors were generated using Moran's Eigenvector Maps (hereafter MEMs; Dray et al., 2006). These spatial predictors are orthogonal vectors maximising the spatial autocorrelation (calculated by Moran's coefficient), acting as spatially explicit multiscale predictors that combine fine and broad scales (Dray et al., 2012). First, geographical distance MEMs were obtained using the geographical coordinates projected to WGS1984 and Gabriel's neighbour network and W standardisation (sums over all links of the number of plots). Second, we created MEMs by including the latitude coordinate of each point to the corresponding elevation to generate a XY coordinate-like vector. Then, we calculated elevation MEMs using neighbourhood-by-distance method and globally standardisation (sums over all links to plots). Finally, we defined two spatial grains for the whole gradient, MEMs were divided into two groups corresponding to different spatial scales ranging from the broadest (MEM with highest associated eigenvalues) to the finest spatial scales (MEM with lowest associated eigenvalues; Menegotto et al., 2019). This approach suits our nested design with strong positive spatial autocorrelation.
The intrinsic spatial autocorrelation between distance and environmental variables constitutes an important challenge for ecological studies. However, we retained distance (MEMs for geographical distance and elevation) as a proxy for isolation because evidence shows that orchid seeds tend to land <5 m from their mother (Acevedo et al., 2020;Kindlmann et al., 2014;Tremblay, 1997), a pattern resulting from the rapid in-situ speciation process in the Andes . Furthermore, ecological research shows that short-distance dispersal is common in epiphytic communities, with epiphytism the dominant habitat syndrome in Andean orchids (Zotz, 2016), dwelling in pasturelands and protected areas (Cascante-Marín et al., 2009;Einzmann et al., 2021;Einzmann & Zotz, 2017a, 2017b).

| Statistical analyses
We ran analyses based on our two expectations: (i) to quantify the patterns of compositional turnover across regional and local scales, and elevations; and (ii) to quantify the roles of niche-based (natural and anthropogenic drivers) and neutral-based processes (i.e. geographical and elevational distance) in shaping compositional turnover. We used two different spatial extents (whole gradient and two elevational bands), two grains of analysis for the first expectation (regional and local) and two spatial grains for the second expectation (broad and fine scales based on MEMs).

| Compositional turnover
We tested for the patterns of compositional turnover across regional and local scales, and elevations using species occurrence by site matrices and calculate the Simpson dissimilarity index referring to the spatial turnover component (species replacement) of the Sørensen dissimilarity index (Baselga, 2010). Compositional turnover ranges from zero (no turnover or all species are shared between sites), and one when no species are shared.
To explore the general pattern in compositional turnover, we used a generalised additive mixed model (Wood, 2017). We fitted the model with the average turnover of each sampling plot as the response variable, elevation as the fixed-effect and sampling point as a random effect, with a binomial family. To account for spatial correlations, we added a spatial covariance structure using the geographical coordinate of each sampling plot.
We further quantified compositional turnover within and between habitats at regional and local scales. At the regional scale, we quantified the pairwise dissimilarity of each habitat type as a single metric of compositional turnover of each habitat type (betapart::beta.multi). For this analysis, we averaged the spatial turnover of each community at each sampling plot within the same forest fragment (n = 23), within the same pastureland (n = 22), and within the same paramo (n = 8). The mean pairwise distance is a robust estimate of species turnover regardless of the number of samples (Pfeifer et al., 2017;Olivares & Kessler, 2020). We fitted generalised least-squared models with an exponential spatial structure class with the geographical coordinates of the centroid of each habitat because we detected spatial autocorrelation between communities and predictors (Moran I test = 0.529, p-value = 0.001; At the local scale, we established the sampled plots within each sampling habitat as units of analysis (n = 171 out 332 sampling plots) and quantified compositional turnover within habitats (i.e. turnover within forests, within pastures and within paramos) and among habitat types (forest vs. pastures, forest vs. paramo, and pasture vs. paramo). We computed the mean value of all pairwise dissimilarities of each sampling point: within-forests, within-pastures, forest versus pasture, within-paramo, and paramo versus forest. We fitted generalised linear-mixed models including compositional turnover at each sampling plot as the response variable, habitat type as the fixed effect, sampling plot and habitat type as random effects, and a spatial exponential covariance structure using the geographical coordinates of each sampling plot to account for spatial autocorrelations (formula in glmmTMB notation in Supporting Information; Brooks, Kristensen, et al., 2017;Brooks, Mittermeier, et al., 2017).
Our turnover data was one-skewed with few zero values (85.7% of species appear in less than five sampling plots, Table S2), so we subtracted 1 and multiplied by −1 each turnover value to make the data zero-skewed and fit zero-inflated models while preserving the structure of the data. The zero-inflated model was fitted with beta distribution family (suitable for proportional data between bounded limits, zero to one; Cribari-Neto & Zeileis, 2010), the dispersion error log link, and geographical coordinates as spatial exponential covariance (zero-inflated formula notation in Supporting Information). To further control for the effect of elevation across local landscapes (local scale) using the hump-shaped pattern seen in epiphytes and orchid richness and total abundance (E. Parra-Sanchez, unpubli. data; Figure S2), we divided our dataset into a low-mid elevational band where orchid richness peaks ('band-1'; 1150-2500 m) and midhigh elevational band where orchid richness decreases ('band-2'; 2501-3500 m). This design allows us to decompose the effects of elevation and spatial scales in compositional turnover.

| Quantifying the roles of niche-based and neutral-based processes
To determine the relative influence of niche-based (environment) and neutral-based (distance) processes, we used variation partitioning jointly with the Moran spectral randomization (MSR) method. MSR accounts for spurious spatial correlations (Borcard et al., 1992;Clappe et al., 2018;Peres-Neto et al., 2006;Wagner & Dray, 2015).
Variation partitioning is a widely implemented statistical tool used to disentangle the roles of niche-based and neutral-based processes (via dispersal limitation) by quantifying the relative importance of variation in compositional turnover explained by the environmental and spatial components and their shared effect (Borcard et al., 1992;Smith & Lundholm, 2010 We used constrained null models from MSR to account for spatial autocorrelations in variation partitioning. Given the ubiquity of spatial autocorrelation in environmental data, spurious spatial autocorrelations can bias direct gradient analyses and inflate the variation explained by the environmental component, underestimating the role of spatial structure (Clappe et al., 2018;Wagner & Dray, 2015). MSR allows for the separation of the contributions of spatially structured communities via environment, by confronting the observed result with what is expected under neutral dynamics due to spurious relationships between autocorrelated species and environmental dynamics (Clappe et al., 2018;Wagner & Dray, 2015).
We ran three sets of variation partitioning analyses: (i) whole gradient, (ii) by habitat type and (iii) partitioned by elevational bands (low-mid band-1 and mid-high band-2). Each analysis was performed for geographical distance (GD + environment), elevation (Ele + environment) and both distances together (GD + Ele; except for the habitat analysis). Each variation partitioning model was fitted with the principal component analysis of the orchid communities as response variable, the generated spatial predictors (i.e. MEM's of geographical and elevational distance) and the environmental predictors. Prior to analyses, the species by site matrix was Hellinger transformed (Legendre & Gallagher, 2001). This standardised matrix was analysed via centred principal component analysis using adjusted
Eleven new orchid species to science were found in the forest habitat, with published descriptions forthcoming (E. Parra-Sanchez, unpubli. data), while no species were found in high-elevation pastures.
Compositional turnover was also high across elevational bands.

| Quantifying the roles of niche-based and neutral-based processes
Geographical distance (proxy of niche-based process) was the single best predictor that consistently explained the variation of compositional turnover across scales, elevation and habitat types, albeit the residual (unexplained) fraction was the largest fraction across models (p-value < 0.05; Supporting Information 2). The principal component models of community composition were significant across analyses (geographical distance, adjusted R 2 = 0.18, p-value < 0.01; elevation adjusted R 2 = 0.26, p-value < 0.001).
We found that the spatial fraction (neutral-based process) in the form of geographical distance accounted for the largest proportion of the variance explained in orchid communities against the environmental fraction (GD = 10.9%-15.5% vs. Ele = 0.3%-1.1%) and the shared fraction (GD + Environment; GD = 10.4-14.9 vs. GD + Ele = 0.06-0.90; Figure 5). Elevation was more important than the environmental fraction in explaining compositional turnover only when analysing the whole gradient, albeit with very low variance explained (2.2% vs. 1.4% environment). Likewise, the spatial fraction (neutral-based process) was important across the two spatial grains of analyses, at the fine-scale (variance explained 2.4 %; F I G U R E 2 Elevation pattern of orchid compositional turnover. Tick marks on the x-axis represent sampled plots along the elevation gradient and the Y-axis represent the smooth function of the centred values of compositional turnover. The grey area shows 95% confidence intervals around the prediction.

F I G U R E 3
Compositional turnover across forest, pasture and paramos at regional (a) and local scale (b). The shaded area indicates the kernel density of the turnover range of each habitat type. The violin plots show the median turnover value, the 25th and 75th percentile range for habitat type. The circular dots represent the outliers in observed records of compositional turnover at each sampled plot.
The environmental fraction explained more variance for band 1 and band 2 than elevation per se, although explaining a very low proportion of the variance (band 1, 1.9 vs. 0.03 elevation; band 2, 2.6 vs.
1.9 elevation). The shared influence of the GD + Environment fraction was negligible in explaining community structure across models and elevation bands (0.11%-1.39%, Supporting Information). Across habitat types, communities were also more responsive to GD as the major process involved shaping forest habitats (15.6% vs. 0.03% environment) and pasture habitats (5.82 % vs. 0.01% environment; Supporting Information 2).

| DISCUSS ION
We found a strong pattern of high compositional turnover driven by dispersal limitation (neutral-based process) shaping communities along elevations, habitats and spatial scales. We provide evidence of highly distinct communities even within the same habitat, where communities are mainly composed of a large pool of rare species that might be more responsive to dispersal limitation and much less to environmental (natural and anthropogenic) conditions across spatial scales.

| Compositional turnover
Our expectation of high turnover was validated because a large proportion of the sampled species (85.7%) were found in fewer than five plots even within the same habitat. Although several studies have shown the same pattern across different taxa, our study captures it at the finest scale (within the same habitat; Malizia et al., 2020;Muñoz Mazón et al., 2021;Tello et al., 2015). Wind-dispersed organisms are expected to attain larger geographical distribution ranges and be less sensitive to geographical distances (Qian & Guo, 2010;Tuomisto et al., 2003). However, our findings are consistent with previous studies in which local epiphyte communities have few highly abundant and several rare species (Benavides et al., 2011;Janzen et al., 2020;Laube et al., 2006). Communities show extremely local diversity-clusters with populations of some species clumped within the same forest (Acevedo et al., 2020;Kindlmann et al., 2014;Olaya-Arenas et al., 2011). Likewise, as reported for epiphytic vascular communities, orchids might show some degree of near-neutral structuring, potentially due to their high degree of ecological equivalence among species (Catchpole & Kirkpatrick, 2010;Janzen et al., 2020). Evidence also suggests that populations within forests have very low dispersal ranges (4.8 m on average) and their dispersion is asymmetrically driven by wind direction (Murren & Ellison, 1998), with low colonisation success (Acevedo et al., 2020;Kindlmann et al., 2014;Olaya-Arenas et al., 2011).
The main pattern on mountainous communities consists of high turnover at low elevations and lower turnover at high elevation across taxa (Du et al., 2021;Kraft et al., 2011;Peters et al., 2019).
Our results mirror this pattern, revealing decreasing compositional turnover as elevation increases, possibly due to lower availability of habitats (e.g. absence of trees in paramo shrublands and grasslands for epiphytic plants), evolutionary isolation, and lack of more stable eco-evolutionary adaptations to high elevation habitats such as UV-protection (Arellano et al., 2017;Trujillo et al., 2019).
There is evidence that dispersal limitation, not elevation, affects colonisation-extinction dynamics in a well-studied orchid species at F I G U R E 4 Local compositional turnover at (a) band 1 between 1163 and 2500 m; and (b) band 2 between 2500 and 3763 m within forest, pasture and paramos across all elevational bands. The shaded area indicates the kernel density of the turnover range of each habitat type. The violin plots show the median turnover value, the 25th and 75th percentile range for habitat type. The circular dots represent the outliers in observed records of compositional turnover at each sampled plot. mid elevations (Acevedo et al., 2020). However, the role of elevation on groups with large numbers of rare species would require more refined studies on its role in compositional turnover.
We found that compositional turnover among paired forest and pasture plots was very high (79%-100% in average of the total beta diversity is given by species turnover across scales), implying that forest species cannot tolerate pasture conditions. This pattern across scales and elevational bands shows that harsh conditions in pastures across human-modified landscapes, including lower number of phorophytes, higher temperature (vs forest), and lower humidity (vs forest), have strong effects on species composition (de la Rosa-Manzano et al., 2014;Einzmann & Zotz, 2017;Zotz & Bader, 2009). In contrast, empirical and theorical studies reveal that community convergence (i.e. the decrease in compositional turnover between communities, biotic homogenization) is associated with human-modification of landscapes, pointing to small subsets of species capable of surviving under harsh conditions (Brooks, Kristensen, et al., 2017;Brooks, Mittermeier, et al., 2017;Jamoneau et al., 2012;McKinney & Lockwood, 1999).
We found that pastures had one of the highest compositional turnovers across habitats (88%-100% of species turnover), showing that few species are shared within the habitat with the highest degree of human intervention. The high orchid turnover in pastures implies that community convergence is only possible if species can disperse and successfully colonise these less-suitable F I G U R E 5 Variation partitioning of the relative effects of elevation (Ele), environmental variables and geographical distance (GD) on orchid community composition. Community variation was partitioned using principal component analysis on Hellinger-transformed community data. (a) shows the relative loads of the environmental variables (cloud: percent of cloud cover, precipitation (mm), tree structure (PCA axis 1); elevation; forest cover (%)) across communities for the geographical distance model of the whole gradient (GD in the barplot). (b) depicts the relative loads of geographical distance MEMs predictors across communities for the geographical distance model of the whole gradient (GD in the barplot). (c) shows the variance explained of all models from the variation partitioning models after MSR corrections (adjusted coefficient of multiple determination). From top-down, barplots show models based on geographical distance for the whole gradient (GD), GD of band 1 (GD B1), GD of band 2 (GD B2), GD of broad spatial scales (GD broad refers to models using broad spatial scales MEMs), GD of fine spatial scales (GD fine refers to models using fine spatial scales MEMs); and then, distance models based on elevation, for the whole gradient (elevation), elevation of band 1 (elevation B1), elevation of band 2 (elevation B2), elevation of broad spatial scales (elevation broad refers to models using broad spatial scales MEMs) and elevation of fine spatial scales (elevation fine refers to models using fine spatial scales MEMs). Finally, barplots of geographical and elevational MEMs in the same models for the whole gradient (GD + Ele), across band 1 (GD + Ele B1) and band 2 (GD + Ele B2).

| The roles of niche-based and neutralbased processes
Our study supports neutral-based processes shaping mountainous community turnover across spatial scales, elevations, and habitat types with low to almost negligible effects of the environment and its shared proportion (Environment+Space). This aligns with research suggesting dispersal limitation as the dominating process in hyperdiverse forests with high levels of taxonomic rarity (Du et al., 2021;Muñoz Mazón et al., 2021;Myers et al., 2013), with the pattern potentially emerging when present-day tree species' dissimilarity responds to low landscape connectivity in the past Hurtado-M et al., 2021). For instance, 58% of endemic orchids in Ecuador are restricted to one forest type and 26% to two forest types (Endara et al., 2009), and 78% of orchid species from three genera in the Eastern cordillera of Colombia have been found only once or twice in its natural habitat (E. Parra-Sanchez, unpubli. data). Thus, dispersal-limitation processes shape our communities structured largely by uncommon species, consistent with many Neotropical orchid species that have narrow distributions (Jost, 2004), exhibit high mortality rates (Zuleta et al., 2016), or have only one successful migrant per generation (Tremblay et al., 1998).
We found several species of the same genus coexisting and new species of the most species-rich genus in the same local habitat (genera Lepanthes and Epidendrum, E. Parra-Sanchez, unpubli. data), signalling potential parapatric evolution (i.e. speciation when populations overlap geographically). Following the rapid origination of new montane ecosystems after to the fluid orogenetic history of Northern Andes (Pérez-Escobar et al., 2022), orchid lineages experienced rapid in-situ diversifications in these newly available ecosystems , specialising in micro-habitats (e.g. tree branches, boles). Such rapid diversification was also fuelled by the evolution of specific interactions with pollinators (e.g. Ackerman et al., 2023;Ramírez et al., 2011) and the abundance of mycorrhizal fungi necessary for germination and survival (McCormick et al., 2018;Romero-Salazar et al., 2022), thus enabling speciation under either parapatric or peripatric processes (Pérez-Escobar, Gottschling, et al., 2017). Since mycorrhizae are adapted to very specific substrate chemistry and conditions, as well as pollinators, the distribution of orchids might be strongly linked to those of their corresponding mycorrhiza and pollinators. It is tempting to hypothesize that a niche-based process might be acting at an extremely narrow spatial scale promoting niche partitioning (Mccormick & Jacquemyn, 2014), undetectable at the scale we studied here. Thus, our findings align with the idea of regional neutrality as a consequence of local adaptive niche evolution and niche conservatism (Linck et al., 2021).
Overall, the influence of geographical distance consistently was the best predictor across analysis, but the unexplained component dominated the observed variance. The unexplained fraction should not immediately translate into complete stochasticity. In variation partitioning analyses, the unexplained fraction varies with sampling (e.g. length and location along the gradient, or under sampling), and the spatial structure of unmeasured environmental variables.
Rapid assessments can overestimate compositional turnover because of 'sampling effects' and confound observed compositional turnover when quantifying the relative importance of local driving mechanisms (Tuomisto et al., 2012;Wagner & Dray, 2015). Also, our sampling largely included understory communities where conditions might act as an environmental 'prefilter' for shade-preferring species (vs canopy epiphytes). Nonetheless, we found eleven new orchid species to science (E. Parra-Sanchez, unpubli. data) and three other species never found before in their natural habitat (E. Parra-Sanchez, unpubli. data), and 49.6% of species sampled here have been recorded only a few times in nature (<= 10 records; E. Parra-Sanchez, unpubli. data). Our findings highlight that most species are spatially restricted suggesting that our sampling did not overly bias our results. Hence, the full spectrum of community turnover process might require the inclusion of more refined predictors related to biotic interactions at lower spatial resolution (Taylor et al., 2019), the inclusion of canopy species (Parra-Sanchez & Banks-Leite, 2020), and assessing the scale-of-effects of landscape predictors (Jackson & Fahrig, 2015).

| Management implications and conclusion
Our findings of the role of dispersal limitation in structuring orchid communities draws attention to the need to value all forest fragments in the landscape to preserve high levels of structural connectivity. For instance, private protected areas connected to the network of national parks, that are usually large unconnected protected areas, may improve the probability of an adequate number of propagules to disperse across landscapes (Acevedo et al., 2020;Arroyo-Rodríguez et al., 2020).
The roles of niche-and neutral-based processes are one of the most debated concepts in modern ecology because they can explain the prevailing forces shaping community assembly. Across the most hyperdiverse mountainous region worldwide, we showed that dispersal limitation is a key driver of compositional turnover across wind-dispersed communities. The mechanisms influencing species turnover might be controlled by direct and indirect dynamics in communities across spatial gradients. For example, species dispersal might be strongly dependent upon macroevolutionary processes and species life-history traits that structure communities over multiple scales.

AUTH O R CO NTR I B UTI O N S
Edicson Parra-Sanchez and David P. Edwards designed the conceptual framework of the study. Edicson Parra-Sanchez collected the data and led the analysis and interpretation of the data with support of David P. Edwards and Oscar A. Perez-Escobar. Edicson Parra-Sanchez wrote the initial drafts of the manuscript. Oscar A.
Perez-Escobar and David P. Edwards contributed critically to the drafts and gave final approval for publication. Foundation. This study was funded by the Natural Environment

ACK N O WLE D G E M ENTS
Research Council (grant no. NE/R017441/1). We express our gratitude to the very constructive comments that the anonymous reviewers made to the manuscript and to the diligent work that the handling editor put in place. Muchas gracias. This is publication #28 of the Biodiversity, Agriculture, and Conservation in Colombia (Biodiversidad, Agricultura, y Conservación en Colombia [BACC]) project.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors express no conflict of in interest.

PEER R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-revie w/10.1111/ 1365-2745.14140.

DATA AVA I L A B I L I T Y S TAT E M E N T
Community composition matrix and predictor data are available from  (species present in one plot), doubletons (species present in two plots), species present across 3-5 plots, 5-10 plots, or >10 sampling plots. Table S3. Beta diversity average (mean), maximum (Max), minimum (Min) and standard deviation (SD) across orchid communities within each habitat type at regional scale.