Globally, tropical rain forests comprise some of the most diverse and functionally rich ecosystems but are increasingly degraded by human impacts. Protected areas have been shown to conserve species diversity, but their effectiveness at maintaining functional diversity over time is less well known, despite the fact that functional diversity likely reveals more ecological information than taxonomic diversity. By extension, the degree to which species loss decreases functional diversity within protected areas is also unknown; functional redundancy may buffer communities from loss of functional diversity from some local extinctions. Using eight years of camera trap data, we quantified annual functional dispersion of the large mammal community in the Volcán Barva region of Costa Rica and tested for changes in functional dispersion over time in response to environmental and anthropogenic predictors. We quantified functional redundancy based on modeled declines in functional dispersion with species loss. Functional dispersion did not change significantly over time and was not associated with measured environmental or anthropogenic predictors. Quantitative modeling of change in functional traits over time did not identify significant changes. We did however find qualitative trends in relative trait proportions, which could be indicative of functional change in the future. We found high functional redundancy, with average functional dispersion declining significantly only after 9 out of 21 large mammal species were lost from the community. We cautiously suggest that protected tropical rain forests can conserve functional diversity over the course of a decade even in heavily fragmented landscapes.

METHODS

             

1. Study taxa

Large mammals have a disproportionate impact on their ecosystem becuase of their large body size, dietary requirements and home ranges (Bakker et al. 2016). For these reasons, large mammals also have the potential to respond strongly to environmental change based on their functional traits (Chiarello 1999, Newbold et al. 2014). For millenia, large-bodied mammals have consistently exhibited the highest extinction rates (Dirzo et al. 2014), particularly within fragmented tropical landscapes (Crooks 2002; Jorge et al. 2013). Tropical mammal communities may therefore be particularly vulnerable to functional diversity loss in a changing environment. 

            We analyzed a terrestrial large mammal community, which had 21 species (Beaudrot et al. 2016a). We included all predominately ground-dwelling mammal species over 100 grams in average body mass. There were several additional co-occuring mammal species but they either fell outside this definition or were excluded due to very low detection (Table S3). Functions performed by large mammal species were considered functionally distinct in this study, but may also be provided by other taxa (e.g. arboreal mammals, birds, etc.).

 

2. Study site 

Volcán Barva in Costa Rica consists of tropical rain forest that encompasses La Selva Biological Research station and Braulio Carrillo National Park; La Selva has a long, rich history of ecological research (Pringle et al. 1984). With camera trap points ranging in elevation from 49 to 2569 meters above sea level (Ahumada et al. 2013), and an elevational coefficient of variation of 1.01, the survey area encompasses considerable habitat heterogeneity. The terrestrial mammal community at this site has been monitored annually with camera traps since 2007 as part of the Tropical Ecology Assessment and Monitoring (TEAM) Network. TEAM was established to monitor ground-dwelling mammal and bird communities in protected tropical forests worldwide using a standardized camera trapping protocol (Jansen et al. 2014). We selected Volcán Barva because it is the longest running TEAM site. There is no evidence of large mammal species loss at Volcán Barva since 2007 (Beaudrot et al. 2016b), but populations of some mammal species have declined (Ahumada, Hurtado & Lizcano 2013). 

            Human disturbance is also prevalent in and around Volcán Barva. Poachers are a threat to mammal species in the area, and deforestation occurs along the border of Braulio Carrillo National Park (Rovero & Zimmerman 2016; Schelhas & Sánchez-Azofeifa 2006). Over 50% of the park’s border is classified as fragmented landscape (Ahumada et al. 2011). Nevertheless, little forest cover change has been observed within the protected area in recent years (Beaudrot et al. 2019), which suggests that human disturbance from illegal logging has been minimal.

 

3. Data collection

 

1. Camera trap data and occupancy estimates

We used camera-trap data collected by TEAM between 2007 and 2014. TEAM surveys terrestrial (i.e., ground-dwelling) tropical mammal populations on an annual basis, using a standardized protocol with large-scale arrays of permanent camera-trap points (Jansen et al. 2014). Sixty camera traps were deployed at a density of 1 camera per 1 - 2 km², encompassing a survey area of 21,049 hectares out of the 49,317 hectares within the boundaries of the protected area. Each camera trap was activated for 30 consecutive days annually at the same time every year during the dry season. This was done to account for seasonality across TEAM sites and across years. While mammal community compositions can shift in response to seasonal differences such as plant productivity and fruit availability (Ramírez-Bautista & Williams 2019, Wen et al. 2014, Marshall et al. 2014), the consistent temporal camera trap deployment each year likely reduced any potential bias. Camera trap images were identified by TEAM personnel following the standard IUCN Red List (IUCN 2014).

Species-specific annual occupancy values for the Volcán Barva community were obtained from a previously published study (Beaudrot et al. 2016b), which used single species Bayesian dynamic occupancy modeling that provided a posterior distribution of 1000 occupancy values for each species for each year based on the TEAM camera trap data. We selected the median species-specific posterior occupancy value for each species for use in this study. 

 

2. Functional trait data

We obtained mammal functional trait data through an extensive literature search in which all monitored terrestrial mammal species were assigned ranked values for six functional traits: average body mass, diet, average social group size, habitat type, activity period, and average litter size (Table 1; Table S1; Table S2). We selected these functional traits for their association with aspects of individual species ecology relevant to how species utilize their environment and impact their ecosystem (Weiss & Ray 2019). For example, body mass affects the quality of resources needed for survival (Jarman 1974) and approximates the degree of impact that species will have on its ecosystem in terms of quantity of nutrients dispersed, quantity of food consumed and spatial range of impact (Wolf, Doughty & Malhi 2013). The selected traits are used extensively in studies of mammal functional traits (Flynn et al. 2009, Hempson, Archibald & Bond 2015, Jones et al. 2009). Nevertheless, functional traits studies are highly dependent on the selection of relevant functional traits (Petchey & Gaston 2006), and although we cover a wide suite of ecological attributes, it is possible that important traits have been inadvertently omitted. In addition, although these ordered traits cover a breadth of functional aspects, their categorical nature may obfuscate patterns that would be revealed by continuous trait measurements (Kohli & Rowe 2019). 

 

3. Environmental data

We collected data on climatic, biological and anthropogenic factors that could impact mammal functional diversity (Table S4). We selected four variables that were not strongly correlated with other considered variables (r  < 0.6) as potential predictors of mammal functional diversity. The four predictors were 1) annual precipitation, 2) area of new canopy gaps within Volcán Barva each year, 3) mean area of forest fragments within the Zone of Interaction (Defries et al. 2010a) each year and 4) area of forest loss within the Zone of Interaction each year. The Zone of Interaction is the spatial extent believed to most likely affect biodiversity within the sampling area and is systematically quantified based on human settlements, watersheds and migration corridors (Defries et al. 2010a; Beaudrot et al. 2016b).

We selected the above variables for the following reasons. Annual precipitation can affect plant and animal functional traits because water is required by all organisms for metabolism and is a limiting resource for many (Wright et al. 1999; Dwyer & Laughlin 2017). Canopy gaps are critical components of vegetation dynamics and represent important aspects of plant community diversity (Denslow 1987). For example, approximately 75% of the tree species in La Selva are dependent on canopy gaps for seed germination and growth (Hartshorn 1978). Vegetation dynamics, in turn, have the potential to strongly affect mammal community functional traits by altering resource availability and habitat structure (Laurance 1991; Laurance et al. 2008). Edge effects and isolation from fragmentation have been shown to affect community composition in other systems (Newmark 1987; Malcolm 1994; Krishnadas et al. 2018). Deforestation has also been shown to be one of the primary drivers of defaunation in the tropics (Canale et al. 2012) and has the potential to disturb community structure.

Mean annual precipitation data were collected from the NASA POWER project at a 0.5°x0.5° resolution (POWER data access viewer). We used remotely sensed vegetation classification data (Hansen et al. 2013) to calculate canopy gaps, mean fragment size, and forest loss over the study period (package ‘SDMTools’ in R; VanDerWal et al. 2014).

 

4. Analysis

 

1. Overview

First, we calculated annual functional diversity for the mammal community in the Volcán Barva region of Costa Rica for an eight-year period and assessed these values for linear trends over time. We then used linear regression to examine the relationship between temporal change in functional diversity and our four environmental predictor variables. We tested for temporal trends in occupancy-weighted trait values to examine quantitatively how individual functional traits changed over time within the Volcán Barva large mammal community. Lastly, we used bootstrapping to quantify functional redundancy and segmented linear regression to analyze how functional diversity changed with simulated removal of species from the community. 

 

2. Functional traits

We selected functional traits for assessment and ranking based on established methods for examining effects of anthropogenic change on mammalian functional traits that reflect species responses to environmental conditions (response traits) (Flynn et al. 2009, Díaz et al. 2013). The selected traits were also doubly valuable as they are associated with functional impacts of mammals (effect traits) on the environment (Hempson, Archibald & Bond 2015). Two of the trait variables were continuous (i.e. average body mass, litter size) and one was an ordered category (i.e. average social group size). For the remaining three traits, we imposed ordered categories following previous work on mammal functional traits (Hempson, Archibald & Bond 2015, Jones et al. 2009, Flynn et al. 2009). Specifically, we ordered diet from the lowest quality food (grass - grazers) to the highest quality food (vertebrate meat - carnivore), habitat from the most horizontally oriented (terrestrial) to the most vertically oriented (scansorial), and activity period from highest light (diurnal) to lowest light (nocturnal). The traits with natural or imposed ordered categories (diet, social group size, habitat, activity period) were analyzed as ordinal variables in the functional diversity calculations and in our assessment for temporal trends. No traits were strongly correlated with each other (r < 0.6).

 

3. Functional diversity            

We used the abundance-weighted functional dispersion metric (package `FD` in R; Laliberte & Lengendre 2010) to calculate annual functional diversity in Volcán Barva over the eight-year study, using occupancy values as a proxy for abundance weights. We determined that a functional diversity metric based on functional dispersion would be the most effective for understanding the functional effects of this community on its ecosystem based on extensive work showing the link between community functional dispersion and ecosystem functioning (Cadotte 2017, Frainer et al. 2014). Functional dispersion is the mean distance of individual species to the centroid of all species in a community along functional trait dimensions (Laliberte & Lengendre 2010). Because species occupancy changed and species richness remained constant over the eight-year study period (Beaudrot et al. 2016b), we determined that an investigation into a potential trait-abundance shift measured with an occupancy-weighted functional dispersion metric would give the most meaningful results for assessing functional diversity within a single site (Boersma et al. 2016). Laméris and colleagues (2019) found differences in functional dispersion in large mammal communities in Cameroon based on conservation efforts, further justifying our selection of this metric. 

 

4. Linear modeling

To test for change in functional dispersion over time, we ran a simple linear regression model with year as the predictor. We evaluated the estimate, standard error and p-value of the year coefficient from this model. We also constructed linear models using Gaussian distributions and performed model selection (package `MuMIn` in R; Barton 2019) to evaluate the predictive power of environmental and anthropogenic variables for functional dispersion. The environmental variables were annual precipitation, new canopy gap area within Volcán Barva, mean fragment size within the Zone of Interaction of the protected area, and annual rate of deforestation in the Zone of Interaction. We used AICc model selection to compare models and defined the best model as the model that had the lowest AICc value by a margin of 2 or more (Anderson & Burnham 2004).

 

5. Functional trait distributions over time 

To test for changes in individual functional traits over time, we ran linear regression models for quantitative traits (i.e., body mass, litter size) and ordinal regression models for ordered traits (i.e., diet, group size, habitat, activity period) with year as a predictor. We weighted the response variable for each trait based on the occupancies of the species with the given trait. We used the estimates, standard deviations and p-values to assess the statistical significance of temporal trends.

            

6. Functional redundancy

To estimate functional redundancy, we used occupancy and functional data to simulate how community functional dispersion changed as species were removed from the community. We calculated this change with bootstrapping, randomly drawing subsets of species from the community and calculating their functional dispersion (see above). Species richness values ranged from 2 to N-1, where N was the number of species in the full community (N = 21). We randomly selected 1000 species combinations without replacement of each species richness value to generate a distribution of functional dispersion values and used the mean functional dispersion value in our calculations. For higher species richness, repetition of species combinations was necessary in the analysis to generate 1000 functional dispersion values, but the combinations were still selected randomly without replacement. We performed this analysis for each year using occupancy values from the respective year of the eight-year study period.

We fit a segmented linear regression model (package `strucchange` in R; Zeilies et al. 2002) to the bootstrapped functional dispersion models from all eight years and performed model selection to determine the number of break points (0-4) in the most robust regression. The best-fit model was determined by the lowest AIC value. We used break points to identify functional redundancy, or the number of species lost from the community before rate of functional dispersion loss increased, if ever. 

Code for using these data to calculate values found in the paper is available on Github