Mechanisms of diversity maintenance in dung beetle assemblages in a heterogeneous tropical landscape

Study area

The study was carried out at the REBISO, located in the municipalities of Ocozocoautla de Espinosa and Cintalapa, Chiapas, Mexico (16°45′42″–17°09′00″N and 93°54′19″–93°21′20″W, Fig. 1). The area is mostly underlain by dolomite rocks and limestone, with a dominance of water-soluble sedimentary rocks (Domenici, 2016). The predominant climate types are warm, humid (climate type Am) and warm, subhumid (climate type Am(f)), with a mean annual temperature of 22 °C and heavy rainfall throughout the year (SEMARNAT/CONANP, 2001).

We produced a vegetation map of REBISO from a multispectral SPOT6 image acquired in 2014, using a supervised classification method in QGIS v2.12.3 (Q GIS Development Team, 2016). The vegetation classes considered were tropical forest, second-growth forest, and pastures (Fig. 1, Table 1).

Sampling design

We established eight 1-km² (100 ha) sampling windows to capture the landscape heterogeneity in the REBISO (Sánchez-de Jesús et al., 2016). Each window was separated from each other by at least 2 km to ensure spatial independence between them (Sánchez-de Jesús et al., 2016). The landscape composition in each window was described by estimating the percent coverage of each vegetation class and evaluating the evenness of their distribution (landscape compositional diversity - Shannon diversity). The spatial configuration of the vegetation classes in each window was assessed with the splitting index and edge density metric (McGarigal, Cushman & Ene, 2012). Edge density is computed as the length (m) of the edges of each vegetation class divided by the window area (ha). The splitting index describes the degree of fragmentation of a landscape and is equivalent to the effective number of patches. Thus, as a landscape becomes increasingly sub-divided, the splitting index increases (Jaeger, 2000; Fahrig, 2017). The landscape composition, edge density, and splitting index metrics (Table S1) were obtained with FRAGSTAT v4.2.1 (McGarigal, Cushman & Ene, 2012). Based on the percent cover of tropical forest (F), windows were classified as intact (W1, W2; F >90%), variegated (W3, W4, and W5; 60% <F <90%), or fragmented (W6, W7, and W8; 10% <F <60%).

We sampled dung beetles (Scarabaeidae: Scarabaeinae) during the dry (March to May) and the rainy (July–August) seasons of 2016 using pitfall traps. This sampling scheme allowed us to integrate the seasonal activities of dung beetles (Cajaiba et al., 2017). Each trap consisted of a 1 L cylindrical plastic container with 300 mL of ethylene glycol as preservative, buried at ground level, and covered with a plastic lid to protect the bait from rain and sun radiation. Pitfall traps were baited with 70 g of either an 80:20 mixture of pig and human feces (copro-traps) or squid flesh (necro-traps) in order to obtain a representative sample of the dung beetle assemblages in the area.

Seven sampling sites were established in each window, separated 250–360 m from each other, to proportionally adjust the number of pitfall traps per vegetation class according to the vegetation class composition of each window (Table S2). Proportional sampling is a suitable method for detecting changes in beta diversity in heterogeneous landscapes (Schoereder et al., 2004). In each sampling site, three copro-traps and three necro-traps (42 traps/window), were placed in a rectangular area separated 50 m from each other to minimize interference between them (Larsen & Forsyth, 2005). The rectangular layout of some trap sets was modified in some cases due to the topographic characteristics of the sites. The pitfall traps were left active for 48 h.

The specimens collected were counted and identified to species. To estimate the dung beetle biomass, we randomly selected ten specimens of each species and dried them at 70 °C for 72 h. We weighed each specimen to the nearest 0.1mg with an analytical balance (Explorer Pro) and calculated the average biomass for each species. Finally, we multiplied the mean biomass of each species by its abundance in each window and vegetation class. The dung beetle specimens were deposited in the entomological collection of El Colegio de la Frontera Sur, San Cristóbal de Las Casas. Field sampling in the REBISO was carried out under permit SGPA/DGS/14214/15 issued by the Secretaria de Medio Ambiente y Recursos Naturales, Mexico.

Data analysis

We followed a spatial and structural approach (sensu Rös, Escobar & Halffter, 2012) to analyze the data. Windows were the sampling units for the spatial approach (n = 8), while vegetation classes within windows were the sampling units for the structural approach (n = 21). The sampling completeness of each window and vegetation class was determined using the coverage estimator of Chao & Jost (2012), which allows comparing species diversity across multiple sites.

Alpha diversity in each sampling unit (window or vegetation class) was evaluated using the ⁰D and ¹D diversity numbers. ⁰D is equivalent to species richness and is insensitive to the species abundance (Jost, 2006); ¹D is equivalent to the exponential of Shannon diversity index and accounts for the most abundant species in a community (Jost, 2006).

We examined differences in species richness between windows by constructing and comparing their 95% bootstrap confidence intervals. Non-overlapping confidence intervals denote significantly different species richness (Gotelli & Colwell, 2011; Chao et al., 2014). Differences in species richness between vegetation classes were determined using interpolation–extrapolation curves (Chao et al., 2014). The sampling coverage, ⁰D and ¹D diversity numbers, confidence intervals, and the interpolation–extrapolation curves were obtained with the software iNEXT v2.0.11 (Hsieh, Ma & Chao, 2016).

Generalized linear models (GLM) were used to assess differences in abundance and biomass between windows and vegetation classes. The abundance and biomass data approached a normal distribution after logarithmic transformation and were analyzed assuming a Gaussian error distribution (Crawley, 2013). Pairwise comparisons using Tukey’s test were carried out, with the multcomp package (Hothorn Bretz & Westfall, 2008) whenever significant differences were detected.

GLMs were also used to assess the effect of the landscape composition and configuration on the species richness (⁰D), exponential of the Shannon diversity (¹D), abundance, and biomass of the dung beetle assemblages. These data were first tested for normality and were then analyzed assuming a Gaussian error distribution. Since only eight observations were available to fit these models, separate models containing only one predictor variable were constructed to avoid overfitting (Kelley & Maxwell, 2003). The best-fit models were selected based on the Akaike′s information criterion corrected for small samples (AICc) and the deviance explained (D²). The model with the smallest AICc (ΔAICc >2) and the largest D² values was selected as the best-fit model (Burnham & Anderson, 2002). Based on the results from the Moran I test (as implemented in the package LetsR), no significant spatial structure was detected in the response variables (Table S3) (Vilela & Villalobos, 2015).

True beta diversity (i.e., the effective number of distinct communities) was estimated for species richness (⁰β) and Shannon diversity (¹β) using the multiplicative partitioning method (Jost, 2007). The multiple-site Sørensen dissimilarity was partitioned as β_Sor  = β_Sim + β_Sne using the package Betapart v1.3 (Baselga & Orme, 2012) to determine whether the ecological differences between sampling units resulted from species turnover (β_Sim) or nestedness (β_Sne). Turnover measures the replacement of species between sites caused by environmental differences, disturbance, or competition. Nestedness is a loss of species between sites, usually due to differences in local conditions or ecological niches, where the species-poorer site contains a subset of the species present in the species-richer site (Legendre, 2014).

Null models were used to determine whether beta diversity patterns resulted from either random changes in alpha and gamma diversity, or from underlying deterministic mechanisms in communities or the landscape (Chase et al., 2011). We constructed null models for the beta Raup-Crick index (β_R−C) using the algorithm developed by Chase et al. (2011) with 9999 randomizations. β_R−C compares the observed versus expected beta diversity under the null model, scaling the results to a range between −1 and 1. This value indicates whether the beta diversity observed between windows, or vegetation classes, is more similar (values close to −1), equal (values close to 0), or less similar (values close to 1) than the one expected by chance (β_R−C null model). We built a dendrogram and a nonmetric multidimensional scaling (NMDS) plot based on β_R−C values for windows and vegetation classes, respectively (Chase et al., 2011). The dendrogram was constructed using the complete linkage method, as it produces clusters with ecological discontinuities (Legendre & Legendre, 2003). We compared the dendrogram and NMDS plot based on β_R−C with homologous plots based on Sørensen dissimilarity to examine whether deterministic mechanisms are underlying the observed beta diversity across the landscape (Chase et al., 2011). All statistical analyses and models were carried out using R v.3.3.1 R Development Core Team (2015).

Diversity, Abundance, and Biomass Patterns in Windows

Species richness (⁰D) in the windows sampled ranged from 22 (W1, intact window) to 37 (W4, variegated window), whereas the exponential Shannon diversity index (¹D) ranged from 4.9 (W2, intact window) to 17.6 (W5, variegated window) species (Table 2). Species richness in W4 and W5 was significantly higher than in the other windows (Fig. 2A).

Deltochilum mexicanum, E. maya, and C. vazquezae were the most abundant species in intact windows W1 and W2 (Fig. S1A); C. vazqueze, E. maya, and Eurysternus angustulus were the most abundant ones in variegated windows W3, W4, and W5; and O. batesi, O. corrosus, and Copris lugubris were the most abundant species in fragmented windows W6, W7, and W8 (Fig. S1A). Biomass patterns in the dung beetle communities differed from those observed in their abundance values. Deltochilum mexicanum, E. maya, and Ontherus mexicanus were the dominant species, in terms of biomass, in intact windows W1 and W2 (Fig. S1A); D. mexicanum, E. maya, and Dichotomius amplicollis were the dominant species in W3; D. mexicanum, D. amplicollis, and Dichotomius annae in W4; and Coprophanaeus corythus, Deltochilum sublaeve, and D. amplicollis in W5. Coprophanaeus corythus, C. lugubris, and D. amplicollis were the species with the highest biomass in fragmented windows W6, W7, and W8 (Fig. S1A).

The highest abundance values (44–52 individuals per trap) were recorded in the intact windows W1 and W2 (Table 3), followed by fragmented windows W6, W7, and W8 (30–36 individuals per trap), and variegated windows W3, W4, and W5 (15–30 individuals per trap). However, these differences were not statistically significant (χ² = 8.923; df = 7; P = 0.26; Table 3). By contrast, there were significant differences in mean biomass between windows (χ² = 45.143; df = 7; P =  > 0.001; Table 3). Mean biomass per trap was significantly higher in windows W1 and W2 (8.2–9.1 grams per trap, Table 3), but no significant differences were found between fragmented (2.2–3.7 grams per trap, Table 3) and variegated windows (2.4–3.9 grams per trap, Table 3).

All the landscape variables had a significant positive effect on the species richness (⁰D) in each window. However, the splitting index was the variable that best explained variations in species richness (Table 4). Although the exponential Shannon diversity (¹D) values were positively related to the splitting index and edge density, the splitting index was the best predictor for variations in Shannon diversity between windows. Edge density and forest cover were the best predictor variables for dung beetle abundance and biomass, respectively; edge density was negatively correlated with abundance, and forest cover positively correlated with biomass (Table 4).

True beta diversity of orders 0 and 1 indicated two effective communities between windows, ⁰β being slightly smaller than ¹β (Table 2). The multiple-site Sørensen value calculated for all windows was 0.65 (Fig. 3A “W Total”); 85% of this dissimilarity was due to species turnover (β_Sim) and 15% to nestedness-resultant component (β_Sne). Sørensen dissimilarity for intact windows W1 and W2 was lower than 0.4, mainly due to nestedness (β_Sne) (Fig. 3A). Dissimilarity ranged from 0.3 to 0.45 in variegated windows (W3, W4, W5), and from 0.3 to 0.58 in fragmented windows (W6, W7, W8). In most cases (except for W8), the observed Sørensen dissimilarity values were primarily due to species turnover (Fig. 3A).

The dendrogram based on the Sørensen distance revealed two main groups (Fig. 3B). The first group includes the fragmented windows (W6, W7, W8), while the second group reveals a gradient of increasing similarity ranging from the intact (W1, W2) to the variegated (W3, W4, W5) windows. The null-model analysis showed that the difference between fragmented windows with respect to the variegated and intact windows was higher than expected by chance (β_RC Value: 1.0, Fig. 3C). However, the dissimilarity between variegated and intact windows did not exceed the null expectation of beta diversity (0<β_RC<0.3, Fig. 3C).

Diversity, abundance, and biomass in vegetation classes

The dung beetle species richness (⁰D = 34 species) in pastures was significantly lower than in the other vegetation classes, but there were no significant differences between second-growth and tropical forests (44 and 45 species, respectively) (Fig. 2B). We recorded the lowest exponential Shannon diversity values (¹D) in the tropical forests (7.75) and the highest in second-growth forests (15.54) (Table 2).

Deltochilum mexicanum, C. vazquezae, and E. maya were the most abundant species in tropical forest sites (Fig. S1B); O. corrosus, D. mexicanum, and E. maya in the second-growth forest; and O. batesi, O. corrosus, and C. lugubris in pasture sites (Fig. S1B). Deltochilum mexicanum, E. maya, and O. mexicanus contributed with the highest biomass in tropical forests; D. mexicanum, C. corythus, and E. maya in the second-growth forests (Fig. S1B); and C. corythus, C. lugubris, and D. amplicollis in pasture sites (Fig. S1B).

No significant differences in abundance were observed between vegetation classes (χ² = 3.701; df = 2; P = 0.16, Table 3). The average number of individuals per trap was 37.8 in second-growth forests, followed by tropical forests (30.2) and pasture sites (17.7) (Table 3). By contrast, there were significant differences in mean biomass between vegetation classes (χ² = 10.829; df = 2; P = 0.004). Pasture sites had a significantly lower mean biomass per trap (2.09 g) than tropical forest (5.05 g) and second-growth forest (4.53) sites, which showed no significant differences between them (Table 3).

According to the multiplicative partition of diversity, there were 2.9 effective communities for ⁰β and 2.5 communities for ¹β in the three vegetation classes combined. Two effective communities were estimated for both the tropical forest and second-growth forest classes, with ¹β higher than ⁰β in both cases. Only one effective community was estimated for the pasture class, with ⁰β higher than ¹β (Table 2). The Sørensen dissimilarity between vegetation classes was 0.85, with 88% of this value accounted for by species turnover (β_Sim) and 12% by nestedness processes (β_Sne) (Fig. 4B).

Tropical forests and second-growth forests showed higher Sørensen values (0.71 and 0.70, respectively) than pastures (0.54) (Fig. 4A). The NMDS plot based on the Sørensen distance formed a compact cluster of pasture sites, whereas most of the tropical forest and second-growth forest sites overlapped between themselves and with the pasture sites (Fig. 4B). The NMDS plot based on the beta Raup-Crick null model index (β_R−C) separated the tropical forest sites from pastures, whereas second-growth forest sites overlapped with tropical forest and pasture classes (Fig. 4C).

Local patterns of species richness and assemblage structure

Dung beetle communities in variegated windows showed the highest richness values in the REBISO. Rös, Escobar & Halffter (2012) and Costa et al. (2017) also found a higher richness of dung beetle species in variegated landscapes. Landscape variegation can be a significant environmental driver of local diversity as it increases the range of habitats available for species by creating a complex composition and configuration (Tscharntke et al., 2012; Ramírez-Ponce et al., 2019). The high species richness and diversity found in variegated windows in the REBISO can be attributed to the convergence of multiple dung beetle assemblages including forest specialists (e.g., Eurysternus caribaeus, Sulcophanaeus chryseicollis, and Uroxys boneti), forest-pasture edge specialists (O. landolti, Canthon cyanellus), open habitat specialists (D. annae, C. lugubris, O. corrosus), and generalist beetles (O. batesi) (Favila, 2005; Navarrete & Halffter, 2008).

Intact and fragmented windows showed lower diversity values than variegated windows. This diversity pattern is consistent with the intermediate disturbance theory (Grime, 1973). Sites with little or no disturbance favor the predominance of highly competitive forest specialists such as D. mexicanum, C. vazquezae, and E. maya, which accounted for 85% of the total abundance and 90% of the total biomass in intact windows, thus preventing a higher local diversity. On the other hand, the intense landscape changes caused by livestock production in fragmented windows reduce the local species richness of dung beetles since many native-forest species are unable to adapt to the new open habitat conditions (Halffter, Favila & Halffter, 1992; Silva, Storck-Tonon & Vaz-de Mello, 2016; Alvarado et al., 2018).

The presence of the exotic African species Digitonthophagus gazella (Montes de Oca & Halffter, 1998) in the REBISO is worth mentioning. Although D. gazella was only recorded in pastures of fragmented windows (W6, W7), and contributed with only a small fraction of the community abundance and biomass (six and four percent, respectively), they may pose competitive pressure on native species inhabiting open areas (Lobo & Montes de Oca, 1994). Further studies are needed to assess how this invasive beetle might affect native species in the REBISO.

Dung beetles are involved, among other ecological processes, in the recycling of organic matter, soil bioturbation, and secondary seed dispersal (Nichols et al., 2008). The amounts of soil removed, dung buried, and seed dispersed are significantly and positively influenced by the species richness and biomass of dung beetle assemblages (Nunes et al., 2018; Alvarado, Dáttilo & Escobar, 2019). The tropical forest sites showed the highest dung beetle species richness and biomass values. Besides, forest coverage was positively related to dung beetle biomass in the REBISO. Both results indicate that the tropical forest sites likely contain the most functionally efficient dung beetle assemblages, thus emphasizing the importance of forest conservation in the REBISO.

Effects of landscape composition and configuration on dung beetle assemblages

Previous studies conducted in tropical ecosystems have identified landscape composition as the main predictor of the diversity of dung beetle assemblages (Sánchez-de Jesús et al., 2016; Alvarado et al., 2018). However, in our study, landscape fragmentation was the primary explanatory variable of variations in species richness and diversity. These findings are likely due to the variegated structure of the REBISO landscape and its “biodiversity-friendly” matrix of second-growth forest (see Perfecto & Vandermeer, 2008; Melo et al., 2013). First, second-growth forests in the REBISO are structurally similar to forest habitats (Ramírez-Marcial et al., 2017). Therefore, while many dung beetle species are restricted to forest patches, others may persist and use the second-growth forest matrix to move between forest patches (Díaz, Galante & Favila, 2010). Second, fragmentation in variegated environments creates conditions that allow the coexistence of species from different habitat types (e.g., forest species, pasture species, edge specialist species), thereby increasing the diversity of dung beetles at the landscape scale (Villada-Bedoya et al., 2016; Fahrig, 2017).

Not all fragmentation effects are beneficial since a higher edge density can have adverse effects on the abundance, biomass, and even the physiological condition of tropical dung beetles (Portela Salomão et al., 2018). We found the lowest abundance of dung beetles in variegated windows, where the highest edge density occurs. A higher edge density is coupled with less habitat area, limiting the capacity of the landscape to support medium and large-sized mammal species. Pozo-Montuy et al. (2019) observed that medium- and large-sized mammals are significantly less abundant and diverse in the REBISO buffer zone (i.e., where the variegated windows are located). Such reduction in mammal density can cause a marked decrease in dung quantity and availability, thus limiting the growth of dung beetle populations (Nichols et al., 2009). Also, microclimatic conditions such as temperature and relative humidity are more variable in forest edges, which might negatively affect the reproduction and survival of dung beetles (Klein, 1989; Feer, 2013). Our findings suggest that fragmentation processes in variegated windows foster a high dung beetle diversity, but might also limit their population growth due to insufficient resources, reduced habitat area, or sub-optimal microclimatic conditions. Future studies should assess the strength and extent of this trade-off between dung beetle diversity and abundance, and its functional consequences across the REBISO.

Beta diversity patterns and mechanisms of diversity maintenance

Species turnover is the primary driver of the high diversity and complementarity of the dung beetle communities found in the REBISO. There were between 3 and 27 species not shared between windows, and from 4 to 35 species not shared between vegetation classes. Each window and vegetation class contributed two or three unique species to the overall diversity. The largest turnover values were found between the fragmented windows (W6, W7, and W8) vs. the variegated and intact windows (W1 to W5), and between the forested vegetation classes (tropical forest, second-growth forest) vs. the pasture sites. The anthropic disturbances and the heterogeneous landscape of REBISO favor this high beta diversity since dung beetles are especially susceptible to environmental variability (Arellano, Leon-Cortes & Halffter, 2008; Costa et al., 2017).

The differences observed in the species assemblages of fragmented windows (W6, W7, and W8) and those in the other windows (W1 to W5) are not random. Likewise, the differences between tropical forest and pasture assemblages are not random. Significant deviations from random expectations of beta diversity indicate niche-structured assemblages in which environmental filters determine species membership in a community (Chase et al., 2011; Püttker et al., 2015). Intensive anthropic disturbances such as deforestation can act as an environmental filter in fragmented windows and pastures, selecting stress-tolerant dung beetle species able to survive in open habitats (Halffter, Favila & Halffter, 1992; Spector & Ayzama, 2003; Gardner et al., 2008). Deltochilum mexicanum, C. vazquezae, S. chryseicollis, Canthon femoralis, and E. maya probably are the species most sensitive to the environmental filters caused by anthropic disturbance. Although these forest species are widely distributed in the biosphere reserve (Sánchez-Hernández et al., 2018), their abundance was drastically reduced in fragmented windows.

We found signs of biotic homogenization in the pasture sites. For instance, the lowest alfa and beta diversity values were recorded in pastures, and their species assemblages were more similar to each other than expected by chance, regardless of the windows where they were located, indicating shared environmental filtering processes (Chase, 2010). Anthropogenic environmental filters are one of the main drivers of biotic homogenization, eroding alfa and beta diversity and diminishing ecosystem resilience and viability (Gámez-Virués et al., 2015). Hence, the advance of the agricultural frontier in the REBISO landscape should be monitored closely to prevent further biotic homogenization processes among the dung beetle species assemblages.

In our study, ¹β between the intact and variegated windows (W1 to W5), as well as between tropical forests and second-growth forests, was higher than ⁰β. Thus, the true beta diversity is mainly due to differences in the abundance of shared species rather than to differences in richness (Jost, 2007). Besides, the overall beta diversity between these windows and vegetation classes was not different from that expected by chance. Most species in neutral communities are considered ecologically equivalent since, in the absence of any factor limiting their dispersal, they can appear at random in any of the null assemblages (Püttker et al., 2015; Ortega-Martínez et al., 2020). Both results suggest that the REBISO still holds sufficient vegetation cover to maintain a continuous flow of dung beetles between forested landscape sections (W1 to W5).

Given the significant stochasticity of beta diversity between intact and variegated windows, and between tropical forests and second-growth forests, we can conclude that the landscape variegation in the REBISO does not affect dung beetle diversity negatively. However, it is essential to conserve the forested patches to maintain a high dispersal between sites, thereby increasing the resilience of dung beetle populations to habitat loss and isolation (De Castro Solar et al., 2015; Socolar et al., 2016). Landscape variegation can be an effective strategy to buffer the impact of intense anthropic disturbances (Rös, Escobar & Halffter, 2012; Costa et al., 2017). Variegation can be achieved by maintaining the forest cover and incorporating more biodiversity-friendly production systems, such as agroforestry practices, in the landscape matrix (Perfecto & Vandermeer, 2008).