1. Introduction

Biodiversity distribution patterns and their causes have been central to natural history and ecology since they were first conceived (Caldas, 1803; Humboldt and Bonpland, 1805) until the present (Storch et al., 2005; Sanmartín, 2012). Mountains provide invaluable insight for developing a general theory on biological diversity to explain these patterns (Körner, 2004; Sanders and Rahbek, 2012). Compared to latitudinal gradients, elevational gradients encompass more climate regimes, greater topographic complexity, and wider environmental oscillations over shorter geographical distances (Antonelli et al., 2018; Rahbek et al., 2019a). Thus, mountains are exceptionally biodiverse: where they comprise only 25% of Earth's terrestrial surface but harbour ca. 87% of terrestrial biodiversity, with a large proportion of endemicity (Rahbek et al., 2019b; Wambulwa et al., 2021). In particular, Andean tropical montane forests (TMFs) host one-sixth of all plant species, where 44% of which are endemic (Myers et al., 2000; Hughes and Eastwood, 2006). Therefore, Andean TMFs are among the biodiversity hotspots with highest conservation priority (Mittermeier et al., 2005). Lastly, from a biogeographic perspective, Andean TMFs are unique natural laboratories that provide rigorous tests of diversity trends along altitude gradients and comparisons of these trends among sites (Lomolino, 2001; Guo et al., 2013; Graham et al., 2014) because they allow to replicate ample elevational ranges over a large continental north-south corridor that has witnessed historic processes over millions of years (Antonelli et al., 2009; Malhi et al., 2010; Tito et al., 2020).

Biodiversity has traditionally been explored in terms of taxonomic diversity by treating all Linnaean binomials as ecologically equivalent and evolutionary independent. This approach cannot yield much information about the ultimate drivers of biodiversity distribution patterns (Pavoine and Bonsall, 2011). A more accurate depiction of diversity should also consider functional and phylogenetic diversity (McGill et al., 2006; Swenson, 2011). Functional diversity uses species ecological, morphological, and physiological features to understand where species can thrive, how they interact, or how they contribute to ecosystem functioning (Petchey and Gaston, 2006; Violle et al., 2007; Cadotte et al., 2011). Phylogenetic diversity tracks the origins of taxa lineages through deep historical and evolutionary processes that shaped their distribution patterns and thereby provides insights into the key events in biodiversity distribution (Cavender-Bares et al., 2009; Gerhold et al., 2018). Although these three diversity facets are intrinsically related (e.g., closely related taxa tend to share more traits; Moles et al., 2005) their responses to the environment can differ (Devictor et al., 2010; Cisneros et al., 2014). However, few studies have simultaneously considered the three types of diversity, particularly for plants under changing environmental conditions, such as those found along elevational gradients (but see Tanaka and Sato, 2015; Chun and Lee, 2017). Ignoring one diversity facet not only restricts our understanding of the causes of biodiversity distribution but also hinders conservation strategies, as focusing only in one dimension may lead to the neglect of unique species, functional roles and/or evolutionary histories (Brum et al., 2017).

Previous studies on biodiversity patterns along elevational gradients in different biomes and distinct taxonomic groups have arrived at contrasting results regarding taxonomic diversity, but largely agree on functional and phylogenetic trends. Historically, the prevalent pattern for taxonomic diversity is a decrease with elevation, although some studies have described species richness peaks at medium elevations (Rahbek, 1995; McCain, 2005; Dani et al., 2023). The most common pattern for functional diversity also comprises narrowing towards fewer functional strategies to cope with harsher environmental conditions as elevation increases (Swenson et al., 2012; de Bello et al., 2013; Thakur and Chawla, 2019). The elevational pattern for phylogenetic diversity is specific to the key biogeographic historical events of each site (Ricklefs and Schluter, 1993). In the particular case of the Andes, two hypotheses have been proposed to explain the historical formation of biological communities, reflecting distinct patterns of phylogenetic diversity: a) the establishment of Amazonian clades originated under warm conditions at tropical lowlands –typically exhibiting strong conservatism of ancestral narrow thermal tolerances– which, during the Pleistocene, adapted to the comparatively cooler temperatures of the newly raised highlands resulting from the Andean uplift (Latham and Ricklefs, 1993; Kerkhoff et al., 2014); and b) the arrival of cold pre-adapted lineages –originated much earlier in temperate latitudes– via migration along the Andean high-elevation corridor (Jablonski et al., 2006). The former (a) is the cornerstone of the Tropical Niche Conservatism (TNC, Wiens and Donoghue, 2004) hypothesis, which predicts an elevational decrease in phylogenetic diversity (Tolmos et al., 2022), whereas the latter (b) is the core of the Out of the Tropical Lowlands (OTL, Qian and Ricklefs, 2016) hypothesis, which predicts an elevational increase in phylogenetic diversity (González-Caro et al., 2020).

In this study, we investigated plant diversity in Andean TMFs at four transects ranging elevational gradients from ca. 800–3100 m a.s.l. and spanning ca. 1750 km in latitude across Ecuador, Peru, and Bolivia. We used census data from 114 plots (containing ca. 38, 000 individuals and 2000 taxa) and analysed them using Hill numbers to quantify different types of abundance-based diversity indices. We aimed to 1) assess the elevational trends in taxonomic, functional, and phylogenetic diversity, 2) examine the relationships between different diversity facets along elevation and 3) explore their underlying ecological and evolutionary causing processes. To the best of our knowledge, this is one of the most ambitious elevational studies of the three facets of plant diversity and a similar study has yet to be conducted along the Andean range.

2. Materials and methods 2.1. Study area

We studied four transects ranging elevational gradients of Andean TMFs located in remote regions: Podocarpus National Park (Ecuador), Río Abiseo National Park (Peru), Madidi National Park (Bolivia), and Pilón-Lajas Biosphere Reserve (Bolivia) (Fig. 1). All four transects contained a similar elevational range (ca. 2200 m), within which we defined three elevational belts: premontane (900–1300 m), montane (1800–2250 m), and upper (2700–3100 m). Within each belt, we established nine to ten plots of 0.1 ha (50 m × 20 m) as described by Arellano et al. (2016), with a total of 114 plots (Tables 1 and S1). Plots were located in mature and apparently undisturbed forests at distances of at least 300 m. In each plot we inventoried all woody stems (including trees, shrubs, lianas, hemiepiphytes, palms, and tree ferns, but excluding woody Poaceae) with a diameter at breast height ≥2.5 cm, with a total of 37, 869 individuals belonging to 2021 taxa (71% individuals identified at species level, 18% at genera and 11% at family) (Table 1, Appendix A Supplementary data).

Environmental characterization of the plots was conducted using bioclimatic variables from the CHELSA climatological dataset (Karger et al., 2017) (Table S2). Mean annual temperature is the most representative variable for studies along elevational gradients because it represents the fundamental environmental change associated with elevation (Körner, 2007). In our study sites, mean annual temperature exhibits high correlation with elevation (r > 0.98) and the thermal ranges of transects vary slightly: 12 ℃ in Podocarpus, 9 ℃ in Río Abiseo, and 10 ℃ in both Madidi and Pilón-Lajas. Alternative climatic variables related to precipitation are not as relevant for our study because 1) moist TMFs, characterized by persistent fog and high precipitation throughout most of the year, are not subjected to relevant seasonal water stress, and 2) at each transect, plots were located along the same river basin to minimize spatial rainfall variability inherent to the complex topography of the Andes, thus leaving temperature as the key source of abiotic variation along our elevational gradients (Laraque et al., 2007; Malhi et al., 2010). Small changes in mean annual temperature are expected to have a major influence on species distributions, especially in the tropics where species have narrow thermal tolerances (Janzen, 1967).

2.2. Taxonomic and functional characterization

We collected voucher specimens in the field, which were later identified and deposited in local herbaria (HUTPL and LPB [Thiers, 2022]). Taxonomic standardization was conducted using the 'Taxonstand' R package (Cayuela et al., 2012; v.2.2). A summary of taxonomic composition for elevational belts can be found in Table 1 (specific details for each plot are available at Table S1).

Functional characterization was only conducted for two transects: Podocarpus (Ecuador) and Río Abiseo (Peru). For each taxon, we measured the specific leaf area, leaf thickness, and branch wood density as described by Cornelissen et al. (2003). These functional traits have been widely used for assessing the key woody plant functional strategy axes (Chave et al., 2009; Wright et al., 2004). We calculated mean trait values for each taxon (Báez et al., 2022). On average for each plot, functional data for the three traits were available for 65 taxa (86%) in Podocarpus and 63 taxa (92%) in Río Abiseo.

2.3. Phylogenetic tree generation

We generated a phylogenetic tree (Appendix A) using the 'V.PhyloMaker2' R package (Jin and Qian, 2022), which allowed us to prune our species list from the largest currently available dated mega-phylogeny for vascular plants.

For phylogenetic diversity calculations, we only considered taxa identified to the species or genus level, i.e., 703 taxa (88%) in Podocarpus, 416 taxa (71%) in Río Abiseo, 453 taxa (97%) in Madidi, and 457 taxa (98%) in Pilón-Lajas. To include taxa absent from the original mega-phylogeny in the tree, we joined them at the half-way point of the genus/family branch (for species/genera) using the parameter 'scenarios' = S3 in the 'phylo.maker' function.

2.4. Diversity measurements

Attribute diversity comprises a unified framework for assessing biodiversity by integrating three diversity facets: taxonomic, functional, and phylogenetic diversity (Table S3) (Chao et al., 2014a). Each type of diversity is measured in different entities (species for taxonomic, species pairwise Euclidean distances for functional, and species tree branch length for phylogenetic), but all three can be transformed into directly comparable units by using generalized Hill numbers. Hill numbers (Hill, 1973) quantify abundance-based taxonomic diversity using the parameter, q, which defines the sensitivity to species abundances (Jost, 2006; Chao et al., 2014b). For q = 0, all taxa are treated equally regardless of their abundance (i.e., taxa richness), whereas larger q values give greater importance to the relative abundances of taxa (i.e., more weight to abundant/dominant taxa and less weight to rare taxa). The attribute diversity framework extends the usage of Hill numbers to functional (Chiu and Chao, 2014) and phylogenetic (Chao et al., 2010) diversity by weighting each entity based on either its relative functional distance compared with other taxa (functional diversity) or the relative length of the phylogenetic tree branch for that taxon (phylogenetic diversity). We calculated the attribute diversity indexes for Hill numbers of the orders q = 0, q = 1 and q = 2 at the plot level using the 'renyi' function of the 'vegan' R package (Oksanen et al., 2019; v.2.5–6) for taxonomic diversity, the code from Chiu and Chao (2014) for functional diversity and the 'ChaoPD' function of the 'entropart' R package (Chao et al., 2010; v.1.6–10) for phylogenetic diversity.

2.5. Statistical analyses

First, to acknowledge the influence of sample size on diversity metrics we conducted Pearson's correlations between number of trees per plot and the three facets of diversity (taxonomic, functional, and phylogenetic) for different Hill orders.

Next, we explored Pearson's correlations among the three types of diversity for Hill numbers of different orders.

Then, we used generalized linear models (GLMs) to investigate the influence of elevation on the three facets of diversity for these Hill orders. For each diversity facet and Hill order, we fitted five models that included all possible combinations of effects: 1) no effect (i.e., null model), 2) only elevation, 3) only transect, 4) similar effect of elevation across transects (i.e., without interaction between elevation and transect or parallel slopes model), and 5) differential effect of elevation among transects (i.e., with interaction between elevation and transect or different slopes model). We used a negative binomial error distribution (and a log-link function) to account for overdispersion. Alternative models were compared based on the corrected Akaike information criterion (AICc) and we considered the models within two AICc units of the model with the lowest AICc value. Among these models, we selected the most complex and calculated the explained deviance (D²) as a measure of goodness of fit for GLMs. Analyses were conducted using 'MASS' (Venables and Ripley, 2002; v.7.3–51.6) and 'MuMIn' (Barton, 2013; v.1.43.17) R packages.

All the R Statistical Software analyses were conducted using v.4.0.2 (R Core Team, 2021).

3. Results

Pearson's correlations coefficients between number of trees per plot and any of the three facets of diversity (taxonomic, functional, and phylogenetic) for different Hill orders was very low, ranging between −0.22 and 0.08 and hence there was no need to consider the effect of sample size on diversity metrics.

When all transects were analysed together, taxonomic, and functional diversity were always significantly and positively correlated for q = 0, q = 1, and q = 2 (Fig. 2 and Table S4). However, taxonomic, and phylogenetic diversity, and functional and phylogenetic diversity were significantly and positively correlated only for q = 0 and q = 1. After analysing the transects independently, we found the following.

• Taxonomic vs. functional diversity: always significantly and positively correlated for q = 0, q = 1, and q = 2.

• Taxonomic vs. phylogenetic diversity: significantly and positively correlated for q = 0 and q = 1 (except in Podocarpus for q = 1).

• Functional and phylogenetic diversity: significantly and positively correlated for q = 0 and q = 1 (except in Podocarpus for q = 1).

Overall, results were consistent between transects. The three types of diversity values decreased from lower to higher elevations for q = 0, q = 1, and q = 2. This trend was strongest for q = 0 and weakest for q = 2, as indicated by both the slope of the regression curve (Fig. 3) and the explained deviance (Table 2). The only exception was phylogenetic diversity for q = 2, which appeared to increase slightly with the elevation. The best fitted models always included both elevation and transect (with or without interaction) as predictors, except for phylogenetic diversity for q = 2 (Fig. 3; Tables 2 and S5).

4. Discussion 4.1. Relationships between the three facets of diversity along elevation

Taxonomic, functional, and phylogenetic diversity responded to elevation in a similar manner (Fig. 3), and they were significantly positively correlated in most cases (Fig. 2 and Table S4). These positive correlations may be expected due to changes in species richness along elevational gradients, where more species are likely to exhibit more functional strategies and belong to more phylogenetic lineages (Losos, 2008). However, this is not always necessarily true and empirical studies often show that two communities with equal number of species can differ greatly in terms of number of strategies and/or lineages (Devictor et al., 2010; Safi et al., 2011) and that a community exhibiting similar strategies may foster distantly related lineages (trait convergence) (Freckleton and Jetz, 2009). Few studies have quantified the three types of diversity simultaneously along elevational gradients for plant communities (e.g., Tanaka and Sato, 2015; Chun and Lee, 2017). In the tropical Andes, this type of study has only been conducted for vertebrates (e.g., Dehling et al., 2014 for birds; Cisneros et al., 2014 for bats; and Dreiss et al., 2015 for rodents). Congruence between the three types of diversity along elevational gradients has not been reported for plants, although pairwise correlations between two facets of diversity have been found: e.g., functional and phylogenetic diversity (Tanaka and Sato, 2015; Chun and Lee, 2017). Regardless of the existence or absence of correlations between the different facets of diversity, all three provide complementary insights into the mechanisms that drive biodiversity patterns –as will be discussed in the following–, which may be obscured if only one is considered (Forest et al., 2007; Swenson, 2013).

4.2. Taxonomic diversity

Our results indicated a monotonic decrease in taxonomic diversity along elevation across all study transects that was consistent for different Hill orders (Fig. 3a), thereby suggesting that this decrease occurred regardless of whether taxa were weighted based on their abundance (q = 1 and q = 2) or not (q = 0).

General elevational decreases in taxonomic diversity of woody plants similar to those found in the present study have been widely reported (Kessler, 2001; Sandoya et al., 2017; Nottingham et al., 2018), being the elevational decrease in temperature the ultimate cause of this trend (Grytnes and McCain, 2007). However, studies have also detected a non-monotonic change in diversity along elevation, specifically a hump-shaped pattern where richness peaks at middle elevations (Gentry, 1995; Rahbek, 1995; Girardin et al., 2014; Salazar et al., 2015). This hump-shaped pattern in taxonomic diversity with elevation can be explained by the relatively stable condensation zone (cloud belt) that characterizes the middle elevations of tropical elevational gradients, which provides favourable conditions for most organisms despite the lower temperatures relative to premontane elevations (Kessler and Kluge, 2008; Huaraca Huasco et al., 2014). If this hump was present along our transects, our sampling design would have been unable to detect it because ⅰ) we did not sample below 1000 m a.s.l. (except for a very few plots in the Río Abiseo transect that did show an indication of a hump), and ⅱ) there was a gap in our sampling around 1500 m a.s.l., (where the highest diversity could be expected). These sampling limitations may have concealed a potential hump in species richness (Rahbek, 2005; Nogués-Bravo et al., 2008).

4.3. Functional diversity

Our results indicated a monotonic decrease in functional diversity along elevation across all transects that was consistent for different Hill orders (Fig. 3b). Theory predicts that the range of successful functional strategies narrows along an environmental stress gradient (Weiher and Keddy, 1995; Cornwell and Ackerly, 2009). In the case of elevation, the upslope increasingly harsher and more restrictive conditions –mainly in terms of temperature and resource availability (Stevens, 1992; Sundqvist et al., 2013)– is expected to lead to the selective survival of taxa with conservative strategies –promoting storage and defence– that enable them to cope with the environment, and to the cull of these with acquisitive ones –enhancing photosynthesis and growth– (Poorter et al., 2009; Machac et al., 2011). This pattern has been detected previously in TMFs both for individual functional traits –including leaf size decreases (Schneider et al., 2003; Salinas et al., 2011; Llerena-Zambrano et al., 2021) and wood density increases (Swenson and Enquist, 2007; Slik et al., 2010)–, and for overall community functional diversity (Duivenvoorden and Cuello, 2012; Schellenberger Costa et al., 2017). This evidence has been used in support of environmental filtering as an essential driver of plant community assembly (Swenson et al., 2012; Bañares-de-Dios et al., 2020).

The decrease in functional diversity was less pronounced for Hill numbers of higher order. Our results showed that Hill numbers that gave more weight to the most common entities (q = 2) accounted for two times less diversity than when the entities were strictly weighted by their abundance (q = 1) and four times less than when all entities received the same weight (q = 0) (Fig. 3b). These findings reflect how dominant/abundant species tend to be located at the core of the functional hyperspace whereas rare species tend to be peripheral. Dominant/abundant species exhibit strong traits–environment relationships enabling them a good performance in a given environment, thereby leading to greater demographic success. Conversely, rare species display sub-optimal traits –these species are more reliant on occasional migration events or stochastic environmental fluctuations–, which may explain their poorer demographic success (Umaña et al., 2015, 2017). In highly diverse ecosystems, such as tropical rainforests or coral reefs, rare species exhibit the less common traits (functional rarity), and thus disproportionately contribute to functional diversity, being irreplaceable and their loss not susceptible to being cushioned by the high species richness of the ecosystem as usually assumed. The reason for this is that the most common species have been shown to only support a relatively limited functionality (functional redundancy) (Mouillot et al., 2013, 2014; Leitão et al., 2016). Indeed, functional rarity sustains key ecosystem processes and services and provides an invaluable buffer against climate change and human disturbance (Cadotte et al., 2011; Violle et al., 2017). Therefore, the justification for conserving rare taxa should extend far beyond their taxonomic uniqueness, past, charisma or cultural significance (Tucker et al., 2012).

4.4. Phylogenetic diversity

Our results indicate that phylogenetic diversity decreased along elevation for q = 0 and q = 1 but increased slightly for q = 2 (Fig. 3c). In general, the overall decrease agrees with the predictions of the Tropical Niche Conservatism (TNC) hypothesis, which suggests that biological communities in the highlands are phylogenetically a subset of the communities found in the lowlands, and thus a decrease in phylogenetic diversity is expected with elevation (Wiens and Graham, 2005; Donoghue, 2008). This pattern has been observed for other biota in Andean TMFs (e.g., birds (Graham et al., 2009) and moths (Brehm et al., 2013)), but previous studies of Andean woody plants seem to give more support to the Out of the Tropical Lowlands (OTL) hypothesis (González-Caro et al., 2014, 2020; Tiede et al., 2016; Ramírez et al., 2019). Under OTL, the presence of temperate extratropical lineages within high elevation communities would be responsible for an increased phylogenetic alpha-diversity and turnover (Qian and Ricklefs, 2016), which is the pattern we found when considering phylogenetic diversity for the most abundant species (q = 2).

Therefore, our results are not in full agreement with either hypothesis and endorse neither the support nor the rejection of the TNC or OTL (see also Tolmos et al., 2022). In this sense an alternative intermediate hypothesis termed Environmental Crossroads (Neves et al., 2020) has been proposed, which posits that floras from different origins –tropical and temperate– co-exist at middle elevations, at the edges –lower and upper, respectively– of their thermal tolerances (Griffiths et al., 2021). In any case, further research is needed, including determining the phylogenetic turnover and nestedness, as well as biogeographical origin of clades (Linan et al., 2021) and changes in their age along the elevational gradient (Qian, 2014; Griffiths et al., 2020). A sampling schema that replicates transects along a broader latitudinal range, such as that used in the present study, would also help to assess the relative importance of the TNC vs. OTL. For example, several extratropical lineages such as Quercus prevail at high elevations in Colombia, which may have partly skewed the results reported by González-Caro et al. (2014, 2020) or Ramírez et al. (2019) based on studies that were latitudinally restricted to the northern Andes and which found evidence in agreement of OTL. Quercus is a Holarctic genus that migrated to South America in the last ca. 0.5 My (Hooghiemstra and Van Der Hammen, 2004) and there may have not been sufficient time for it to reach the central Andes, where niche conservatism and the upslope migration of clades originated in tropical lowlands could have been the predominant mechanisms responsible for population by Andean flora.

In general, extratropical lineages that have migrated from temperate to tropical regions along the Andes, such as Quercus or Alnus in Colombia or Araucaria and Nothofagus in the southern Andes, are better adapted to cold conditions than strictly tropical lineages, and thus they are more likely to thrive in the uppermost TMF communities (Segovia et al., 2020). These taxa were absent from the transects that we studied but we found other floristic elements of Laurasian (e.g., Prunus and Cornus) and Gondwanan (e.g., Podocarpus and Retrophyllum) temperate origin, some of which tended to be dominant locally in the past and nowadays are far less abundant, with the only exception of certain fragments of relict forests (see Vicuña-Miñano, 2005; Yaguana et al., 2012; Huamantupa-Chuquimaco et al., 2017). In this regard, the fact that our transects were located on isolated, protected areas difficult to access may have favoured the preservation of the native flora, allowing us to detect a slight increase in phylogenetic diversity at the highest elevations when giving more weigh to the most abundant taxa (q = 2) (Fig. 3c). Indeed, this slight increase resembles the phylogenetic pattern expected under the OTL. Therefore, as one moves far away from the central towards the northern or southern Andes, it would not be surprising to find an increase in the positive relationship between phylogenetic diversity and elevation, which may be detected by simply considering the relative abundances of taxa (q = 1) or equally considering all taxa (q = 0) given the greater prevalence of extratropical lineages.

Finally, it is important to note that the metrics used in this study to calculate phylogenetic diversity are strongly correlated with the number of lineages in recent evolutionary time but have weaker correlations with the number of lineages deeper in the evolutionary history of an assemblage (Dexter et al., 2019). Thus, our results might differ with the use of different metrics or the inclusion of different phylogenetic depths (Bose et al., 2019).

5. Conclusions

We found an overall decrease in taxonomic, functional, and phylogenetic diversity with elevation in Andean tropical montane forests (TMFs) (the only exception was the phylogenetic diversity when we over-weighted dominant species). Our experimental design is particularly robust because we analysed four elevational gradients over a broad latitudinal range, thus avoiding the potential bias of other studies that focused on the much better known but geographically constrained northern Andean TMFs (Pérez-Escobar et al., 2022). Declines in taxonomic and functional diversity with elevation are consistent with the expected outcomes of an environmental filtering process in which temperature is the main driver of environmental variation. Differences in the steepness of patterns of functional diversity between different Hill orders suggest that rare species play key and unique roles in the ecosystem, thereby supporting the conservation of rare taxa for reasons that extend far beyond their taxonomic or historical singularity. In addition, these findings suggest that niche conservatism and the upslope migration of tropical clades may be the predominant mechanisms responsible for shaping the flora of the central Andean highlands, given the limited arrival of temperate taxa in comparison to northern Andes. Further research across the entire Andean range, by integrating analyses of phylogenetic turnover and nestedness, as well as the biogeographical origin of clades and changes in their age along the elevational gradient will be important for elucidating the biogeographic history of the Andean slopes. The fact that the three types of diversity were positively correlated demonstrates that preserving TMF areas that harbour high taxonomic diversity will also preserve high functional and phylogenetic diversity in one of the most threatened biodiversity hotspots on Earth, which is still poorly understood compared to tropical lowlands (ForestPlots.net 2021) and threatened by anthropogenic intervention and global warming (Mathez-Stiefel et al., 2017; Feeley et al., 2020; Salinas et al., 2021).

Acknowledgements

Guillermo Bañares was funded through grants from the Spanish Ministry of Education (FPU14/05303), Escuela Internacional de Doctorado - Universidad Rey Juan Carlos (Doctor Internacional 2017) and the Education and Research Department of Madrid Autonomous Region Government (REMEDINAL TE; S2018/EMT-4338). The study was supported through three grants from the Spanish Ministries of Economy and Competitiveness and Science and Technology (CGL2013-45634-P, CGL2016-75414-P, and PID2019-105064 GB-I00), and a grant from Centro de Estudios de América Latina (CEAL) at Universidad Autónoma de Madrid and Banco Santander. We are indebted to those who helped with fieldwork: Stalin Japón, Wilson Remacha, Percy Malqui, "Rosho" Tamayo, Reinerio Ishuiza, Manuel Marca, Carlos Salas, José Sánchez, Anselmo Vergaray, Gonzalo Bañares, Ángel Delso, Julia González de Aledo, Mara Paneghel, Melecio Sullca, Raúl Huasurco, Honorato Pinto, Santiago Mamani, Aníbal Mamani, Armando Torrez, Justo Salas, Juan Mamani, Beto Apaza Coronel, Marcelo Reguerín, Óscar Quiroga, Juan Carlos Mamani, Rogelio Mamani, Vladimir Chura, Martín Chacón, and many others. We are especially thankful to Alex Nina, Jorge Armijos, Pablo Soliz, Leslie Cayola, Alfredo Fuentes, Peter Jørgensen and Tolentino Cueva, and the people of Los Alisos (Pataz, Peru) for their invaluable help. In addition, we are very grateful to Ministerio del Ambiente (MAE) in Ecuador, Servicio Nacional de Áreas Naturales Protegidas (SERNANP) in Peru, particularly Hugo Macedo, Vladimir Ramírez, Octavio Pecho, and Jhonny Ramos, and Ministerio de Medio Ambiente y Agua (MMAyA) in Bolivia. We extend our thanks to all the national park rangers who helped us, especially Rafael Galán, Guillermo Aguilar, Tito Heras, Percy Franco, Berardo Rojas, and Grover Benites. Permits to work in protected areas were granted by national authorities: Ecuador (MAE-DNB-CM-2015-0016; N° 001-2019-IC-FLO-FAU-DPAZCH-UPN-VS/MA), Peru (001-2016-SERNANP-PNRA-JEF), and Bolivia (MDRAyMA-VBRFMA- DGBAP-UAPVS N°2869/08).

CRediT authorship contribution statement

Guillermo Bañares-de-Dios: Writing - review & editing, Writing - original draft, Visualization, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. Manuel J. Macía: Writing - review & editing, Writing - original draft, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Data curation, Conceptualization. Gabriel Arellano: Writing - review & editing, Writing - original draft, Validation, Investigation, Funding acquisition, Data curation, Conceptualization. Íñigo Granzow-de la Cerda: Writing - review & editing, Writing - original draft, Validation, Investigation, Data curation, Conceptualization. Julia Vega-Álvarez: Writing - review & editing, Visualization, Formal analysis. Itziar Arnelas: Writing - review & editing, Validation, Investigation, Data curation. Carlos I. Espinosa: Writing - review & editing, Project administration. Norma Salinas: Writing - review & editing, Project administration. Luis Cayuela: Writing - review & editing, Writing - original draft, Visualization, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

Declaration of competing interest

Authors declare no conflict of interests.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2024.03.005.