The study was conducted in the Serranía del Perijá, Cesar department, Colombia, in two villages of La Victoria district (La Jagua de Ibirico; 9°33′40″N, 73°20′11″W) and two villages of Estados Unidos district (Becerril; 9°42′11″N, 73°16′39″W) (Fig. 1). Sites ranged from 800 to 1500 m a. s. l., with vegetation typical of transitional zones between tropical dry and sub-humid forests. Mean annual temperature is 24 °C, with 1000–2000 mm of rainfall. The climate has two rainy seasons: moderate (April–June) and intense (September–November) (Rangel et al., 2019).

Fig. 1.

Study area in the Serranía del Perijá, northern Colombia. The four sampling windows are shown (W1–W2 in La Jagua de Ibirico; W3–W4 in Becerril, black polygons), along with sampling points in forest (green), regeneration (blue), and crops (red).

To represent the regional variation of the Serranía del Perijá and locate the sampling sites, four windows (W1–W4) of approximately 4 km² each were delineated, separated on average by 1.5 km (Fig. 1). A high-resolution GeoEye-1 satellite image (spatial resolution: 50 cm; spectral: 4 bands—Blue, Green, Red, NIR1) was used to identify land-cover composition and configuration in each landscape window. Land-cover classes were visually interpreted and manually digitized following CORINE classification criteria (IDEAM, 2010). From the original classification, three categories were grouped for sampling: (1) Forest (For) –only the forest category; (2) Crops (Cro) – temporary crops, permanent crops, and heterogeneous agricultural areas; and (3) Regeneration (Reg) – secondary vegetation dominated by herbaceous and shrub species. The area of each land-cover type is presented in the supplementary material (Table S1).

Landscape composition was assessed using the percentage of area (PLAND), which indicates the proportion of the landscape occupied by each cover class. Configuration metrics included Shape index (Shape), which quantifies patch shape complexity relative to a standard shape of the same area, and nearest-neighbor distance (DistE), which measures the degree of isolation among patches of the same class. All landscape metrics were calculated as class-level metrics within each landscape window using the landscapemetrics package (Hesselbarth et al., 2019) to characterize the composition and configuration of the dominant land-cover types across the heterogeneous mountain mosaic. Additionally, at each sampling point we recorded canopy cover (CanCov) using a spherical densitometer and altitude (Alt, m a. s. l.) with a GPS to characterize local conditions.

In each window, 12 sampling points were established (four per land-cover type), totaling 48 points (Fig. 1). Points were distributed to represent the dominant land-cover types and ensure accessibility for repeated sampling. Three sampling campaigns were conducted between March 2021 and May 2022.

At each point, we set up a transect with four baited traps, spaced 50 m apart, and one funnel-type light trap. The baited traps consisted of 1 L containers filled with decomposing fruits (Musa paradisiaca L., Carica papaya L., Mangifera indica L.) mixed with beer (attractant), suspended 2−5 m above ground and active for three days. Light traps were automated with photocells and powered by solar panels, operating for two nights.

Collected specimens were preserved in 95% ethanol, counted, and identified to species level when possible, or otherwise to genus. Identifications were verified by Jhon César Neita (Instituto Alexander von Humboldt). Voucher specimens were deposited in the entomological collections of the Instituto Alexander von Humboldt and Universidad del Atlántico (UARC).

Analyses were conducted at two spatial scales. At the local scale, diversity was estimated within each landscape window using the four sampling points per land-cover type to characterize spatial variation among landscape units. At the regional scale, analyses were based on pooled assemblages combining all windows within each land-cover category. Thus, regional estimates represent landscape-scale patterns rather than comparisons among independent geographic regions. All analyses were performed in R v.4.3.1 (R Core Team, 2023). Sampling representativeness was assessed using sample coverage (Ĉn) with the iNEXT package (Hsieh et al., 2022), which estimates the proportion of species recorded relative to the total expected (Chao et al., 2014).

Alpha taxonomic diversity was estimated using Hill numbers expressed as effective species: °D (richness), 1D (diversity of common species), and 2D (diversity of dominant species) (Jost, 2006), calculated with iNEXT. Hill numbers were compared among land-cover types at both scales. At the local scale, comparisons were performed within each landscape window, whereas regional comparisons were based on pooled assemblages across all landscape windows. Differences were assessed with Kruskal–Wallis tests, followed by Mann–Whitney comparisons.

Beta diversity was calculated among land-cover types at both spatial scales using Sørensen’s index, partitioned into turnover (βtur) and nestedness (βnes) components (Baselga, 2010), using the betapart package (Baselga et al., 2023).

At least three individuals per species and land-cover type were selected to measure six functional traits: eye diameter (ED), hind leg–body length ratio (LeBoRat), biomass (BIO), wing aspect ratio (WiHeRat), trophic group (TrGr), and wing loading (WL) (see values in Table S2). These traits were selected because they are associated with resource use, dispersal capacity and habitat use in phytophagous Scarabaeidae beetles (García-Atencia et al., 2024). Functional richness (FRic), evenness (FEve), and dispersion (FDis) were calculated (Villéger et al., 2008; Laliberté and Legendre, 2010) using the “FD” package. To control for richness-related biases, standardized effect sizes (SES) were computed using null models (Mouchet et al., 2010; Swenson, 2014). Differences were tested with Kruskal–Wallis + Dunn (FRicSES) and ANOVAs + Tukey (FEveSES, FDisSES).

Functional beta diversity was calculated among land-cover types at both spatial scales using Sørensen’s index and its partition into turnover and nestedness components (Villéger et al., 2013). In this case, turnover reflects the replacement of trait composition among assemblages, whereas nestedness indicates that communities differ mainly by containing subsets of the functional trait space present in richer assemblages. A Gower distance matrix was constructed from the functional traits, followed by a principal component analysis (PCA) to reduce trait dimensionality while preserving the main variation in trait composition. The first three axes, summarizing the main variation in functional traits, were retained as synthetic traits to construct the multidimensional functional space and calculate functional diversity indices (Villéger et al., 2011). Functional space was estimated using convex hulls, and all calculations were performed with the betapart package.

To identify which landscape composition and configuration variables were associated with taxonomic and functional diversity patterns across land-cover types, landscape variables (PLAND, Shape, DistE, CanCov, and Alt) were summarized with PCA. Axis 1 was mainly associated with landscape composition and patch shape, whereas Axes 2 and 3 were primarily related to altitude and patch isolation metrics (see contributions in Table S4). The retained PCA axes were subsequently used as predictors in Generalized Additive Mixed Models (GAMMs) (Wood, 2017). Alpha diversity metrics (°D, 1D, 2D, FRicSES, FEveSES, FDisSES) were modeled as a function of these axes, with land-cover type as a fixed factor and window as a random effect. A Poisson distribution was assumed for °D and Gaussian for the other metrics. Preliminary model evaluation indicated no evidence of overdispersion, supporting the use of the Poisson error distribution.

The relationship between environmental dissimilarity and beta diversity (taxonomic and functional) was assessed using GAMMs relating Sørensen dissimilarity matrices to Euclidean distance matrices calculated from the original environmental variables, using a beta regression family (“betar”) appropriate for proportional response variables bounded between 0 and 1.

At the regional scale, based on pooled assemblages across all landscape windows, taxonomic alpha diversity did not differ significantly among land-cover types (p > 0.05). At the local scale, species richness (0D) varied among land-cover types in some landscape windows, but not consistently. In these cases, values of 1D and 2D tended to decrease, indicating that assemblages were dominated by a few abundant species (Fig. 2).

Fig. 2.

Taxonomic and functional diversity of phytophagous scarab beetles across different land-cover types. (Top) Taxonomic alpha diversity expressed as effective numbers of species for Hill numbers 0D, 1D, and 2D in the four landscape windows (W1–W4), each considered a distinct landscape unit, and at the regional scale, represented by pooled assemblages across all windows within each land-cover type. (Middle) Standardized effect size (SES) values of functional diversity (FRic), functional evenness (FEve), and functional dispersion (FDis); red dashed lines indicate the expected range from –1.96 to 1.96. (Bottom) Taxonomic and functional beta diversity among land-cover types (Cro = crops, For = forest, Reg = regeneration), partitioned into turnover (βturn, blue) and nestedness (βnes, red). Points and error bars represent means and standard errors.

Functional evenness (FEveSES) was significantly higher in crops than in forests (p < 0.05), while functional dispersion (FDisSES) did not show significant differences but tended to be higher in forests and regeneration. Standardized functional richness (FRicSES) did not vary among land-cover types (p > 0.05). At the local scale, FRicSES values fell below expected intervals in some forest and crop sites, while FEveSES in regeneration showed high variability among points (Fig. 2).

Taxonomic beta diversity among land-cover types, based on pooled assemblages across all landscape windows, was below 20% and was dominated by species turnover (βtur). The highest local dissimilarity (55%) was observed between crops and regeneration in W1, followed by forest and regeneration in the same window, whereas the lowest (25%) corresponded to the forest–crop pair in W2. In this latter window, nestedness explained a considerable fraction of dissimilarity, indicating that regeneration represented a subset of the other land-cover types. Turnover predominated in the remaining comparisons (Fig. 2).

Functional beta diversity was entirely driven by nestedness. High local dissimilarity occurred between forest and regeneration in most windows, and between forest and crops in one of them (Fig. 2). This suggests that communities share dominant functional traits but lose attributes in less diverse habitats.

The PCA showed that Axis1 (42.6% of variance) was associated with canopy cover, percentage of area, and patch shape; Axis2 (23.8%) with altitude and distance between fragments; and Axis3 (15.4%) with altitude and shape (Table S4).

GAMMs indicated that altitude and distance between fragments (Axis2) influenced species richness (0D) across all three land-cover types (Table S5, Fig. 3). In forests, richness increased at lower altitudes and shorter distances; in crops, it showed a nonlinear relationship with both variables; and in regeneration, richness was higher at intermediate altitudes and in more irregular patches. The 1D and 2D metrics showed similar trends, with increases in crops and decreases in forests at higher altitudes.

Fig. 3.

Relationship between diversity metrics and landscape variables. (A–B) PCA biplots for axes 1–2 and 2–3, showing the contribution of landscape composition and configuration variables (PLAND, CanCov, Shape, Alt, DistE) (see contributions in Table S4). Points represent sampling sites in forest (For), crops (Cro), and regeneration (Reg) across the four windows (W1–W4). (C–N) Results of significant GAMMs (p < 0.1), where the y-axis corresponds to the smoothed function estimated for each taxonomic (0D, 1D, 2D) and functional (FRicSES, FEveSES, FDisSES) diversity metric across the three land-cover types. Detailed results by spatial scale are shown in Table S5.

Regarding functional diversity, FRicSES in crops was negatively related to altitude and distance between fragments, whereas FEveSES was lower at higher altitudes and in more regular patches. FDisSES increased in forests with altitude and in regeneration with more irregular patch shapes (Fig. 3, Table S5).

Fig. 3. Relationship between diversity metrics and landscape variables. (A–B) PCA biplots for axes 1–2 and 2–3, showing the contribution of landscape composition and configuration variables (PLAND, CanCov, Shape, Alt, DistE) (see contributions in Table S4). Points represent sampling sites in forest (For), crops (Cro), and regeneration (Reg) across the four windows (W1–W4). (C–N) Results of significant GAMMs (p < 0.1), where the y-axis corresponds to the smoothed function estimated for each taxonomic (0D, 1D, 2D) and functional (FRicSES, FEveSES, FDisSES) diversity metric across the three land-cover types. Detailed results by spatial scale are shown in Table S5.

Beta-diversity patterns were evaluated using environmental dissimilarity matrices. Taxonomic beta diversity was mainly associated with altitudinal dissimilarity, patch shape, and fragment isolation, whereas taxonomic turnover (βtur) was primarily related to altitude (Fig. 4A, Table S6). Functional beta diversity was strongly associated with PLAND, canopy-cover, altitude, and fragment isolation (Fig. 4B, Table S6). Functional beta diversity was predominantly driven by nestedness, whereas turnover values remained consistently low across comparisons.

Fig. 4.

GAMM responses of beta diversity to landscape dissimilarity variables. (A) Taxonomic turnover vs. altitudinal dissimilarity. (B) Functional beta total vs. PLAND dissimilarity. Solid lines represent fitted smooth functions and dashed lines 95% confidence intervals.

Communities were dominated by Melolonthinae, a group widely associated with disturbed and heterogeneous tropical landscapes (García-Atencia and Amat-García, 2021). This dominance at the landscape scale has been rarely reported and suggests that landscape structure favors species with rhizophagous and saprophagous feeding habits. We also recorded rare and poorly known species in Colombia, restricted to forests and regenerating vegetation, underscoring the value of these habitats for conserving low-abundance and narrowly distributed taxa (Mouillot et al., 2013).

Altitude and fragment isolation were significant predictors of both taxonomic and functional diversity. Higher diversity at higher elevations may reflect environmental conditions that favor Melolonthinae, as previously reported in Colombian mountain agroecosystems (Pardo-Locarno et al., 2005). The positive effect of fragment isolation suggests that greater distances among patches contribute to community differentiation by limiting dispersal and reducing landscape connectivity (Fahrig, 2003).

Patch shape also influenced diversity. Small and irregular patches may intensify edge effects (Murcia, 1995; Riva and Fahrig, 2022), particularly in regenerating vegetation. These findings are consistent with studies showing that habitat loss has a stronger influence than landscape configuration on scarab beetle richness (Sánchez-de-Jesús et al., 2016; Rivera et al., 2021).

Turnover was the main component of taxonomic beta diversity, particularly along altitudinal gradients, suggesting that environmental filtering may promote species differentiation among land-cover types, together with dispersal limitation in fragmented landscapes (Rivera et al., 2021; Marquez-Peña and Domínguez-Haydar, 2023). At local scales, this high turnover among land-cover types indicates that forest, regeneration and crop assemblages differ in species composition, possibly due to variation in vegetation structure and resource availability. At regional scales, the increase in nestedness suggests compositional simplification and species loss with disturbance. Similar patterns have been associated with habitat fragmentation and environmental filtering in scarab beetles from other Neotropical landscapes (Rivera et al., 2021). However, additional studies evaluating dispersal and resource use across landscape gradients would help clarify the relative contribution of these mechanisms.

Functional beta diversity was dominated by nestedness, suggesting reduced functional differentiation among communities. This pattern may reflect functional homogenization associated with the persistence of generalist species and the reduction of specialists linked to more complex vegetation structures (Rivera et al., 2021). Consequently, ecological functions associated with phytophagous scarab beetles could become simplified under disturbance, potentially reducing ecosystem resilience and functional complementarity (Anderson et al., 2011). These results highlight the importance of incorporating a functional perspective into studies of phytophagous scarab beetles, as changes in functional structure may not be detected through species richness alone.

Our results highlight the importance of maintaining structural and functional connectivity in fragmented landscapes. Patch distance and shape are crucial factors to consider in management plans, prioritizing the preservation of small patches that contribute disproportionately to biodiversity (Arroyo-Rodríguez et al., 2020; Riva and Fahrig, 2022). Ecological restoration strategies at the farm scale, such as conserving riparian forests, establishing live fences, and promoting low-intensity production matrices, can generate biological corridors that facilitate movement and reduce homogenization (Fahrig, 2002).

The predominance of functional nestedness further suggests reduced functional differentiation among communities under disturbance. Maintaining functionally diverse assemblages may help sustain ecological processes associated with phytophagous scarab beetles in agricultural systems (Morón et al., 1997) and reduce the risk of functional simplification across fragmented landscapes. Overall, our findings reinforce the importance of landscape configuration for maintaining both taxonomic and functional diversity in heterogeneous mountain landscapes.