We conducted the study in Guatemala, Central America. The sampling was performed in 8 pairs of sites located in the departments of Quetzaltenango, Suchitepéquez, Sololá, Sacatepéquez, Guatemala, Santa Rosa, Jalapa, and El Progreso, one pair in each department, within an elevation range of 748 to 1802 masl ( Supplementary Table S1 ). In addition, the study area is embedded into two life zones within the pre-mountain stratum, sensu Holdridge criteria: Tropical pre-mountain humid forest, with an annual precipitation range of 1000 to 3125 mm, and temperature range of 18 to 24°C ( Supplementary Table S1 ), and Tropical pre-mountain very humid forest, with an annual precipitation range of 2000 to 4850 mm, and a temperature range of 18 to 24°C (same range of the former life zone) ( Supplementary Table S1 ) (IARNA-URL, 2018).

Throughout this region, the species Coffea arabica L. is cultivated and harvested in shaded conditions, under the canopy of forest remnants and plantations, which are mostly comprised of native Inga spp. Mill (Fabaceae) trees and the exotic Grevillea robusta A.Cunn. ex R.Br. (Proteaceae) (ANACAFE, 2011). In addition, coffee is usually cultivated in association with other crops, such as rambutan Nephelium lappaceum L. (Sapindaceae), mangosteen Garcinia mangostana L. (Clusiaceae), tea Camellia sinensis (L.) Kuntze (Theaceae), and macadamia Macadamia integrifolia Maiden & Betche (Proteaceae), where coffee is dominant. Semi-natural vegetation is characterised mostly by young trees, bushes, and weeds and areas of bare soil, which primarily constitute non-cultivated and poorly vegetated areas. Human settlements are also present around the zones dedicated to coffee cultivation. The removal of small bushes and herbs is a common practice across all land covers, especially in semi-natural and human settlement areas neighbouring coffee crops (Casiá-Ajché et al., 2023). Altogether, crop growing areas (including coffee and other associated plantations), forest, semi-natural vegetation areas, and human settlement areas are the principal land-uses in our study sites.

Our fieldwork was carried out in 8 pairs of sites in semi-natural areas, each consisting of two different managed sites (one agroecological and one conventional), with a distance of at least 1 km from each other ( Figure 1 ; information and coordinates of sites in Supplementary Table S1 ). “Agroecological management” refers to the management in which the nearby crops have higher proportions of tree shadow cover, low agrochemical use, use of organic fertilisers, high crop biodiversity and more traditional sustainable practices, while “conventional management” refers to a management where coffee crops have lower tree shadow cover, rely on agrochemicals, have lower crop biodiversity, and less use of traditional sustainable practices. For more details see Escobar-González et al. (2024).

At each site, we recorded bee flower visitors in semi-natural areas across the 8 pairs of sites, recording the identity and frequency of bee-visited plants, and their interactions. Based on these observations, we calculated the composition and diversity of both bee and plant communities and constructed their interaction networks. Fieldwork was carried out during the rainy season, from the final part of March 2018 (onset period of the first rains) to October 2021 ( Supplementary Table S2 ). Intra-seasonal periods were defined as high-precipitation peaks (period A: March to June) and (period B: July to October), separated by a short absent or low precipitation period (see Supplementary Figure S1 ; Supplementary Table S3 ).

We classified land uses into five categories: (i) forest, (ii) semi-natural vegetation, (iii) crops, (iv) bare soil, and (v) human settlements ( Supplementary Table S5 ). Since satellite images for each specific year were not available, the closest available images were used ( Supplementary Table S2 ). Using the coordinates of the sampling sites, we calculated the proportion and edge density of each land use within radii of 1000, 700, 500, and 250 meters, surrounding each sampling sites, for the corresponding years ( Supplementary Figure S2 ; Supplementary Table S4 ). This was done using Google Earth images version 7.3.4.8642 (Google Earth/Maxar Technologies, 2020) with a resolution of 15 meters, the supervised classification algorithm from the ArcGIS software (ESRI, 2023), and the library landscape metrics 1.5.5 (Hesselbarth et al., 2019) from RStudio (R Core Team, 2022). To define the spatial scale at which bee and plant composition, diversity, and their interaction networks respond to the landscape variables, we evaluated the correlation between the landscape variables and the network and community metrics using Pearson correlation, as in Landaverde-González et al. (2017). The spatial scale with the highest correlation for all metrics was 700 meters, so we used this scale in the subsequent analyses ( Supplementary Table S5 ).

Bees and their interactions with flowering plants were collected through transect walks in semi-natural vegetation areas surrounding coffee plantations. The observations were exhaustive and conducted over 2.5 hours per site (total of 40 hours), covering 5-minute observation on each of 30 different plants, per sampling event at each locality. To ensure consistent sampling effort across sites, we used 30 individuals—a statistically recommended sample size (Zar, 1999; Mendenhall and Sincich, 2016)—regardless of the species of flowering plants observed, accounting for the diversity of plants around plantations. Bee visits were recorded for individuals that touched or searched for the reproductive parts of the observed flower. Bee visitors and herbarium specimens of observed plants were collected when these could not be identified in situ, for later taxonomic identification. Samplings were performed during the high activity hours for flower visitors (9am–1pm), and during sunny weather conditions. We conducted a total of 42 samplings –transects- (mean transects per locality = 2.8; SD = ± 0. 34) across 8 pairs of sites from 2018 to 2021 ( Supplementary Table S2 ). For each sampling –transect-, we calculated the diversity and composition metrics of both bee and plant communities, as well as the metrics of their interaction networks.

Both bee and plant individuals were identified to the lowest taxonomic level (genus or species level), using the available identification keys for the region. Bees: Michener et al. (1994); Ayala (1999); Ayala and Griswold (2012); Ayala-Barajas (2016); Balboa et al. (2017); Mérida-Rivas et al. (2022). Plants: Standley and Steyermark (1946–1977). Specimens that could not be identified to species level were defined as morphospecies and classified using morphological criteria. To determine the diversity and composition of both bee and plant communities and to analyse their response to land use types, coffee crop management, and intra-seasonal climatic variation, we calculated: 1) diversity metrics: visit frequency, Chao1 corrected richness (Chao and Shen, 2003), and Shannon diversity index (Shannon and Weaver, 1949); and 2) composition metric: local contribution of beta-diversity (LCBD), based on Sorensen dissimilarity index (Legendre and De Cáceres, 2013). Large LCBD values indicate sites that have strongly different species compositions compared with sites with lower or mean values (Legendre, 2014), which also indicates the ecological uniqueness of sites according to their contribution to global beta-diversity. The diversity metrics were calculated using the Past statistics software, version 4.09 (Hammer et al., 2001), while the LCBD values were calculated using the library adespatial v0.3-8 (Dray et al., 2012), in RStudio software (R Core Team, 2022).

We constructed species/morphospecies-based interaction networks in which bees were considered the higher trophic level (81 taxa identified to species level; 62 taxa identified to morphospecies), and plants, the lower trophic level (160 taxa identified to species level; 15 taxa identified to morphospecies). From each matrix, we computed the metrics: specialisation (H2′), nestedness (NODF), modularity (Q), and connectance (C), described in detail in Escobedo-Kenefic et al. (2020, 2022). We also compute network generality for both trophic levels (bees and plants). Bee generality quantifies the weighted mean effective number of high-level species (bees) per lower-level species (plants) (Bersier et al., 2002). Conversely, plant generality quantifies the weighted mean effective number of lower-level species (plants) per higher-level species (bees); this metric was calculated following vulnerability by Bersier et al. (2002). The network from La Montañita locality in El Progreso department, northeastern of Guatemala (Site code 4B, network 26, Conventional management, and Tropical pre-mountain humid forest, visited in September 2019) was too small to compute the metrics, so we excluded it from the following analyses ( Supplementary Table S6 ).

Due to the dependence of network metrics on network dimensions, we performed delta (Δ)-transformations of the metrics. On the basis of 1000 random networks, Δ-transformed metrics were calculated as Net _(obs) −Net _(rand-mean) in which Net _(obs) is the raw value of the metric and Net _(rand-mean) is the mean value of 1,000 random networks, as presented in Escobedo-Kenefic et al. (2020); Landaverde-González et al. (2021), and Escobedo-Kenefic et al. (2022). All the analyses were performed using the bipartite package (Dormann et al., 2009) in RStudio (R Core Team, 2022). In addition, we used t-tests to evaluate whether Δ-transformed metrics are consistently different than zero; Supplementary Table S7 .

Regarding our question about analysing the effects of land use and coffee crop management on the composition and diversity of both bee and plant communities, and their interaction networks, we used linear mixed models (LMMs) and generalised linear mixed models (GLMs). The composition metric (LCBD values) and diversity metrics (visit frequency, Chao1 corrected richness and Shannon diversity index) of bee and plant communities, and bee-plant network descriptors (Δ-specialisation (H2), Δ-nestedness, Δ-modularity, Δ-connectance, Δ-bee generality, Δ-plant generality), were defined as response variables. To assess the influence of land use categories and their edge density on the response variables, landscape variables were defined as fixed factors. To evaluate the influence of crop management on the response variables, a categorical variable with two levels for coffee management—agroecological and conventional—was defined as a fixed factor; and to measure the impact of plant diversity metrics on bees, plant community metrics were set also as fixed factors. We set site (n = 8) and sampling year (2018 to 2021) as random factors in all the models used to assess our first question.

To address our second question about analysing the effects of intra-seasonal variation on response variables, we used linear mixed models (LMMs) and generalised mixed models (GLMs). Precipitation period, a categorical variable with two levels for each of the precipitation periods A and B, was defined as a fixed factor. For these last models, we set site (n = 8) and sampling year (2018 to 2021) as random factors.

Using the fitdistrplus package 1.2-1 (Delignette-Muller and Dutang, 2015), we confirmed that our continuous response variables had a normal distribution, except for Δ-plant generality. We then normalised this variable using the Log₁₀ (x) transformation. Finally, we tested all model residuals for normality using a Shapiro-Wilk test, using the Past statistics software, version 4.09 (Hammer et al., 2001). The other model assumptions (linearity and homogeneity of variances) were verified visually. These models were run following the Gaussian distribution and using the nlme package 3.1-131.1 (Pinheiro et al., 2018), in R 2.3-1 (R Core Team, 2022). Additionally, visit frequency of bees adjusted to Poisson distribution and Chao1 corrected richness of bees followed a negative binomial distribution. This was verified using the JMP statistical software (SAS Institute Inc, 2012). Both variables were modelled using their respectively fitted distribution using the lme4 package 1.1-35.5 (Bates et al., 2015) in R 2.3-1 (R Core Team, 2022).

In models with landscape variables as fixed factors, we did a model selection by evaluating the effects of land use variables and their edge density on the response variables in separate models, by applying the Akaike Information Criterion (AIC), and its second order value (AICc), in the AICmodavg library 2.3-1 (Mazerolle, 2020). We repeated the analyses with the selected variables, with AICc values < 2, and reported the most parsimonious models (AICc values < 2), as in Escobedo-Kenefic et al. (2022) ( Supplementary Table S8 ).

Additionally, we compared the composition of bee and plant communities between two intra-seasonal precipitation periods, using a paired permutational multivariate analysis of variance (PERMANOVA) using the Wd*-test available in the MicEco package (Russel, 2021). We used non-metric multidimensional scaling (NMDS) to visualise the variation of composition of bee and plant communities between intra-seasonal precipitation periods. Bray-Curtis’s distance matrix was the distance metric with the best fit (stress = 0.186 for bee communities; stress= 0.135 for plant communities) and was the response variable in the PERMANOVA analysis. In addition, precipitation period, a categorical variable with two levels for each of the precipitation periods A and B, was set as the independent variable. The strata (block) argument was set to site, so randomisations were constrained to occur within each site and not across all the sites.

We sampled 1702 bees (more details in Supplementary Figure S3 ; Supplementary Table S9 ) and identified 81 species and 62 morphospecies of bees from 60 genera. Overall, 1151 specimens were Apidae (67.63%), 343 Halictidae (20.15%), 139 Megachilidae (8.17%), 54 Andrenidae (3.17%) and 15 Colletidae (0.88%). The most frequent species recorded throughout the study sites was the stingless bee Trigona fulviventris Guérin-Méneville, 1845 (14.33%). The remaining taxa did not exceed 6% of the total sampled specimens ( Supplementary Table S3 ). Bee taxa visited 160 species and 15 morphospecies of flowering plants. Bees visited more frequently the following plants: 194 visits registered in Bidens pilosa L. (Asteraceae) (11%) and 93 in Sida acuta Burm.f. (Malvaceae) (5%). The rest of the taxa did not exceed 5% of their frequencies of the total plants ( Supplementary Figure S3 ; Supplementary Table S10 ). We observed 1702 interactions within bee–plant networks, with Asteraceae family and stingless bees being the most linked taxa ( Supplementary Figure S4 ). Among our study sites, the highest values of local contribution to beta-diversity of both bees and plants were found in: Finca Buenaventura (conventional management and Tropical pre-mountain humid forest) (Site code 8B, LCBD _bees = 0.15; LCBD _plants = 0.18); Finca El Bosque (Site code 3B, LCBD _bees = 0.15; LCBD _plants = 0.11) and Finca Panamá (Site code 2B, LCBD _bees = 0.16; LCBD _plants = 0.11) (both localities with conventional management in Tropical pre-mountain very humid forest). ( Supplementary Figures S7 , S8 ).

The proportion of land use types (forest, semi-natural vegetation, crop, bare soil, and human settlement covers) and their edge density did not exhibit a significant relationship with the local contribution to beta-diversity (LCBD), visit frequency, Chao1 corrected richness, and Shannon diversity of bees (P > 0.05). For the plant communities, the local contribution to beta-diversity (LCBD), visit frequency, Chao1 corrected richness, and Shannon diversity index were not significantly affected by the proportion of land uses or their edge density (P > 0.05). None of the bee and plant diversity and composition metrics had any significant relationship with the management systems (P > 0.05). None of the plant diversity metrics showed a significant relationship with the bee diversity metrics (P > 0.05).

Regarding bee–plant interaction networks, there were no significant relationships between the proportion of land uses or their edge density with Δ-connectance, Δ-specialisation (H2’), Δ-nestedness, Δ-modularity, and Δ-plant generality (P > 0.05). We observed that Δ-bee generality was not related to any of the edge density of any land use, but concerning the proportion of land use, it did show a marginal negative relationship with human settlement cover (LMM; t = −1.87; P = 0.08; Table 1 ; Figure 2 ). Additionally, Δ-bee generality (across all the sites) was significantly higher than 0 (2.40 mean ± 0.65 SD, P < 0.05) ( Supplementary Table S7 ; Supplementary Figure S10E ). None of the Δ-network metrics were significantly different between crop management systems (P > 0.05) ( Supplementary Figure S11 ).

Neither plant local contribution to beta-diversity nor diversity showed changes between periods of intra-seasonal climatic variation (P > 0.05). The visit frequency and Chao1 corrected richness of bees were significantly different between the intra-seasonal precipitation periods; both descriptors were higher in the second part of the rainy season, encompassing the months of July to October (GLM; z = 0.26; P < 0.001; Table 2 ; Figure 3 ; and GLM; z = 3.17; P = <0.01; Table 2 ; Figure 3 ; respectively).

Plant community composition did not vary between the periods of intra-seasonal climatic variation (PERMANOVA, P > 0.05); Supplementary Figure S9 . However, bee community differed between the intra-seasonal periods (PERMANOVA, W = 2.30; P < 0.01; Figure 4 ).

The network is more generalist during the second period of the rainy season, with bees visiting a more diverse range of plants. However, in both periods, the Asteraceae family remains the most visited plant taxa ( Supplementary Figures S5 , S6 ). Similarly, Δ-plant generality was the only metric of interaction networks that showed differences among intra-seasonal precipitation periods; Δ-plant-generality was significantly higher in the second part of the rainy season (t = 3.27; P < 0.01; Table 2 ; Figure 3 ). Additionally, Δ-plant generality was significantly higher than 0 in both intra-seasonal precipitation periods: Δ-plant generality period A (March to June), 3.77 mean ± 1.39 SD, P < 0.05; Δ-plant generality: period B (July to October); 5.31 mean ± 1.82 SD, P < 0,05) ( Supplementary Table S6 ; Supplementary Figure S12F ).

The number of individuals and species richness of flower visitors recorded in our study are lower than those observed in previous studies conducted in semi-natural areas and coffee crops, with more forest cover in Sololá department in Guatemala (Armas-Quiñonez et al., 2020) and Costa Rica highlands (Cristóbal-Perez et al., 2024). However, unlike our study, both studies conducted monthly sampling, which may explain the greater number of species they observed. In contrast, we observed greater richness of species and abundance compared to other studies conducted in tropical areas. This included studies in coffee crops, in the same sites of this work (Escobar-González et al., 2024), semi-natural areas in the highlands of Guatemala (Escobedo-Kenefic et al., 2020), semi-natural areas and Cucurbita pepo L. crops in Guatemala and Costa Rica (Enríquez et al., 2015; Lobo and Bravo-Méndez, 2021), and semi-natural areas and habanero chili crops in Yucatán, Mexico. In our study, the majority of species belonged to the Apidae family and the second most abundant was the Halictidae family, consistent with other studies in the region (Enríquez et al., 2015; Landaverde-González et al., 2017; Escobedo-Kenefic et al., 2020, 2022; Armas-Quiñonez et al., 2020; Lobo and Bravo-Méndez, 2021; Cristóbal-Perez et al., 2024), indicating a consistent bee family composition across the region.

Our study shows a similar composition of bee-visited plants as reported by Escobedo-Kenefic et al. (2020, 2022) and Flórez-Gómez et al. (2020), where the Asteraceae family was the most visited. In our study, Malvaceae was the second most bee visited family. In contrast, Escobedo-Kenefic et al. (2020) reported Brassicaceae as the second most visited family, while Flórez-Gómez et al. (2020) reported Fabaceae as the second most visited during the rainy season. However, during the dry season, Flórez-Gómez et al. (2020) found that the Fabaceae family was the most visited. This variation seems to depend on the ecosystem and the type of semi-natural area surrounding the study sites.

We observed no significant difference in community composition of bee and plant communities in areas surrounding coffee plantations with different coffee management. Previous works have found no changes in composition of bee communities (Brosi et al., 2007; Serralta-Batun et al., 2024) and plant communities (Dominschek et al., 2021; Serralta-Batun et al., 2024) between semi-natural vegetation remanent and agricultural settings, suggesting that the influence of the management of nearby crops does not cause a change in bee and plant communities. Similarly, crop management did not cause a change in floral visitors on coffee at the same sites of this work (Escobar-González et al., 2024). This pattern may be because many coffee plantations are embedded within Guatemala’s national protected areas system, which includes extensive forested and seminatural areas ( Supplementary Table S4 ; Supplementary Figure S2 ). These natural habitats likely moderate the effects of different management practices, maintaining relatively low-intensity land use across the landscape. While there are differences in management intensity, they may not be extreme enough to significantly impact bee and plant communities of the surrounding areas, as discussed by Escobar-González et al. (2024).

Furthermore, the diversity and composition of bees from vegetation surrounding coffee plantations were not significantly related to any land use variable. Previous investigations have showed that bee visitation frequency and diversity can be affected by local factors such as the abundance and richness of flowering plants surrounding tropical agricultural lands (Landaverde-González et al., 2017; Laha et al., 2020; Escobedo-Kenefic et al., 2020, 2022; Assunção et al., 2022). Semi-natural and agricultural lands have been also found to reduce bee taxonomic and phylogenetic diversity by 16% (Tsang et al., 2024). In addition, we did not observe any relationship between the visit frequency and richness of flowers and bees in the areas surrounding coffee, which might have been caused by a transient dilution of pollinators, likely due to increased competition for pollinators between mass flowering crops such as coffee and wild plants in the surrounding areas (Holzschuh et al., 2016; Grab et al., 2017). We also observed three sites with large values of LCBD for both bee and plant communities, Supplementary Figures S7 , S8 , indicating that these sites had high species turnover, and were ecologically unique and distinct from the others. This may suggest that they are degraded and species-poor, as discussed by Legendre (2014). Interestingly, all three of these sites are conventional farms with a higher degree of human settlements compared to the other farms as shown in Supplementary Figure S2 , suggesting that the type of management and human intervention, at least for these farms, has an impact on the surrounding vegetation, and consequently on pollinators. Intentional removal of wild weeds is a common practice in our study sites, especially in human-dominated land uses such as crop and human settlement cover (Casiá-Ajché et al., 2023). This translates directly into reduced plant diversity available to bees, which may have negative consequences for bee populations surrounding coffee plantations (Machado et al., 2021). However, the extensive natural areas around these farms appear to buffer the unique ecological characteristics of these sites when viewed within the broader landscape context. Previous studies indicate that urbanisation decreases plant richness, with animal-dependent species being more susceptible than those species dispersed by abiotic variables, particularly in tropical regions (McKinney, 2002; Kolbe et al., 2016; Theodorou et al., 2020; Ruas et al., 2022). However, studies in the Neotropics thus far have been limited, and insufficient to fully comprehend which urban drivers, and how the complex relationship between different drivers may have a greater impact on plant diversity and pollinators in tropical urban areas (Ruas et al., 2022; Solís-Montero et al., 2023). In this regard, future studies need to consider the impact of each variable within tropical urban environments on the diversity of plants and pollinators.

We observed that bee generality decreased with the increase of human settlement cover. Our results align with Amado De Santis et al. (2023) in Argentina and Suni et al. (2022) in California, who reported an increase in specialisation in urban areas, possibly due to the presence of invasive plants, which may be leading to changes in pollinator fidelity for exotic plants in urban areas. Opposite effects have been found in temperate regions, where the heterogeneity of intermediate succession in disturbed semi-natural areas led to the formation of robust ecological networks, with consistent generality despite the varying levels of habitat loss (McWilliams et al., 2019; Deprá et al., 2022). Silva et al. (2020) showed that urbanisation may pose a significant threat to pollinator–plant networks in tropical regions, however, these ecosystems remain understudied, warranting further research. Furthermore, bee generality exceeded 0 across all sites, indicating that bees are more generalised than expected by chance. This suggests that although human settlements tend to reduce the number of flower visitors visiting plants, the system remains generalist, possibly due to the availability of diverse floral resources in the surrounding semi-natural and forested areas. Considering that our sampling areas predominantly encompass agricultural regions, with high forest and semi-natural cover, and with urban areas comprising about 5% of the landscape, our findings suggest that while urbanisation does influence bee generality, it may be present in low enough levels to not significantly disrupt the overall stability of the interaction network. The strong representation of semi-natural and forested areas in our study, as shown in Supplementary Figure S2 , likely explains the lack of variability in other network metrics, as these areas may act as a buffer against the effects of crop management on diversity of bees and their interactions with plants. This underscores the robustness of the structure of bee-plant networks at the current intensity of land use and human disturbance from settlements.

We observed a different composition of bee flower visitors during the two periods within the rainy season. Furthermore, visit frequency and richness of bees showed a significant increase in the second period of the rainy season (i.e. July to October). It has been outlined that tropical plant communities often synchronise their flowering with the peak of seasonal rainfall (Lasky et al., 2016). In our study sites, the period from July to October marks the second annual peak of rainfall. The effect of precipitation period could have been caused by precipitation-driven blooming of plant communities (Lasky et al., 2016; Cortés-Flores et al., 2017, 2023). Similarly, several studies have shown that pollinator diversity varies according to floral resource abundance. In natural systems, a strong correlation between pollinators and their services with flower abundance has frequently been observed, which favours pollination during early flowering, although pollinator abundance increases as flowering progresses (Elzinga et al., 2007). Other studies have demonstrated that bee and wasp diversity correlate with the number of plant species and flowers (Abrahamczyk et al., 2011). This correlation is suggested to be due to the availability of flowers and nectar, which increase with rainfall (Ballarin et al., 2022). Other studies observed a significant influence of floral richness and precipitation on both abundance and richness of pollinators. However, despite these influences, the interaction networks remained stable across months and years (Stewart and Waitayachart, 2020) or across varying levels of habitat disturbance (McWilliams et al., 2019; Deprá et al., 2022). Similarly, we observed that most of the network metrics did not vary between the precipitation periods, except for plant generality, which increased during the second period of precipitation. This indicates that with the increase in rainfall during the second period, the number of plant species (lower trophic level) visited by each bee species significantly increases compared to the increase in the number of bee species (high trophic level). This effect is notably more pronounced in the second precipitation period, where heightened precipitation likely supports an increase in the blooming of plant communities. Furthermore, plant generality exceeded 0 in both intra-seasonal precipitation periods, indicating that plants are more generalised than expected by chance during both periods. Altogether the increase in visit frequency and richness of bees may be specifically linked to rainfall-triggered blooming that positively increases plant generality, resulting in a more robust network during the second part of the rainy season.