INTRODUCTION

Species loss and biotic homogenization have been recognized as major consequences of human-induced disturbance on natural ecosystems (^{Gámez-Virués 2015}). Mammals are recognized as one of the faunal groups most sensible to habitat loss fragmentation, and alteration (^{Fahrig 2003}; ^{Michalski & Peres 2005}; ^{Dotta & Verdade 2011}). Conversion and alteration of original vegetation cover may change mammal assemblages’ composition and diversity, and species behavior, consequently altering the community structure (^{Fahrig 2003}; ^{Lyra-Jorge et al. 2008}; ^{Dotta & Verdade 2011}). These processes also act as a selective filter on mammal communities, which favors generalist species in detriment of specialist ones (^{Fahrig 2003}; ^{Dotta & Verdade 2011}). This ongoing species impoverishment would impact ecosystem properties and functioning, limiting the response-diversity and resilience of ecosystems against disturbances (^{Pringle et al. 2007}; ^{Fritz & Purvis 2010}; ^{Gámez-Virués 2015}).

The Middle Magdalena Valley in Colombia hosts a variety of unique ecosystems that harbor high levels of biological diversity and support a wide range of ecosystem services (^{Garzón & Gutiérrez 2013}). Even though this region has high levels of biodiversity, it has been categorized as a deforestation hotspot (^{Garzón & Gutiérrez 2013}; ^{Sanchez-Cuervo & Aide 2013}). Increasing land-use intensification processes (e.g., cattle, oil palm plantations, petroleum exploration and extraction industry, mining, and illegal crops) has led to the loss of more than 85% of its native vegetation, with further deforestation since 1920 (^{Etter et al. 2006}; ^{Garcia-Ulloa et al. 2012}; ^{Garzón & Gutiérrez 2013}). Besides, the reconfiguration of natural habitats promotes diverse pervasive processes, including edge effects, barrier effects, and anthropogenic disturbances (^{Saunders et al. 1991}; ^{Kupfer et al. 2006}; ^{Holland & Bennett 2009}). Considering these threats, it is critical to understand how wildlife persists in this human-dominated landscape to guide conservation actions.

Despite all modifications caused by human activity, forest patches in the Middle Magdalena Valley still retain high levels of mammal biodiversity (see ^{Boron & Payan 2013}; ^{Boron et al. 2016}; ^{Figel et al. 2015}, ²⁰¹⁹), including threatened species (e.g., Myrmecophaga tridactyla Linnaeus, 1758) and species performing critical ecological functions (e.g., top predators such as Panthera onca [Linnaeus, 1758] and Puma concolor [Linnaeus, 1771] and long- distance seed dispersers like Pecari tajacu [Linnaeus, 1758]). These patches also vary greatly in area, shape, isolation, and fragment dynamics, offering a unique opportunity to evaluate the effects of habitat loss and fragmentation on mammal diversity and their func tional roles. Understanding how species richness responds to environmental disturbance and how it relates to ecosystem functioning is critical to design effective conservation strategies, especially those aimed to restore or conserve healthy ecosystems (^{Cadotte et al. 2011}).

In human-dominated landscapes, species responses to forest fragmentation and inter-patch isolation may vary depending on surrounding matrix composition, configuration, and quality (^{Perfecto & Vandermeer 2010}; ^{Fahrig et al. 2011}). Responses to fragmentation can also vary among species due to variation in species-specific traits related to use of habitat, acquisition of resources, niche breadth, and movement (^{Antongiovanni & Metzger 2005}; ^{Holland & Bennett 2009}; ^{Lyra-Jorge et al. 2010}; ^{Brodie et al. 2015}). Landscapes with heterogeneous and high-quality matrices, such as those with secondary forests and diverse types of land covers, are expected to support more species as they may (i) provide resources for both forest specialist and habitat generalist species, and (ii) exhibit enhanced landscape connectivity (^{Antongiovanni & Metzger 2005}; ^{Lyra-Jorge et al. 2008}; ^{Perfecto & Vandermeer 2010}; ^{Fahrig et al. 2011}). Consequently, effective conservation in fragmented landscapes requires a deep understanding of how species respond to landscape changes.

Here, we evaluated whether the spatial configuration of forest patches predicts the species richness and functional diversity of medium and large-bodied mammals in a fragmented landscape in the Middle Magdalena Valley in Santander department, Colombia. We predict that species richness, functional diversity, and the number of functional groups are increased in less-fragmented and well-connected forest patches, irrespectively of their size and shape, as they can be colonized by both specialist and generalist species from adjacent areas and may provide more opportunities for immigrating species to become established. On the contrary, they will be reduced in highly fragmented and isolated patches probably due to habitat filtering where some species (e.g., specialists) do not have the traits that allow them to arrive and/or survive in such patches. We also discuss some management recommendations to improve the conservation value of this highly fragmented, human-dominated, unprotected landscape.

MATERIALS AND METHODS

Study area

The landscape of study is located in the Middle Magdalena Valley, Santander department, Colombia (central coordinate: 06°59’41.6" N, 73°44’51.5" W, elevation range between 70 and 130 m). The area is classified as Tropical Moist Forest and is characterized by a bimodal rainfall regime with an average annual rainfall of 2 917 mm, a mean annual temperature of 27.9°C, and relative humidity of 80% (^{IDEAM 2019}). Forest patches in the study area are surrounded by a matrix of contrasting characteristics in terms of composition and structure. For this study, we sampled eight forest patches that differ in a variety of ways, including size (17.6 to 141.2 ha), shape, isolation, and fragmentation degree (Fig. 1). It was not possible to select these forest patches randomly owing to insurmountable access difficulties (e.g., large-scale oil palm plantations did not allow us access to their lands).

Field surveys

We conducted field surveys for medium-and large-sized mammals between June 2015 and May 2017, using two complementary non-invasive sampling methods: (1) active search and (2) camera-trapping. We selected these methods as they are appropriated to register species with a wide range of ecological and behavioral characteristics (^{Beca et al. 2017}).

We collected data through active search by walking on transects (e.g., natural and human-made trails) at an average speed of 1 km/h, searching for direct evidence of mammals (e.g., sightings, vocalizations) and signs of mammalian activity (e.g., tracks, burrows, feces, food leftovers). Transect length varied between 0.6 and 2.5 km depending on the size of each forest fragment. Surveys started in the morning (8:00 h) and at night (18:00 h), lasting for at least five hours per fragment, depending on the number of records and the environmental conditions. We carefully measured, photographed, and compared mammal signs with those from ^{Emmons & Feer (1997)}, ^{Aranda (2012)}, and ^{Suárez-Castro & Ramírez-Cháves (2015)}. Sightings, vocalizations, and signs for each species were considered as a single record for each sampling session. We conducted active search during 24 sampling campaigns, each consisting of eight days of fieldwork, for a total sampling effort of 192 days (i.e., 24 days per forest patch).

We also installed 16 camera traps (Bushnell Trophy Cam model 119435) in sites presumed to increase the probability of capturing wild mammals (targeted placement), including natural trails, near water bodies, or where signs of mammalian activity were found. The cameras were set up along with terrain conditions to heights between 50 and 70 cm above the ground and were programmed to take either 10- seconds videos (720 x 480) or three photos per trigger (8 mp). Cameras were baited with fruits or a mix of tuna and sardines to increase the likelihood of detection. We assume that images of the same species at the same station within a 1-h period correspond to the same individual. We carried out camera-trapping surveys over an 8-day period during the rainy season (November 2016), using two cameras per forest patches, for a total sampling effort of 128 camera trap days (i.e., 16 trap days per forest patch).

We recorded the specific locality information, geographic coordinates, and checked the taxonomic identification for all registered specimens. We followed the taxonomic nomenclature by ^{Wilson & Reeder (2005)} for taxonomically stable species. Besides, we followed ^{Rylands & Mittermeier (2009)} and ^{Ruiz-García et al. (2018)} for primates. Species weighing between 1 and 7 kg were categorized as medium-sized (^{Chiarello 2000}), while species weighing over 7 kg were categorized as large (^{Emmons & Feer 1997}). It is important to note that we included in our species list two relatively small-sized species that were recorded during the surveys (e.g., Caluromys lanatus and Metachirus nudicaudatus).

Landscape metrics

We estimated four landscape metrics for each study site, which correspond to the explanatory variables used in the regression models: patch size, patch shape, landscape connectivity, and forest fragmentation. Patch size and patch shape were estimated based on the land cover map for the study area using the Polygon Complexity tool in the software QGIS 3.8.3 (^{QGIS 2019}). Mean current density for each forest patch was estimated based on the multi-scale structural connectivity model generated by ^{Meza-Joya et al. (2019)}, using the zonal statistics tool in SDMtoolbox 2.0 (^{Brown et al. 2017}). This connectivity model combines several landscape features, including normalized difference vegetation index, surface temperature, land cover classes, fragmentation index, river proximity and road proximity (For additional details see ^{Meza-Joya et al. 2019}). Forest fragmentation in the vicinity of each sampling site was measured by forest area density (FAD), calculated as the percentage of forest pixels in a fixed-area neighborhood (^{Riitters et al. 2012}; ^{Riitters & Wickham 2012}). Because the most appropriate spatial scale for the local landscape was unknown, we evaluated five measurement scales defined by square neighborhood sizes equal to 7, 13, 27, 81, 243 pixels (see ^{Riitters & Wickham 2012}). We then estimated the average FAD value for each sampled forest patch, which ranged from 0 (highly fragmented) to 100% (not fragmented or intact).

Fig. 1 Study landscape located at the Middle Magdalena Valley, Colombia. The studied forest fragments are numbered. 

Data analysis

We first estimated mammal species richness for each forest patch as the number of species recorded in each site. As the performance of many functional diversity metrics is highly sensitive to undersampling (^{Cardoso et al. 2014}; ^{van Der Plas et al. 2017}), we assessed the completeness of our sampling by computing the sample-size-and coverage-based accumulation curves for each sampled patch. Because sample coverage was somewhat variable among forest patches (see Results), our raw estimates of species richness could be slightly biased by differences in sample completeness. Thus, following ^{Chao & Jost (2012)}, we estimated the species richness (q=0) using coverage-based extrapolations. The curves and estimates were performed using species frequency counts (i.e., number of specimens recorded per species) as a surrogate for species abundance in the R package iNEXT (^{Hsieh et al. 2016}).

We combined species abundance data (i.e., species frequency counts) from each sampled site and species trait information to calculate the following functional diversity indices: functional evenness (FEve) and functional dispersion (FDis). The former index (FEve) measures the distribution regularity of functional traits within a functional space for a specific community (^{Villéger et al. 2008}; ^{Farias & Jaksic 2009}), while the latter one (FDis) measures the dispersion in trait abundances in the multidimensional trait space of a community (^{Laliberté & Legendre 2010}). We compiled functional traits from previously published works (e.g., ^{Carvalho et al. 2010}; ^{Safi et al. 2011}; ^{Magioli et al. 2015}, ²⁰¹⁶) including physical, dietary, behavioral, and environmental sensitivity data (Table S1, for additional details see ^{Magioli et al. 2015}).

We calculated functional indices applying Gower dis- tances (^{Gower 1966}), which allows mixed trait types while giving them equal weight (^{Villéger et al. 2008}), using the Calinski criterion (^{Calinski & Harabasz 1974}) to estimate the optimal number of functional groups. For convenience, we constrained FDis values to a range between 0 and 1 dividing them by 10 (see ^{Laliberté & Legendre 2010}). We also calculated Rao’s square entropy (^{Rao 1982}; ^{R Core Team 2019}), but the results were highly correlated to those found for FDis (r=0.978) and are not presented. All functional indices were calculated with the function dbFD from the FD package (^{Laliberté & Legendre 2010}) in R version 3.6.1 (^{R Core Team 2019}).

We used beta regression models to evaluate if landscape metrics (i.e., patch size, patch shape, current density, and FAD) explain the variation in functional diversity (response variable). Beta regression is appropriate for dealing with continuous response variables that are restricted between zero and one, as is the case of FDis (^{Ferrari & Cribari-Neto 2004}). We chose FDis for this analysis because it is independent of species richness by construction (see ^{Laliberté & Legendre 2010}) so that its use ensures that the number of species does not influence the fit of the resulting models. To avoid multi-collinearity, we evaluated pairwise correlations among explanatory variables using Spearman’s rank correlation coefficients (r). We used log-transformed explanatory variables to improve linearity and reduce data dispersion. Model-averaging started with a global model with all the explanatory variables. Then, we used the dredge() function of the MuMIn package (^{Barton 2019}) to create a set of models with all possible variable combinations, including an intercept-only model, and from these we identified our best models based on comparisons of the second-order Akaike information criterion (AICc). Besides, we calculated the pseudo r2 of the best-ranked models (i.e., the five models with the lowest AICc) as a measure of the models fit.

In order to evaluate the response of species richness to landscape features, we constructed Generalized Linear Models (GLM) with a negative binomial error distribution using the same explanatory variables (i.e., patch size, patch shape, current density, and FAD) and extrapolated species richness as a response variable to avoid potential bias related to our sampling effort. For this, we followed the above procedures for environmental variable transformation and model averaging. Lastly, we estimated the r² of the best-ranked models as a measure of the models fit. We also run the analysis using the observed species richness as a response variable, but results were highly similar to those found using the extrapolated richness and are not presented.

RESULTS

We recorded 25 species of medium and large-bodied mammals, encompassing eight orders and 17 families (Table 1). Carnivores were the most representative order in the study area, with nine species belonging to four families: Canidae, Felidae, Mustelidae, and Procyonidae. Three species of rodents, three primates, three pilosa, three marsupials, two cingulate, one artiodactylid, and one lagomorph were also recorded. Species classified as habitat and diet gener- alists were more frequently recorded (see Table S1). Regarding their conservation status, two species are categorized as Near Threatened (Lontra longicaudis and Panthera onca) and one species (Aotus griseimembra) is listed as Vulnerable in the IUCN Red List (Table 1). Species richness was variable across the sampled forest patches (Table 1), oscillating between 12 and 21 species (mean ± SD=17.25 ± 2.60). Sample coverage (Fig. 2a) was also variable between sites, with values among 0.916 and 0.965 (mean ± SD=0.948 ± 0.02), suggesting that most species were sampled in each site. Similarly, sample-size-based curves (Fig. 2b) showed a little gain in estimates of species richness (from 2 to 4 species, mean ± SD= 2.915 ± 0.93) when doubled the sampling size.

The recorded mammal assemblages showed a wide range of functional diversity values in all the used metrics (Table 2). The number of functional groups per assemblage was variable, ranging from 7 to 9 (mean ± SD=8.13 ± 0.60). Values for FEve ranged from 0.481 to 0.667 (mean ± SD=0.557 ± 0.048), whereas FDis oscillated between 0.401 and 0.503 (mean ± SD=0.469 ± 0.029). There was a strong positive correlation between species richness and the number of functional groups (r=0.730). Species richness was moderately correlated to functional diversity metrics, both positively to FDis (r=0.582) and negatively to FEve (r=-0.561). There was a moderate negative correlation between FEve and FDis (r=-0.523).

Pairwise correlations among explanatory variables (i.e., patch size, patch shape, current density, and forest fragmentation) were between very low and moderate, ranging from -0.047 to 0.601. FDis was weakly positively correlated to patch size (r=0.429), patch shape (r=0.095), current density (r=0.190), and FAD (r=0.433). Species richness was weakly posi- tively correlated to patch size (r<0.422), moderate positively to FAD (r=0.602), and weakly but negatively to patch shape (r=-0.108) and current density (r=-0.024). Based on average modeling values of the group of proposed models, FAD was the most important parameter explaining species richness and functional diversity (measured as FDis) of medium-and large-sized mammals. For species richness, FAD was positively related to the number of species recorded (r² = 0.566, relative importance of overall predictor=0.87, Table 3, Fig. 3a). Similarly, FDis was also positively affected by FAD (pseudo-r² = 0.659, relative importance of overall predictor=0.96, Table 4, Fig. 3b).

DISCUSSION

We show that mammal species richness and functional diversity vary according to forest fragmenta tion in a human-dominated landscape. As predicted, we found that species richness and most functional diversity metrics were negatively related to forest fragmentation, but contrary to our expectations they were not related to landscape connectivity. Overall, our results show that mammal assemblages in less fragmented forested areas display a wider array of ecological traits and support a broader spectrum of ecological functions, but are less resilient to environmental perturbations and less resistant to ecological invasions. We discuss the insights arising from this work and consider the implications for management and conservation.

Table 1 Medium-and large-sized mammals recorded in Middle Magdalena Valley, Colombia, including the IUCN Red List category (IUCN) for each species (https://www.iucn.org/): Data Deficient (DD), Least Concern (LC), Near Threatened (NT), Vulnerable (VU). Record type: burrows (B), camera-trap (C), feces (Fc), food leftovers (Fl), tracks (T), sightings (S), and vocalizations (V). Forest patch identification as in Fig. 1. 

Fig. 2 Sample-coverage-based (A) and sample-size-based (B) curves (with 95% confidence intervals) for the studied medium-sized and large mammal assemblages at Middle Magdalena Valley, Colombia. Continuous line: interpolated data. Discontinuous line: extrapolated data. Forest patch identification as in Fig. 1. 

Fig. 3 Estimated relationships (with 95% confidence intervals) between forest fragmentation and species richness (A) and functional diversity (B) for medium- and large-sized mammal assemblages at Middle Magdalena Valley, Colombia. For the analysis in B, we use functional dispersion (FDis) as a measure of functional diversity. 

Table 2 Species richness (observed (OSR) and extrapolated (ESR)) and functional diversity metrics (Functional dispersion [FD], Functional evenness [FE] and functional groups [FG] for the studied medium-sized and large mammal assemblages at Middle Magdalena Valley, Colombia. Forest patch identification (FP) as in Fig. 1. 

Table 3 Top-ranked candidate GLM models explaining variation in medium-and large-sized mammal species richness. df: degrees of freedom. LLR: Log-likelihood ratio. F: fragmentation. C: connectivity. S: size. SH: shape. 

Table 4 Top-ranked candidate models explaining variation in medium-and large-sized mammal functional dispersion (FDis). df: degrees of freedom. LLR: Log-likelihood ratio. F: fragmentation. C: connectivity. S: size. SH: shape. 

Species richness

Previous studies in adjacent areas to the study region recorded 18 to medium- and large-bodied mammals (^{Boron & Payan 2013}). Species recorded in this study (n=25) represent about 17% of medium-sized and large mammals in Colombia, according to ^{Ramírez- Chaves et al. (2016)}. The difference between the observed and extrapolated species richness values obtained indicates that our sampling effort was sufficient to record most species in the sampled sites (Table 2), even though some species may be present but missed by our sampling effort. Mammal assemblages in the study area are characterized by a predominance of generalist species (see Table S1) categorized as Least Concern by the IUCN Red List (Table 1). Nonetheless, the study area still maintains sensitive species that exhibit habitat and feeding specializations, such as Aotus griseimembra, Lontra longicaudis, and Panthera onca, which may not have been able to adapt to habitat modifications but are forced to move through lower quality habitats when searching for resources. We did not record some sensitive species such as cervids (genus Mazama) and more susceptible species as Ateles hybridus Geoffroy, 1829 and Myrmecophaga tridactyla, which are probably extinct in the sampled forest patches or were missed during our sampling.

Our results indicate that less fragmented forest patches maintain a greater number of medium and large-bodied mammal species. According to the “habitat heterogeneity hypothesis” (^{MacArthur & MacArthur 1961}), structurally complex habitats, such as forest, provide more niches and complementary resources and, hence, increase species diversity. Our findings are consistent with previous mammal studies showing that complex environments provide greater opportunities for resource partitioning and, hence, are expected to support high levels of diversity (e.g., ^{Andrade-Núñez & Aide 2010}; ^{García-Herrera et al. 2015}; ^{Arévalo-Sandi et al. 2018}; ^{Sukma et al. 2019}). In this sense, anthropogenic disturbances in the study area may act as a relevant environmental filter operating to exclude species whit narrow physiological tolerance or stable habitat requirements (^{Gámez-Virués 2015}; ^{García-Llamas et al. 2019}). The loss of some large mammal species (e.g., Panthera onca and Pecari tajacu) from local assemblages in the study area will result in the loss of key functional groups and ecosystem functions (e.g., top predation and long-distance seed dispersal).

Functional diversity

Higher FDis values indicate a wider distribution of species in the functional space, with more individuals occupying the margins of the functional space in relation to its centroid (^{Laliberté & Legendre 2010}). FDis values were higher in less fragmented forest patches, suggesting that mammal species inhabiting them exhibit stronger resource partitioning and exploit the available resources in different ways. Conversely, FEve values were lower, indicating that in those assemblages some parts of the niche space, although occupied, are under-utilized (see ^{Mason et al. 2005}; ^{Villéger et al. 2008}). Accordingly, when values of FEve are higher the use of resources is maximized, while a reduction in these values means that some resources are available to be used (^{R Core Team 2019}). This result suggests that assemblages in less fragmented forested areas may be less resilient to environmental perturbations and less resistant to ecological invasions. The low level of FDis variation between the studied assemblages (mean ± SD=0.469 ± 0.029) is likely due to a similar composition of generalist species (see Table 1). Despite this fact, less fragmented areas support more diverse mammal assemblages with a higher diversity of traits, especially from specialist species (e.g., Panthera onca and Tayassu pecari), resulting in high FDis values.

Forest fragmentation (measured as FAD) was the most important variable explaining FDis. The dispersion of traits increased with increasing FAD values, indicating that forested habitats allow a broader distribution of functional traits in the assemblage. These results suggest that more heterogeneous habi tats offer a wider variety of niches, which might enhance FDis. Moreover, as forest fragmentation decreases a more dispersed distribution of species in the trait space is allowed, probably due to the addition of habitat-specialized species that rely on highly conserved forest for foraging and protection (e.g., Panthera onca and Pecari tajacu). Therefore, the loss of specialist and sensitive species could result in changes in species composition and functional erosion (see ^{Flynn et al. 2009}). Our findings are consistent with other studies where high vegetation structural complexity was related to high values of FDis in terrestrial mammal communities and assemblages around the world (e.g., ^{Flynn et al. 2009}; ^{Ahumada et al. 2011}; ^{Arévalo-Sandi et al. 2018}; ^{Sukma et al. 2019}). Overall, our findings reinforce the importance of retaining and enhance forest land cover in the study area, as landscape modification and forest fragmentation threaten their sustainability and functioning (see ^{Riitters & Wickham 2012}).

Implications for conservation

Most native vegetation in the Middle Magdalena Valley has been transformed at critical thresholds of around 85% of natural vegetation lost (^{Etter et al. 2006}; ^{Garcia-Ulloa et al. 2012}; ^{Garzón & Gutiérrez 2013}). Today, the remaining native forest in the study area are mainly threatened by development activities, such as the expansion of oil palm plantations, petroleum exploration operations, and road networks (see ^{Garzón & Gutiérrez 2013}; ^{Ramos & Meza-Joya 2018}; ^{Meza-Joya et al. 2019}), which may lead to further forest fragmentation, isolation, and loss. Our results suggest that management efforts aimed at protecting the remaining natural forest cover from further loss are urgently needed in this landscape. The most urgent measure is to protect the remnant forest patches to prevent further fragmen- tation, while identifying critical areas for restoration activities aimed to increase landscape connectivity and matrix quality. Increasing the integrity of re maining forest habitat through these conservation actions may be the only available option for conserv ing sensitive and forest interior species.

Overall, our findings reinforce the importance of defining conservation actions for minimizing the impact of human activities on the remaining native forests in the study area, as they are fundamental for maintaining wild mammals and the ecological processes in which they participate. Indeed, heterogeneous landscapes with high-quality matrices within which highly diverse forest communities can persist are essential to establish an integrated and sustainable productive landscape (see ^{Perfecto & Vandermeer 2010}). Interestingly, the landscape composition and configuration of the study area, its potential for agroforestry, and its importance as a national core of economic growth offer an invaluable opportunity for implementing environmental- friendly socio-economic practices aimed to reduce fragmentation and enhance the connectivity of the remaining native forest. However, to reach this goal it is necessary (i) identifying priority conservation areas (e.g., forested areas, ecological corridors) requiring protection and restoration action, (ii) a rigorous, effective, objective, and transparent environmental auditing process to ensure that development projects are complying with environmental policies and compensation measures, and (ii) enforcement of the law and regulations to address illegal deforestation, especially by the expansion of oil palm plantations and cattle.

Supplementary materials