Seeing the wood despite the trees: Exploring human disturbance impact on plant diversity, community structure, and standing biomass in fragmented high Andean forests

Abstract High Andean forests harbor a remarkably high biodiversity and play a key role in providing vital ecosystem services for neighboring cities and settlements. However, they are among the most fragmented and threatened ecosystems in the neotropics. To preserve their unique biodiversity, a deeper understanding of the effects of anthropogenic perturbations on them is urgently needed. Here, we characterized the plant communities of high Andean forest remnants in the hinterland of Bogotá in 32 0.04 ha plots. We assessed the woody vegetation and sampled the understory and epiphytic cover. We gathered data on compositional and structural parameters and compiled a broad array of variables related to anthropogenic disturbance, ranging from local to landscape‐wide metrics. We also assessed phylogenetic diversity and functional diversity. We employed nonmetric multidimensional scaling (NMDS) to select meaningful variables in a first step of the analysis. Then, we performed partial redundancy analysis (pRDA) and generalized linear models (GLMs) in order to test how selected environmental and anthropogenic variables are affecting the composition, diversity, and aboveground biomass of these forests. Identified woody vegetation and understory layer communities were characterized by differences in elevation, temperature, and relative humidity, but were also related to different levels of human influence. We found that the increase of human‐related disturbance resulted in less phylogenetic diversity and in the phylogenetic clustering of the woody vegetation and in lower aboveground biomass (AGB) values. As to the understory, disturbance was associated with a higher diversity, jointly with a higher phylogenetic dispersion. The most relevant disturbance predictors identified here were as follows: edge effect, proximity of cattle, minimum fragment age, and median patch size. Interestingly, AGB was efficiently predicted by the proportion of late successional species. We therefore recommend the use of AGB and abundance of late successional species as indicators of human disturbance on high Andean forests.


| INTRODUC TI ON
High Andean tropical montane forests (herein bosques altoandinos) can be found between ca. 2,700 and 3,300 m in the Northern Andes, extending from Venezuela to Ecuador, with considerable levels of species diversity and endemism (Gentry & Ortiz, 1993;Girardin et al., 2014;Killeen et al., 2007;Still et al., 1999;Young, 1992).
Bosques altoandinos have been subjected to extensive anthropogenic transformation across their natural range. In Colombia, large portions of the forest cover were cleared during the past four centuries and turned into agricultural or residential areas, in order to satisfy the growing demand for resources of an increasing human population (Brown & Kappelle, 2001;Cavelier et al., 2001;Etter et al., 2008;Heath & Binswanger, 1996;Sánchez-Cuervo et al., 2012;Wassenaar et al., 2007). Such a reduction of forest cover can not only lead to loss of biodiversity but also to a lower structural integrity and resilience of the remaining fragments (Mori et al., 2013). Changes in species composition also go along with shifts in functional diversity and biological interactions (Bovendorp et al., 2019;Diaz & Cabido, 2001;Flynn et al., 2011;Petchey & Gaston, 2002, 2007Poos et al., 2009;Swenson, 2014). Eventually, this affects ecosystem services (González et al., 2011;Menon et al., 2007;Rangel, 2000;Torres et al., 2012).
In the recent past, forest cover has increasingly been monitored using remote sensing techniques. For the Colombian high Andean forests, this has shown modest signs of recovery in some areas (Calbi et al., 2020;Etter, 2002;Rubiano et al., 2017;Sánchez-Cuervo et al., 2012; but see Anselm et al., 2018). However, remote sensing cannot detect cryptic forms of forest degradation, such as selective logging or understory grazing. Even plot-based surveys focusing on trees may not reveal such alterations. Yet, cryptic forest degradation has significant impact on soil erosion, successional dynamics, and regeneration, since understory and epiphytic plants are major drivers of ecosystem functioning (Nilsson & Wardle, 2005). Understanding the effects of anthropogenic disturbance on all major forest components, that is, tree, shrub, understory, and epiphyte layers, is therefore essential to elaborate and implement effective strategies for the sustainable management of these forest ecosystems (Battles et al., 2001;Fahey & Puettmann, 2007;Halpern & Spies, 1995;Roberts & Gilliam, 1995). In addition, multiple predictor and response variables should be analyzed simultaneously to properly address disturbance effects within this complex environment.
One of the best areas to study the impact of human-induced alterations on bosques altoandinos in the northern Andes is the area of Bogotá, the capital of Colombia, which is situated at approximately 2,600 m altitude. With a population of around 9 million inhabitants, Bogotá is by far the largest city in the Andean high montane forest belt, putting tremendous pressure on the surrounding ecosystems.
Beyond such floristically oriented approaches, few studies have addressed the effects of disturbance on these forest ecosystems.
Some preliminary research works on forest succession and regeneration were carried out as thesis works (Acuña, 2013;Restrepo Abadia, 2016). In a recent study, Rodríguez-Alarcón et al. (2018) found a negative effect of forest fragmentation on functional diversity and aboveground biomass, a first indication that more complex parameters such as functional diversity are indeed related to ecosystem services such as carbon storage. However, studies that simultaneously consider multiple disturbance predictors and different plant communities response variables were so far lacking.
According to the available literature, the most relevant disturbance factors, which variation proved to be significantly related to differences in forest species composition or diversity metrics, are as follows: age of forest fragment (Köster et al., 2009;Laurance et al., 2006), proximity to houses or roads and people and livestock density (Ribeiro et al., 2015(Ribeiro et al., , 2016, edge effect, and proximity to pastures (Parra Sánchez et al., 2016;Werner & Gradstein, 2009), as well as forest cover fragmentation metrics (Fahrig, 2003;Hertzog et al., 2019;Laurance et al., 2006). Nonetheless, it has not yet been tested whether these factors would be still relevant when a larger number of variables are considered simultaneously. For this reason, we conducted a comprehensive integrated assessment of the potential effects of multiple environmental and disturbance variables on the taxonomic, phylogenetic, and functional diversity of the two main forest layers (tree layer and understory) and on epiphytes cover.
We therefore hypothesized that anthropogenic disturbance as a whole, understood as a composite variable sensu Paine et al. (1998), affects the composition, and aboveground biomass of bosques altoandinos, with impacts on community diversity metrics, that is, taxonomic, phylogenetic, and functional diversity. We also hypothesized that our comprehensive analysis would identify significant predictor and response variables other than those found in previous studies. We specifically set out to answer three questions: (a) Which environmental and disturbance variables best explain K E Y W O R D S aboveground biomass, biodiversity, bosque altoandino, Colombia, cryptic forest degradation, understory species diversity and composition of tree and understory layers?

| Plot setup
Due to the usually small size of forest fragments, we used a plot size of 20 × 20 m (0.04 ha) as established in the framework of the Rastrojos project (Acuña, 2013; Hurtado-Martilletti et al., 2020;Muñoz-Camacho et al., 2017). We complemented the data from the tree layer assessments of 20 plots obtained from the Rastrojos project with data from 12 plots set up and assessed during this study.
In addition to the tree layer data, we also assessed the understory layer, and epiphyte cover in the totaling 32 plots, which are located in six administrative regions of Bogotá D.C. and Cundinamarca ( Figure 1; Appendix A1). We aimed for a widely scattered position of plots in order to represent the landscape (e.g., including differently inclined slopes). Our sampling design was influenced by the distribution of available and accessible fragments. Plot locations belonged to privately owned protected areas and farms, for which we obtained the required permits of entry from the corresponding owners. (http://www.prono stico syale rtas.gov.co/mapas -grafi costiemp o-clima/ indic adore s-clima tolog icos). Mean population density was extracted in two buffers (radius 1 km and 5 km) around the plots from the worldpop database for South America at 1 ha resolution (https://www.world pop.org, Sorichetta et al., 2015). A complete list of all macro-environmental variables can be found in the Appendix A2.

| Tree and shrub layer assessment
Following the protocol of Hurtado-Martilletti et al. (2020), for every woody plant with basal diameter > 5 cm (measured at 5 cm from the ground-DAH: Diameter at "ankle" height), we recorded its DAH, DBH and visually estimated tree height, a method that proved to be quite precise for lower canopies such as the ones studied here (Silva et al., 2012). Plant material was collected and identified with the available literature (Gentry & Vasquez, 1993;Trelease & Yuncker, 1950; or webpages: https://plant asdec olomb ia.com), by comparison with herbarium specimens, digitized specimens available online (JBB: http://herba rio.jbb.gov.co; COL: http://www.biovi rtual.unal.edu. co/en/colle ction s/searc h/plants), or with additional help from local experts. Specimens were deposited in the herbarium of the Jardín Botánico de Bogotá José Celestino Mutis (JBB); high-resolution digital specimen images can be provided upon request; and a plotresolved list of vouchers can be found in the Appendix A3.

| Understory assessment
In each 20 × 20 m plot, eight 1 × 1 m quadrants (with marked 10 cm subgrids) were placed randomly. All vascular plants, including tree seedlings, were recorded, and mean height and total cover (the sum of all individuals cover) were measured for every species in each quadrant. When available, fertile material was collected and deposited in the JBB. Additionally, cover of bare soil, leaf litter, bryophytes, lichens, and coarse woody debris was visually estimated for every quadrant.

| Epiphyte cover
In each plot, we sampled 40 randomly selected trees to estimate the epiphyte cover. Categorical cover classes (ranging from 0 to 3) were assigned to each of five major epiphyte groups (bryophytes, lichens, ferns, bromeliads, and orchids), separately for trunk and canopy branches.

| Functional traits and functional diversity
Three leaf functional traits (specific leaf area: SLA; leaf thickness: LT; and leaf dry matter content: LDMC) were measured for each tree species following the protocols provided by Pérez-Harguindeguy et al., (2013). Five leaves were collected from each of up to three different individuals per species and stored in wet paper for at least 12 hr, then weighted (petiole included). LT was measured with a digital micrometer, and a digital scan of the fresh leaves was taken with a Hewlett-Packard F4280 scanner. Leaf area was calculated with ImageJ 1.8.0 (Schneider et al., 2012). Leaves were oven-dried at 60°C until constant weight and weighted; SLA was then calculated as one-sided area of a fresh leaf divided by its dry mass, expressed in cm 2 /g. LDMC was calculated as the dry mass (mg) divided by its fully hydrated fresh mass (g), and expressed in mg/g. Additionally, wood density (WD) was obtained from Rodríguez-Alarcón et al. (2018) and the global wood density database (Chave et al., 2009) et al. (2008), using the R package FD (Laliberté & Legendre, 2010;Laliberté et al., 2014). We specified "corr + lingoes, m = 3" to reduce dimensionality. Functional diversity (FD) index (Petchey & Gaston, 2002) was calculated as the total branch length of a functional dendrogram generated on a distance matrix of traits with the R function hclust, using the PD function in the R package picante (Kembel et al., 2010). We decided to compute functional diversity according to the framework proposed by Mason et al. (2005) and Villéger et al. (2008). The calculated indices provide independent information about the position and relative abundances of species in a multidimensional functional space, allowing for a more detailed examination of the mechanisms linking biodiversity to ecosystem function (Villéger et al., 2008).

| Landscape metrics
A Landsat 8 raster was downloaded from the US Geological Survey and processed in QGIS with the SCP plugin (Congedo, 2016) to obtain a land cover map. Landscape metrics refer to the size, shape, configuration, number, and position of land-use patches within a landscape and were obtained for the forest class within a 1,000 m diameter buffer zone around the plots with the LecoS plugin (Jung, 2013).
Additionally, fragments of forests were manually vectorized and the area was calculated on a prepared Bing aerial map obtained through the Openlayers plugin (see Figure 1 for an example).

| Community composition and structural variables
Based on the available literature (Cleef, 1981;Cortés, 2008;Cuatrecasas, 1958;Sturm & Rangel, 1985;Van der Hammen, 1998), tree layer species were classified either as late successional slowgrowing, early successional fast-growing, exotic, or "other" (see Appendix A3 for details). Additionally, understory exotic species cover was calculated. The number of species and the relative proportion of individuals (in case of trees) or the percent cover (in case of the understory) of exotic species were used as indicators of disturbance versus conservation. Variance of tree DBH and height was also computed across all trees within each plot, together with the overall number of tree individuals, stems, stems per tree, and the percentage of large trees (DBH > 30 cm). Mean understory height and cover was calculated, as well as mean epiphytes cover.
The Gini coefficient, a measure of inequality within a distribution widely used in forestry (Bourdier et al., 2016;Latham et al., 1998;Lexerød & Eid, 2006), was calculated in each plot for stem basal areas with the gini function in the R package reldist (Handcock, 2016).

| Taxonomic and phylogenetic diversity
Alpha-diversity indices (Shannon's diversity, Simpson's and Pielou's evenness) were computed for each plot with the R package vegan (Oksanen et al., 2013). Phylogenetic community structure was assessed on the basis of a published angiosperm supertree (Phylomatic tree R20120829, available at https://github. com/camwe bb/tree-of-trees/ blob/maste r/megat rees/R2012 0829.new). First, a regional pool tree was generated with the Phylomatic webtool (Webb & Donoghue, 2005), and then, branch lengths were assigned with the bladj algorithm in the software Phylocom 4.2 (Webb et al., 2008), using the wikstrom.ages file (Wikström et al., 2001). Phylogenetic diversity (PD), mean pairwise distance (MPD), mean nearest taxon distance (MNTD), and their standardized counterparts (sesPD, sesMPD, and sesMNTD) were calculated for both trees and understory in the R package picante (Kembel et al., 2010). Moreover, abundance-weighted MPD and MNTD were calculated to account for differences in species abundance (Webb et al., 2011). Four species of the Lycopodiaceae had to be removed from the understory regional pool since the family was not included in the used supertree.
The standardized PD metrics express the difference between observed and average value in units of standard deviation (SD).
Positive values indicate phylogenetic overdispersion (co-occurring species are more distantly related than expected by chance) and negative values phylogenetic clustering (co-occurring species are more closely related than expected by chance). were considered a potential cause of disturbance, whereas the nearest distance to a road was considered a potential facilitator.

| Aboveground biomass
Diversity indices and biomass estimation were included in the level category.
To filter for dominant variables, we first ordinated the plots based on relative species abundances using nonmetric multidimensional scaling (NMDS) in vegan, using the mds function with Bray-Curtis distances, and specifying three as maximum number of axes.
Subsequently, we fitted all variables using the envfit function and examined variable ordination scores, in order to identify the variables most strongly correlated with community composition and to assess redundancy. Both predictor and response variables were included in the same analyses, and NMDS was performed separately for tree and understory layer. Second, using Sørensen distances and flexible beta (set to -0.25) as group linkage method (McCune & Mefford, 2015), cluster analysis and subsequently indicator species analysis for each cluster were carried out in PCORD 7 (McCune & Mefford, 2015), in order to further classify community types and their characteristic elements.
Following this preliminary analysis, we determined a subset of variables that correlated with the main axes above a given threshold (R sq > 0.35; Table 1) and performed either Kruskal-Wallis or parametric ANOVA, depending on the determined conditional distribution, using the clusters as independent variables and the filtered subsets of variables as response variables.
To verify the presence of spatial autocorrelation in our predictors and responses, we calculated a geographical distance matrix between study plots and performed Moran's I test for all calculated variables. We detected spatial autocorrelation for 24 predictors, but none in our response variables (i.e., diversity metrics).
Finally, partial redundancy analysis (pRDA) was performed in vegan, separately for tree and understory layer. To take into account spatial autocorrelation, we fitted a pRDA specifying "locality" as condition, to be able to rule out locality effect on the ordination.
The "condition" argument thereby defines partial terms that are fitted before other constraints and can be used to remove the effects  (Oksanen et al., 2013).
Additionally, we performed Hellinger transformation of our community data as recommended by Legendre and Gallagher (2001).
We further selected predictors from the set obtained with the NMDS screening, by checking for correlation (r > 0.7), performing Variance Inflation Factor (VIF) analysis (setting the threshold to 10), and then using the vegan ordistep function which performs automatic stepwise model building for constrained ordination methods (Oksanen et al., 2013).

| General linearized models between main causes and facilitators of disturbance and main response variables
To select meaningful variables to fit our GLMs, we inspected the NMDS and pRDA graph and selected a set of uncorrelated response variables based on the direction of the arrows in the graphs. We then compiled a set of predictor variables that correlated with each selected response and checked for correlation within each set, removing one of the elements in pairs with r > 0.7. In parallel, we merged all predictor sets and removed highly correlated and spatially autocorrelated variables. Once a set of consensus predictors was obtained, we conducted a VIF analysis (setting the threshold to 10) and obtained a reduced set of primary and secondary predictors (Table 2). We thus reduced the pool of geo-environmental variables to four, that of causes to four, and that of facilitators to seven. In addition, we selected response variables for level, including diversity metrics and indicators.
For each selected response, we identified the best conditional distribution and then performed automated selection of the optimal Generalized Linear Model (GLM) with the regsubsets function in the leaps package (Lumley & Lumley, 2013), unifying all groups of predictor variables. We specified a maximum number of predictors of four. Predictors were scaled, and a "log" link was specified in the family argument. were associated with tree layer group 1 which was also positively correlated with Shannon's landscape diversity and negatively with elevation. Group 4 was defined by lower values of Shannon's landscape diversity and was positively correlated with minimum fragment age. Group 5 had some degree of negative correlation with minimum fragment age. Group 6 had an inverse correlation with elevation and minimum fragment age, and was associated with signs of logging, higher Shannon's landscape diversity, and absence of cattle.
Groups 2 and 3 were not characterized by any particular association with the ordination variables ( Figure 2).

| Understory layer
From the set of 38 variables with R sq > 0.35 (Table 1), after the assessment of redundancy, we limited our analysis to a subset of 10: elevation, number of trees, edge density, Shannon's landscape diversity, mean tree aboveground biomass (mAGBT), presence of cultivated fields in a 500 m buffer, minimum fragment age, distance to roads, people density in a 5 km buffer, and fragment size. The ordistep function selected seven of these: elevation, edge density, Shannon's landscape diversity, mAGBT, presence of cultivated fields in a 500 m buffer, distance to roads, and people density in a 5 km buffer.
The pRDA had an R sq of 0.26 and adjusted R sq of 0.11. The proportional conditional explained variance was 0.29, while the constrained explained variance was 0.26. The unconstrained was 0.45.
The results of the pRDA indicated that group 1 was characterized by higher values of Shannon's landscape diversity, lower values for distance from roads, lower elevation, and lower edge density. In contrast with that, group 2 was linked with higher values for distance from roads, higher elevation, lower Shannon's landscape diversity, absence of cultivated fields in a 500 m buffer, and lower population density. Group 3 had lower values of mean tree biomass and higher values of population density. Group 5 was associated with lower values of mean tree biomass. Group 4 was not characterized by any particular association with the ordination variables ( Figure 3).

| Generalized linear models
A total of 15 primary predictors, 17 secondary predictors, and 17 responses were retained for GLM building (  Table S1. None of the variables associated with epiphytes cover was retained through the variable selection process and analysis. Tree layer Shannon's diversity decreased with slope. Understory Shannon's diversity increased with distance to tracks, mean precipitation, and tree layer functional richness (FRic), but decreased with minimum age and functional dispersion (FDis). Understory abundance-weighted phylogenetic mean pairwise distances (HsesMPDABU) increased with functional divergence (FDiv) and tree layer abundance-weighted phylogenetic mean pairwise distances

| D ISCUSS I ON
Pressure of urbanization on natural environments and its consequences has been the subject of numerous studies. However, high Andean forests (bosques altoandinos) have rarely been investigated in this context. Our study is the first to analyze the role of multiple factors in shaping environmental impact on these forests through urbanization and associated factors in the metropolitan area of Bogotá. However, we are aware of the limitations of this research, which is of rather explorative character and based on data from an area of in total 1.28 ha only. Our sampling design reflects the hurdles F I G U R E 3 pRDA and Cluster analysis convex hulls of the understory. RDA graphs with convex hull volumes of understory groups for axis 1-2 (a), 2-3 (b), and 1-3 (c  Also, having a limited number of plots, we decided to put a stronger emphasis on the variables filtering, to drastically reduce the number of tested hypotheses. Nevertheless, the studied forest fragments belong to several of the localities harboring the highest forest cover within the Capital District and we find the types of high Andean forests covered here to be representative for the hinterland of Bogotá. Using the composition of natural vegetation as a benchmark, our study plots were dominated by Melastomataceae, Ericaceae, and Asteraceae in the tree layer, which is in accordance with previous work (Cuatrecasas, 1934(Cuatrecasas, , 1958Franco et al., 2010;Torres & Marina, 2016). Bromeliaceae and Orchidaceae were the most diverse families in the understory, coinciding with reports by Cuatrecasas (1934Cuatrecasas ( , 1958 and Rangel et al. (2008). Notably, with the exception of Rangel et al. (2008), no recent inventories of the understory were undertaken in the target area prior to this study. The fact that many epiphytic species were found terrestrial in the understory may be due to certain favorable environmental conditions, such as low incidence of light, high humidity, and lower influence of wind than in the canopy (Krömer et al., 2007).
Overall tree species richness of the total area assessed (98)  and conspicuous lianas in the understory. This community corresponds to our tree clusters 1 and 6 and understory clusters 1 and 5.
The pRDA further revealed a lower elevation, higher Shannon's landscape diversity, lower minimum fragment age, presence of logging and lower distance to roads as characteristic for this community, supporting the notion that it represents secondary forest, probably developing on patches of abandoned agricultural areas on the slopes surrounding cultivated and farmed plains (Cortés, 2008). Understory cluster 5 was generally found at medium elevations, on small high plains, with a drier climate (Cortés, 2008;Cortés et al., 1999), and in forest patches with generally low values of aboveground biomass.
The Drimys granadensis-Weinmannia tomentosa community is a second bosque altoandino subtype (Vargas & Zuluaga, 1980), corresponding to our tree clusters 2 and 4. Cluster 2 is similar to the Criotoniopsis bogotana-Weinmannia tomentosa forest subtypes described for elevations between 3,100 and 3,300 m (Cortés, 2008), whereas cluster 4 is found at the slopes and peaks of the watershed of the Río Bogotá between 2,700 and 3,200 m (Cortés, 2008). According to Cortés (2008) and Luteyn (2002), the presence of Macleania rupestris in the lower canopy of these communities points toward recent human intervention. This association is known to prefer humid, cold climates and steep grounds; according to our field observations, it is also associated with high lichen and moss cover in the canopy, which prosper in such a relatively high humidity (Batke et al., 2015;Munzi et al., 2014;Wolf, 1993). As shown in the pRDA ordination, it is also linked to low Shannon's landscape diversity, and higher minimum fragment age, probably representing secondary forest fragments approaching the structure of natural forest communities.
Our tree clusters 3 and 5 did not correspond to previously described communities. Cluster 3 was restricted to bosques altoandinos near Sumapaz, the largest known paramo on Earth. Characteristic elements of this cluster are families of high elevations such as Asteraceae and Ericaceae (Bach et al., 2007;Cuatrecasas, 1958;Sturm & Rangel, 1985), also typically found in areas subjected to fires or selective logging (Cuatrecasas, 1958). The latter notion is opy, and a dense cover of mosses and ferns, which suggests that some small "islands" of mature forest elements were able to persist within the disturbed, secondary forest matrix. Understory cluster 3 did not fit any previously described communities either, but the indicator species of this cluster are known to be either dispersed by birds, for example, Monnina aestuans (Romero, 2002) and Nertera granadensis (Vargas-Ríos, 1997), or by small mammals or birds, for example, in the case of the sticky fruits of Peperomia (Frenzke et al., 2016). Possibly, this cluster represents a successional understory community mainly dispersed by animals, which prosper in previously disturbed areas, as suggested by the high people density within 5 km radius and relatively low mean tree biomass. Tree cluster 5 was found in the Guasca region only and exhibits features of a disturbed, gap-filled forest (azonal páramo) including the presence of invasive Ulex europaeus, which is confirmed by the pRDA correlation with lower minimum fragment age values. Another common species, Cavendishia bracteata, has been associated with secondary growth (Cortés, 2008). This cluster had rather low like adjacencies values and average Shannon's landscape diversity and distances to roads, which point to a somehow continued disturbance regime in the past.
Indeed, this area, up to the 1990s, used to be an open-pit limestone mine (Pèrez Sanz de Santamaría, 2013).
Notably, tree and understory communities found in the same plots did not always correspond to the same community's type, which suggests that different types of intervention act differentially on the tree and understory layers. For instance, cattle grazing, erosion, and expansion of edge species will affect the understory at a different pace than the tree layer (Halpern & Lutz, 2013;Millspaugh & Thompson, 2011;Thrippleton et al., 2016).
Our findings support the notion that bosques altoandinos in the vicinity of Bogotá are floristically and structurally not homogeneous, resulting in overall high species diversity, especially in the understory, with each of the study sites and plots contributing a portion to this diversity (i.e., high beta diversity). The observed differences in species composition between the study sites, and the high proportion of pRDA-explained variance that was linked to the "locality" condition, may be determined by topographic variation, which promotes changes in structure, composition, and dynamics of the vegetation, even at small scales in high Andean ecosystems (Homeier et al., 2010;López & Duque, 2010). Our results are similar to a recent study that found substantial differences in species composition between municipalities in the region (Hurtado-Martilletti et al., 2020), The compositionally based clustering of tree and understory communities was largely correlated with both geo-environmental and disturbance variables, namely, elevation, people density, Shannon's landscape diversity and distance to roads. Mean temperature, relative humidity, logging, and minimum plot age were important factors driving tree species composition, but not the composition of understory species. For the latter, additional variables associated with edge effects, such as the proximity to cultivated fields, edge density, and distance from main roads were relevant. Additionally, mean tree aboveground biomass was a determinant factor in shaping the understory community. These results support the notion of a higher sensitivity of the understory to fragmentation and habitat heterogeneity (Forman & Alexander, 1998;Tyser & Worley, 1992).
Our results show effects of both geo-environmental parameters and disturbance-related variables as predictors of both community structure and diversity. Among the geo-environmental parameters, the negative effects of the increase in slope on tree and understory diversity and aboveground biomass were evident. Slope is related to soil erosion, water drainage, and other unfavorable growth conditions which may act as environmental filters, reducing the number of taxa that can cope with them effectively and may also limit aboveground productivity. Higher mean temperatures were linked to lower tree Pielou's evenness and Understory phylogenetic diversity.
This fact could be linked to the higher density of human activities at milder temperatures/lower parts of our study area, which are associated with highly disturbed forest communities, mostly dominated by species as Miconia squamulosa or Cavendishia bracteata, and host poorer understory communities. Higher precipitation values were linked to higher understory Shannon's diversity, possibly due to increased soil nutrients and moisture and thus by the absence of an environmental filter related to water availability.
With regard to human disturbance predictors, many of the previously identified relevant variables in literature were also selected through our multi-step analysis, such as minimum age of the forest fragment, distance to houses, edge effect, and presence or proximity of cattle and cultivated fields. People density, on the other hand, showed to be too spatially autocorrelated to be used in our GLMs.
Also, among all calculated forest fragmentation metrics, the only one which was selected was (median) forest patch size, already reported to be relevant for plant diversity as an indirect measure of habitat loss in the review of Fahrig (2003). As to the selected responses, tree layer diversity metrics were not particularly sensitive, retrieving only one GLM with a good fit. The correlation between higher distance from houses and forest protection status with lower tree species richness and low phylogenetic diversity was not immediately intuitive, but could be a sign of the deliberate introduction of useful tree species in the vicinity of rural houses, to be harvested for wood or other uses, or of the lack of edge-related tree species in the interior of protected forest fragments. However, the presence of cattle and cultivated fields in the immediate proximity of plots leading to tree phylogenetic clustering, but on the other hand to understory phylogenetic dispersion, illustrates the disrupting, multi-layer impact of landscape-level patchiness and human activities.
Disturbed forests tend to exhibit functional and phylogenetical clustering due to the elimination of entire lineages sensible to disturbance, an effect known as environmental filtering (Chun & Lee, 2018;Gerhold et al., 2015;Kusuma et al., 2018;Mouchet et al., 2010;Ribeiro et al., 2016). Phylogenetic dispersion is expected to be higher in undisturbed, more mature forests than in early successional forests, due to competitive exclusion (Ding et al., 2012;Letcher, 2009;Norden et al., 2012;Purschke et al., 2013). In our study, local, chronic disturbances, such as proximity to farming activities or the presence of cattle in the immediate surroundings, had indeed a negative effect on tree phylogenetic diversity and resulted in phylogenetic clustering, supporting findings by Ribeiro et al. (2015Ribeiro et al. ( , 2016. Likely, the floristic drift associated with this type of disturbance results in the co-occurrence of more closely related taxa by decreasing effects of competitive exclusion. On the other hand, the observed increase of phylogenetic dispersion in the understory in close proximity of cattle or cultivated fields may be the result of opportunistic pioneer or exotic species, which introduce different lineages from those associated with more mature forest fragments (Hill & Curran, 2001;Kupfer et al., 2004).
Identified understory diversity metrics with the highest sensitivity to human disturbance were Shannon's diversity and phylogenetic clustering. As suggested by Forman and Alexander (1998) and Tyser and Worley (1992), the number and diversity of understory species were positively related to disturbance-related variables. Proximity to human activities such as farming and the more recent establishment of forest patches (lower minimum age) fosters generalists or fast-growing, nutrient-, and light-demanding species (Marcantonio et al., 2013). However, at the same time the edge effect promotes less phylogenetic diversity of the understory vegetation, which is in accordance with Ribeiro et al. (2016). This could be explained, in our case, by the fact that ferns and other early diverging taxa diversity tends to diminish toward the edge of a forest fragment to leave place to generalists and agricultural weeds, which can cope better with the site conditions. Larger median forest fragments size also resulted in less understory species, suggesting that recruitment of edge-related species increment the number of species in smaller forest patches.
The observation that increasing tree functional divergence, and tree phylogenetic dispersion were linked to higher understory phylogenetic dispersion, may indicate that higher trait diversity in the upper stratum allows for more species to colonize the understory. 112 Mg/ha for forest plots within a similar elevation range. In regard to our models, AGB seemed to decrease at higher values of slope, which in our study area may relate to eroded soils and drier conditions, supporting a trend that has been reported for relatively moist forests in the Americas (Keith et al., 2009;Stegen et al., 2011), which is perhaps related to the lower soil water content available to sustain photosynthesis (Parton et al., 2012;Stegen et al., 2011), but that can also be a secondary effect of the different rate of agricultural exploitation or forest clearing history between lower and drier and higher and wetter soils in the study area in recent times (Etter et al., 2008;Etter & van Wyngaarden, 2000). Notably, low AGB was linked to the proximity of cultivated fields, suggesting a clear correlation between intervention causing patchy landscapes and lower biomass accumulation. However, the presence of cattle within the plot was linked to higher AGB values. This may be particular to our study area, in which we observed forest fragments with large trees but a much depauperate understory, located in proximity to farms.
This is alarming as grazing may interfere with tree species recruitment and stamping may lead to higher soil erosion which in turn will reduce productivity over time in these last standing carbon stock fragments (Nepstad et al., 2002). In conclusion, the increase of disturbance resulted overall in a negative effect on tree phylogenetic diversity and dispersion.
Notably, disturbance affected aboveground biomass negatively. As to the understory, disturbance was associated with more diversity and more phylogenetic dispersion. The causes and the facilitators category variables were quite efficient in predicting diversity or AGB, among which edge effect, proximity of cattle and cultivated fields, and minimum fragment age appear to be the most important ones.
The plurality of diversity metrics can be difficult to interpret in the light of human disturbance. However, AGB proved to be sensitive to human disturbance and was closely related with the proportion of late successional species. Such indicators could serve as immediate proxies of human disturbance, rather than diversity measures themselves, which have also been shown to react ambiguously to the effects of fragmentation (Fahrig, 2003).

| CON CLUS IONS
In summary, our study on taxonomic, phylogenetic, functional diversity and ABG of high Andean forest underscores the complexity and singularity of interactions between disturbance drivers and plant communities. The main goal of our approach was to test and quantify the alteration of high Andean forest composition, structure, and functioning through human disturbance, testing the effectiveness of known relevant drivers and indicators when a large number of variables are considered simultaneously. We contributed to the characterization of high Andean patterns of tree and understory diversity and local and regional human disturbance, which is usually considered to have a negative effect on native biodiversity and carbon storage.
In our case, this fact was confirmed by lower tree layer diversity and a lower ABG in relation to increasing human disturbance, but was however not always apparent through the score of all the diversity metrics that we employed. Decline of AGB and disappearance of the forest ecosystem's late successional species is a warning signal that should impulse protection efforts and restoration measures. Yet, it is also true that the study area has now undergone anthropic disturbance over centuries, with continuous agropastoral activities and subsequent land cover change. In the context of the recovery of forest cover and ecosystem services, then our findings could be interpreted as a positive sign of resilience at a regional scale. Relatively small isolated fragments of high Andean forests can still host high plant diversity and serve as stepping stones or temporary refuges for the local fauna within the rural modified matrix. In this sense, efforts to implement forest connectivity and corridors and to guarantee land-use continuity even in partially forested areas are priorities that should be taken into account by local decision-makers. Successful conservation strategies require a sound understanding of community and ecosystem dynamics, and we hope that with the predictors and indicators of disturbance that we pointed out, it will be possible to improve the management strategies for the passive or active restoration and protection of the remaining forest fragments in the study area.
Our results contribute to urgently needed but yet missing baseline knowledge on main drivers of disturbance and its effects on the biodiversity in the study area. However, we strongly recommend that future studies should expand further the established plot network and that more investigations test our results on similar ecosystems to further disentangle the relationship between natural and human-induced causes of diversity loss and their underlying mechanisms. As shown here, a first approximation can be achieved through an exploratory approach like the one that we employed.

ACK N OWLED G M ENTS
We thank our colleagues at the herbarium of the Jardín Botánico de Bogotá José Celestino Mutis (JBB) for the logistic support, their help in the identification of plant material, and valuable comments.
We are also grateful for the support from local contacts in Torca, Sumapaz, and Pasquilla. We thank the Sintrapaz association and the Colegio Nuevo Horizonte for their kindness and collaboration. Finally, we thank everybody who participated in field data collection and collation of this study. We also thank Ana Belén Hurtado-Martilletti and Natalia Norden for their help and feedback on earlier versions of this study. We are grateful for the valuable comments from two anonymous reviewers that helped to improve this manuscript.

CO N FLI C T O F I NTE R E S T
None declared. writing-review and editing (supporting).

DATA AVA I L A B I L I T Y S TAT E M E N T
The tree layer and understory sampling datasets and the complete

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.

A PPE N D I X A 5
I n d ic ato r s p e cie s a n a l ys i s v a l u e s (I V I) fo r t r e e a n d u n d e r s to r y l aye r

A PPE N D I X A 5 (Continued)
A PPE N D I X A6 N M DS t r e e l aye r g r a p h s a n d a n a l ys i s of v a r ia n ce b ox p l ot s NMDS graphs of the tree layer. (a) Ordination graph of plots in tree species space for axis 1-2: cluster analysis groups are outlined.
(b) Ordination graph of plots in tree species space with plotted variables for axis 1-2, with only the variables with 'p.max = 0.05' plotted.
The variables that most correlated with the ordination axes (R Sq > 0.35 for any of the 3 ordination axis, Table 1) and therefore with species composition and abundances are depicted clockwise. For full variable names and acronyms please refer to Appendix A2 .
Boxplot of analysis of variance for tree layer groups variables. Only variables that significantly differed among various tree groups ( Figure 2) are shown. From right to left: age; presence of cattle inside the plot; (c) presence of cattle within 50 m from the plot; presence of cultivated fields within 100 m from the plot; elevation; Shannon's landscape diversity in 1 km buffer; Like adjacencies in 1 km buffer; logging; mean annual temperature; relative humidity.

A PPE N D I X A7
N M DS v a r ia b l e s co r r e l at io n w it h o r d i n at io n a xe s fo r t r e e a n d u n d e r s to r y l aye r

A PPE N D I X A7 (Continued)
A PPE N D I X A8 N M DS u n d e r s to r y l aye r g r a p h s a n d a n a l ys i s of v a r ia n ce b ox p l ot s