Environmental determinants and spatial mismatch of mammal diversity measures in Colombia

Environmental determinants and spatial mismatch of mammal diversity measures in Colombia.— Including complementary diversity measures into ecological and conservation studies should improve our ability to link species assemblages to ecosystems. Recent measures such as phylogenetic and functional diversity have furthered our understanding of assemblage patterns of ecosystems and species, allowing improved inference of ecosystem function and conservation. We evaluated spatial patterns of taxonomic, phylogenetic and functional diversity of mammals in Colombia and identified their main environmental determinants, as well as interrelationships and spatial mismatch between the three measures. We found significant effects of elevation and precipitation on species richness, slope and species richness on phylogenetic diversity, and slope and phylogenetic diversity on functional diversity. We also identified a spatial mismatch of the three measures in some areas of the country: 12% of the country for species richness and 14% for phylogenetic and functional diversity. Our results highlight the importance of including species relationships within environmental drivers with biogeographical and distribution analyses and could facilitate selection of priority areas for conservation, especially when mismatch occurs between measures.


Introduction
Measuring biological diversity has become a major research challenge in the increasingly interdisciplinary fields of ecology and conservation biology (Cardillo et al., 2008;DeFries et al., 2010).Traditional measures have emphasized species richness and evenness (Purvis & Hector, 2000), but more sophisticated measures including evolutionary history and species function within ecosystems have recently been used (Buckley et al., 2010;Cadotte et al., 2011).These new measures purportedly allow more precise and comprehensive assessments of ecological and conservation issues related to biodiversity across spatial and temporal scales (Rosenzweig, 1995;Hadly & Maurer, 2001;Cadotte et al., 2011).However, interpreting the results of phylogenetic and functional diversity remains a challenge, particularly in understanding mechanisms potentially responsible for observed patterns (Calba et al., 2014;González-Maya et al., 2016).This situation hinders what could otherwise be considered a substantive advance in comprehensive conservation planning (Devictor et al., 2010;Dalerum, 2013;Dehling et al., 2014).
Phylogenetic diversity is a measure of evolutionary history, such as the diversity within an assemblage (Collen et al., 2011).It provides insights into the evolutionary distinctiveness of an assemblage and can be used to estimate the evolutionary potential for future ecosystem performance (Dalerum, 2013).In contrast, functional diversity incorporates elements of biodiversity that represent how ecosystems function (Dı áz & Cabido, 2001) and is important for understanding current ecosystem dynamics, resilience and services (Tilman et al., 1997;Díaz et al., 2013).Continued environmental degradation places growing pressure on the world´s biodiversity (Flynn et al., 2009), leading some authors to suggest we are facing the sixth mass extinction (Barnosky et al., 2011;Pimm et al., 2014;Ceballos et al., 2015).As the complexities of ecosystems and species interactions are unraveled, comprehensive and interdisciplinary conservation schemes are needed to account for the dynamic and interwoven patterns in the tapestry of the world´s natural capital (Dalerum et al., 2009;Devictor et al., 2010;Dalerum, 2013;Zupan et al., 2014).Incorporating multiple measures of species diversity including species richness, but more importantly evolutionary history and ecosystems functioning into planning, would be a major advance in conservation (Dalerum, 2013).
The Neotropical region has high biodiversity (Mora et al., 2011), but is one of the least biologically known regions of the world (Myers et al., 2000;Cardillo et al., 2006).This knowledge gap is a constraint for comprehensive conservation planning (Boitani et al., 2011;Visconti et al., 2011), exacerbated by our limited knowledge of how environmental drivers influence species diversity (Kerr, 1997;Colwell et al., 2008;Rondinini et al., 2011b;González-Maya et al., 2016).Most previous assessments of priority setting for the Neotropical region have focused on classic measures such as species richness (Gerardo Ceballos, 2007;Jenkins & Giri, 2008;Forero-Medina & Joppa, 2010;González-Maya et al., 2015), but few analyses have explored multi-measure approaches to account for evolutionary history and ecosystem function, and most analyses have been at global scales (Safi et al., 2011).
As a first step to further our understanding of biodiversity patterns, we provide the first nationwide assessment of terrestrial mammal diversity in the mega-diverse country of Colombia using a new generation of measurement tools.Specifically, we assessed i) spatial patterns of species richness, functional and phylogenetic diversity, ii) the influence of environmental factors on these patterns, and iii) the spatial mismatch of the three measures.

Material and methods
We first developed a grid of 127 points located in the centroid of 1 x 1 degree cells over the country, defining these cell centroids as our sampling unit.We selected this coarse resolution to match the scale of global distribution maps for mammal species developed for the IUCN Red List of Threatened Species (IUCN, 2012;Schipper et al., 2008) since the country lacks information for finer resolution.As these species distribution maps were developed for the global scale, we used this spatial resolution to reduce potential bias from overestimated distribution polygons and with illustrative purposes (González-Maya et al., 2015); nevertheless, we acknowledge potential limitations of using such a coarse resolution (Rondinini et al., 2011a).Therefore, our approach is based on centroids, namely sampling points, and regardless of the cell size, the data come from specific localities not related with resolution.Cells are used only for illustrative purposes, since points (e.g., random) do not allow to illustrate spatial patterns.Based on these distribution polygons, we extracted those species present or potentially present in the country (González-Maya et al., 2012, 2015, 2016).We calculated the centroid of each cell (i.e., cell geographic center) and extracted all species overlapping each centroid as a proxy of the species assemblage present in each cell (Safi et al., 2011).Given there is no other detailed information for the country regarding species distributions (Solari et al., 2013) and as these polygons were corrected by national experts (Schipper et al., 2008), we believe this was the best available approach for country-wide mammal analyses (González-Maya et al., 2015, 2016).For cells overlapping the edge of the country, we estimated the centroid after clipping the perimeter of the country over the grid.
For each cell assemblage we calculated species richness and phylogenetic and functional diversity.phylogenetic diversity (PD) was estimated using the Meredith et al. (2011) phylogenetic tree and Faiths phylogenetic diversity index (Faith, 1992), which is based on the sum of the branches connecting all species in phylogenetic space for a given species assemblage.Functional diversity (FD) was calculated using the Petchey and Gaston FD index (Petchey & Gaston, 2002b;Safi et al., 2011;González-Maya et al., 2016), where the FD index is defined as the sum of the branches necessary to connect all species in trait space.To do this, we first estimated a matrix of distances based on the Gower distance (because we used both qualitative and quantitative traits) and built a dendrogram with this distance matrix.We then summed the distance of all branches within the tree for each assemblage (Petchey & Gaston, 2002b).High FD values indicate high complementarity of species functions, therefore low redundancy based on the traits used and low values indicate lower diversity of species functions, thus higher redundancy (Safi et al., 2011).We used Pearson's correlation to assess the degree of lineal relationships of the three measures.Even although new measures have improved Faiths phylogenetic diversity index (Chao et al., 2015), we still used this approach because it is the most similar approach to FD and because our analyses is not based on sampling data.
All traits are available in other sources, and we note that not all species had complete data, especially those known only from a few localities or those with no distribution information available through the IUCN database.We therefore used approximations to trait values based on the closest species relative within the same genus.
To assess the effects of environmental variables on diversity measures, we estimated mean cell annual precipitation and temperature, elevation and slope.Precipitation and temperature (i.e., annual precipitation (BIO12) and annual mean temperature (BIO1)) were calculated from the WorldClim database (Hijmans et al., 2005) based on representative estimates from 1950-2000 at 30 arc-seconds, while elevation and slope were derived from the Hydro1k Digital Elevation Model for South America, also at 30 arc-seconds (US Geological Survey, 2012).We used ordinary least squares regression to assess the influence of environmental variables on species richness (González-Maya et al., 2016), the influence of the environmental drivers and species richness on phylogenetic diversity, and the influence of environmental variables, species richness and phylogenetic diversity on functional diversity.For each measure, we generated all possible variable combinations without interaction of model terms, including 15 models for species richness, 31 for phylogenetic diversity and 63 for functional diversity (supporting information).
We identified the best performing models using Akaike Information Criteria corrected for small samples (AICc) and Akaike weights (Wagenmakers & Farrell, 2004), selecting those models with Δ > 2 as significantly best-performing models (Burnham & Anderson, 1998;Burnham et al., 2011).We used adjusted R 2 values to estimate the proportion of variation in the dependent variable explained by the model variables.After selecting the best models according to AICc, we identified the variables of the best performing models and we calculated the variable coefficients and the Variance Inflation Factor (VIF) to assess correlation among them; models containing variables with scores > 7.5 were considered correlated (O'Brien, 2007).As we ran all main effect model combinations, when this occurred we discarded those models and selected the model with the next lowest AICc value (O'Brien, 2007).
We estimated the Koenker studentized Breusch-Pagan statistic, K(BP), to assess the reliability of standard errors when heteroscedasticity was present.If the K(BP) was significant, we used the robust probability instead of the raw probability estimation (Breusch & Pagan, 1979;Koenker, 1981).Heteroscedasticity and non-stationarity indicate that the relationship between the dependent variable and the drivers would not change with changes in the magnitude of the drivers, and the relationship is not equal across geographic space, respectively.Moran´s I tests were used to test for residual clustering; clustered residuals indicate some variables or terms were missing from the model (Li et al., 2007).
To explore where spatial mismatching occurred and where selected models did not perform adequately (i.e., indicating at least one important variable was missing from the model), we performed a hot-spots analyses using the residuals of the selected models based on the Getis-Ord Gi* statistic, estimating z-scores and p-values for each cell (Getis & Ord, 1992;Ord & Getis, 1995).Hot spots are those where z-cores are significant (p > 0.05), therefore indicating where high clustering occurs.P-values are considered significant when z-scores estimated with the cell and its neighbors differ from expected when compared proportionally to the sum of all features.Significant z-scores in a cell (p > 0.05) indicates clustered residual patterns and suggests one or more explanatory variables are missing in the model for that cell, in turn indicating spatial mismatch of the measures, overall model and the explanatory variables (Getis & Ord, 1992;Ord & Getis, 1995).We mapped cells with high levels of residual clustering that were significant (p > 0.05), therefore the spots of significant spatial mismatch of the terms of the model.All geographic and statistical analyses were performed using ArcGIS 10.2 (Environmental Systems Research Institute, 2013) and Spatrial Analyst extensions, and R software (R Team Development Core, 2008) and VIF package.

Results
We found a heterogeneous distribution of the three diversity measures of mammals in Colombia.We observed that species richness was concentrated near the center of the country, specifically towards the Andes piedmont, while the lowest values were concentrated in the Northern Llanos region (Eastern Colombia) and the northernmost portion of the country (i.e., Guajira Peninsula; fig.1A).Phylogenetic and functional diversity were similarly distributed with greater values towards the southern Andes (FD) and the Eastern and Western Andes cordilleras in the north (FD); lowest values were found for Eastern Colombia (i.e., Llanos) and an apparent homogenous distribution of both PD and FD was also found for Eastern Colombia, with homogeneous PD for the Llanos and homogenenous FD for the Amazon regions (fig.1B, 1C).The three measures were highly related; species richness and phylogenetic diversity were highly correlated (Pearson = 0.96, p < 0.001), as were phylogenetic and functional diversity (Pearson = 0.95, p < 0.001), and species richness and functional diversity (Pearson = 0.90, p < 0.001).
For species richness, the best model (table 1) included elevation and precipitation as the most important influencing variables.Nevertheless, Moran´s I test was positive for clustering of the residuals (Moran´s Index = 0.37, p < 0.001) and the overall proportion of variation in the dependent variable explained was low (R 2 = 0.45).As the K(BP) indicated non-stationarity of the model (K[BP] = 2.36, p = 0.0381), we used robust probabilities (table 1).For phylogenetic diversity, the best model included species richness and slope, with a high proportion of the variability explained (adjR 2 = 0.94).Moran´s I (Moran´s Index = 0.005, p = 0.15) and K(BP) tests were not significant, indicating the change in relationship between the dependent variable and the determinants would not change when the magnitude of determinants change and that this relationship was constant across geographic space (table 1).Stationarity and heteroscedasticity were also found for the selected model.For functional diversity, the best model included phylogenetic diversity and slope, with a high proportion of variability explained by the model (adjR 2 = 0.94), and indicated non-clustering of the residuals (Moran´s Index = -0.032,p = 0.09; table 1).The variable slope in both cases had less influence on phylogenetic and functional diversity than did species richness and phylogenetic diversity (i.e., species richness explained 61.56% of PD and PD explained 62.66% of functional diversity).
Spatial mismatch for the three models identified areas where the models failed to explain the different diversity measures.For species richness, model mismatches occurred in the southern Andean region     and Guajira peninsula in the northernmost part of the country, covering ~12% of the country (fig.2A).A clear mismatch between environmental determinants and phylogenetic diversity occurred along the southern Pacific coast and to a lesser extent in a portion of the Amazon basin (~14%; fig.2B).Greatest clustering and mismatch for functional diversity occurred in the Sierra Nevada de Santa Marta, Paramillo Complex and Guajira peninsula of the Caribbean region (~14%; fig.2C).

Discussion
We provide the first quantification of phylogenetic and functional diversity of terrestrial mammals in Colombia.Our results show a high degree of relatedness between biodiversity measures, including high similarity in spatial patterns, similar to large scale analyses in other latitudes (Pavoine & Bonsall, 2011;Barnagaud et al., 2014;González-Maya et al., 2016).The evaluation of these complementary measures demonstrates how different attributes of biodiversity can be used simultaneously to better understand the assembly of species communities (Pavoine & Bonsall, 2011).Furthermore, previous works have identified how both phylogenetic and functional characteristics can explain species co-occurrence at regional and global scales (Barnagaud et al., 2014).Our results highlight that even when the taxonomic, phylogenetic and functional diversity significantly overlap, some areas will differ between the three measures (i.e., spatial mismatch), suggesting other factors affect assemblage composition and species co-occurrence.
The three diversity measures were strongly correlated, supporting how taxonomic richness drives phylogenetic diversity (Safi et al., 2011;Dalerum, 2013), which in turn is a key driver of functional diversity (Safi et al., 2011).Assemblages with larger or shorter evolutionary histories can markedly affect functional diversity and stability of those assemblages (Safi et al., 2011).However, for global analyses of these measures on mammals (Safi et al., 2011), spatial mismatches occur at different scales and geographic localities where additional factors or distinctive evolutionary histories have an unidentified effect on both measures (Devictor et al., 2010;Zupan et al., 2014).This is clear for areas previously found to have high endemisms and singularity like Sierra Nevada de Santa Marta and Southern Andes (Le Saout et al., 2013;Solari et al., 2013).Our results further indicate that selecting the areas with species-rich assemblages for conservation priorities will typically directly affect other biodiversity measures (Dalerum, 2013), and that areas identified as mismatches require further assessment and stronger consideration for conservation.
Environmental determinants, or environmental filtering, determines a substantial proportion of the basic diversity measure, namely species richness (Messier et al., 2010;Pavoine & Bonsall, 2011), while the effect of these variables on other measures was less influential in our study than expected (González-Maya et al., 2016).Previous studies of trait-and phylogenetic-based diversity measures have proposed environmental filtering as the most likely driver of these measures (Messier et al., 2010;Pavoine & Bonsall, 2011;Safi et al., 2011;Swenson et al., 2012).However, in our study these drivers were not sufficient to fully explain variation in species richness.Safi et al. (2011) evaluated these three measures of diversity at a global scale and found they have a significant degree of surrogacy.
Nevertheless, there appears to be considerable spatial mismatch across geographic scales, suggesting greater support for local and environmental factors influencing species assembly at taxonomic, functional and evolutionary levels.Local drivers, such as assemblage time, environmental constraints and biogeographic scale, have been proposed to explain phylogenetic and functional diversity patterns (Kraft & Ackerly, 2010;Spasojevic et al., 2014), which in our study seems to be a key aspect.Environmental constraints significantly affect species richness, slope, as a locally-defined variable affects functional and phylogenetic diversity, and assemblage evolutionary distinctiveness affects functional diversity.Our identification of slope as an important determinant of phylogenetic and functional diversity has been previously identified for substantially different systems and taxonomic groups (Cardoso et al., 2011;Punchi-Manage et al., 2013).In this study, as a locally-constrained variable, it is likely representative of locally-defined conditions of species assemblage.
Our results can certainly be refined since both phylogenetic and functional diversity measures are highly influenced by the phylogeny and traits used, respectively, and the method used to estimate both measures (Dalerum, 2013).Our selection of both measures was based on using similar dendrogrambased measures, so the comparison and simultaneous evaluation of both measures was congruent (Faith, 1992;Petchey & Gaston, 2002b;Dalerum, 2013).As previously stated for concepts such as redundancy (Naheem, 1998;Bueno et al., 2013;Mouillot et al., 2013), as more information becomes available and we further understand evolutionary and ecological aspects of a group, analyses on diversity measures can be further refined.Nevertheless, as the first effort to map phylogenetic and functional diversity for Colombia, our results could be incorporated into conservation schemes and facilitate further exploration of these topics not only in Colombia but also in other Neotropical countries.Furthermore, we acknowledge that one degree cell size in a topographically complicated region such as the Andes surely average environmental diversity data, blurring the fine-grained ecological mechanisms that might be responsible for the observed patterns.Thus, new studies along elevational or aridity gradients are needed in this area, probably one of the world's richest in mammal species.
Including complementary diversity measures helps us to further understand species assemblages by incorporating not only species richness but also species evolutionary history, function and patterns (Cadotte et al., 2011;Belmaker & Jetz, 2013;Monnet et al., 2014).Furthermore, it allows to explore mechanisms linking species to environment, ecosystem processes and vulnerability (Petchey & Gaston, 2002a, 2002b;Biswas & Mallik, 2010;Kraft & Ackerly, 2010;González-Maya et al., 2016), thus providing valuable information for more effective conservation planning (Barnagaud et al., 2014).Our results help to understand diversity patterns in Colombia, and our preliminary mapping could help define priorities for complementarity and singularity (Devictor et al., 2010;Zupan et al., 2014).Despite the potential surrogacy of diversity measures, identifying areas with high values of any of these would likely result in more comprehensive and integral conservation planning, integrating singularities at taxonomic, evolutionary and ecosystem function levels (González-Maya et al., 2016).The variation in three important biodiversity measures, and the fact that one measure cannot represent the entire reality of species assemblages, suggests species assemblages are best represented using multiple metrics.This approach could increase the likelihood that areas for conservation are adequately identified and provide a basis for integral conservation planning.

Fig. 2 .
Fig. 2. Spatial clustering of residuals graphed as standard deviations indicating spatial mismatch of models to explain the influence of environmental determinants on mammal species richness (A), phylogenetic diversity (B) and functional diversity (C) in Colombia.

Table 1 .
Best performing candidate models (M; s.Selected model) testing the influence for environmental drivers on mammal species richness, phylogenetic diversity and functional diversity in Colombia: R(SE).Robust standard error; R(P).Robust p-value; VIF.Variance Inflation Factor; AICc.Corrected Akaike Information Criterion; K(BP).Koenker's studentizedBreusch-Pagan Statistic and p-value.