We considered the Brazilian territory in the study as its 8.5 million km² makes it the largest country in Latin America (5th in the world), ranging from various climatic zones from the humid tropics in the north to temperate areas in the south. Those climatic differences led to the formation of the following distinct biogeographic regions: the Amazon and Atlantic Forest (rainforest), the Pantanal (floodplain), the Cerrado (Neotropical savannah), the Caatinga (tropical dry forest), and the Pampa (natural grasslands) (Figure 1).

We classified the species within each guild according to the body size classes of small, medium, and large bats. Our approach used mean forearm length as a proxy for body size, which is an adequate classification method since forearm measures are precise and associated with bat’s body length and mass content (Meng et al., 2016; Thiagavel et al., 2017). We gathered the morphological data from published literature (Supplementary Table 2). With a kmeans clustering analysis, we identified intervals of mean forearm lengths that represented small, medium, and large categories for each guild separately. These were then attributed to the values of 1, 2, and 3, respectively.

The hierarchical classification expresses the decrease in the number of species for a potential replacement with increased body size. Further, it also expresses the broader functional range of larger species, which are less limited by morphology in performing their ecological functions (Fleming, 1991; Kalko et al., 1996; Houston et al., 2004). Thus, we attributed greater value to larger species due to their rarity and potential functional uniqueness within their guilds.

We considered all species indicated by the modeling process in each biome separately for this stage. We analyzed the potential contribution of bats to ecosystem service in each major Brazilian biomes, which were divided in a regular grid of cells with a resolution of 0.5 × 0.5 decimal degrees (Figure 1). Within each grid cell we calculated a diversity index of provider richness (hereafter Index of Potential Ecosystem Services Provision or ES index) like the one applied in Ceauşu et al., 2021. This was done by summing the area occupied by a species within the same guild multiplied by their class size (i.e., 1, 2, or 3). The grid cell values were normalized by dividing each cell value by the maximum value found among cells for each biome separately. This resulted in values ranging from 0 to 1 representative of the cell’s importance for providing of each guild’s associated services. We repeated this analysis using all bats independent of the guild in each biome and considered the scenarios with and without natural habitat suppression.

We wanted to test how habitat conversion has affected the potential provision of ecosystem services by bat species in Brazil, identifying regions with different loss of ecosystem service provision.

For the representation of land conversion, we used the MapBiomas version 5 land cover product,³ representing the year 2020. The map is presented at a 30 m resolution and features several landcover classes. We resampled it to the same spatial resolution as our models (ca. 2.5 arc min resolution). We created a binary map of natural vegetation remnants, which included all classes of forest and non-forest natural formations as remnants. All areas occupied by pastures, agriculture, urban infrastructure, and water bodies became zero-valued. Thus, this mask was used to suppress suitable areas within the species’ potential distribution. With it, we could compare the potential for ES provision prior to and after natural area suppression. With this approach, we assume that diversity, its functions, and associated services depend on the existence of natural remnants. This means we followed the notion of “supporting systems” suggested in Kremen (2005) and Tscharntke et al. (2005), i.e., any diversity present in altered landscapes depends upon the support of the surrounding native landscape. All spatial analyses were conducted in the R software (R Core Team, 2021), using the packages raster (Hijmans, 2021) and sf (Pebesma, 2018).

We were interested in identifying which features of diversity work to maintain the ecosystem service provision within grids where natural habitat suppression has occurred. We applied generalized linear models (GLM) with Gaussian distribution using the ES index values after deforestation as our response variable and remnant area occupied by species plus the number of lost species of the same size class as our explanatory variables. We only consider insectivores, frugivores, and nectarivores for this stage because of their number of species. This was done separately for each guild within each biome separately.

We generated 128 generated models with good predictive performance. TSS values ranged from 0.6 to 0.97, with 53 species validated with scores between 0.6 and 0.69 and 60 species validated between 0.7 and 0.97. Five species were removed from the study due to unsatisfactory model results after a visual inspection. Thus, all results presented here are based on 128 bat species, of which 76 were identified as insectivores, 36 as frugivores, 12 as nectarivores, and four as carnivores (Supplementary Table 2).

We calculated for each trophic guild three mean forearm lengths intervals (small, medium, and large) using kmeans clustering (Table 1). Overall, there were 63 small, 38 medium, and 27 large bat species. Size differences were most marked for insectivorous and frugivorous bats (Table 1) and after specialist consultation we decided to keep only small and medium classes for nectarivores and only medium and large for carnivores. Overall, most species within each guild were classified as small (Supplementary Table 2).

For Brazil and for each biome, values of the corresponding potential for the provision of ES were identified (i.e., disregarding natural area suppression). Overall, the values of the ES index varied among biomes, indicating differences in species compositions between regions (Table 2). Cerrado was the only biome that presented mean ES index values greater than the national average for all guilds (Table 2). The potential for insectivores’ associated services is predicted to be greater in the Pantanal, Caatinga, and Cerrado biomes (Table 2). Nectarivores species presented higher values for the Amazon and Caatinga, on average. Frugivory seems to be well distributed across Brazil, with the weakest potentials predicted for Caatinga and Pampa (Table 2). Finally, carnivory potential is predicted to be distributed toward Pantanal, Atlantic Forest, and Caatinga (Table 2).

Knowing the potential provision scores for a country or region is not the only way of assessing a place’s likelihood of providing or sustaining specific ecosystem services. A complementary approach is to see how these values are distributed in space. The Cerrado, for instance, has 60.99% of its original cells, i.e., without considering area suppression, presenting ES index values greater than the national average. The other biomes have 55.52% (Caatinga), 54.44% (Pantanal), 52.83% (Atlantic Forest), 50.81% (Amazon), and 44.44% (Pampa) (Table 3).

With the historical occupation of Brazil, biomes’ ES provisioning potential eroded. The Atlantic Forest biome is an extraordinary case, with a total loss of more than 70% of its original area (ca.1.1 million km²). Another striking case is the one of Cerrado, which had lost more than 1 million km² of its native area. Considering a threshold of 50% of remaining native areas, the percentage of cells below such limit is 74.53% for Atlantic Forest, 59.39% for Pampa, 55.65% for Cerrado, 27.84% for Pantanal, 27.62% for Caatinga, and 16.05% for Amazon.

The land conversion process observed up to 2020 has impacted the provision of ES differently in each biome and, on average, the Cerrado, Atlantic Forest, and Caatinga, were the biomes with the highest losses (Table 3 and Figure 2—Mapa de perda de ES). Most of central-southern Cerrado, central Atlantic Forest, and eastern Caatinga had heavily lost the provision potential of ES. When comparing the absolute loss of ES potential, i.e., the values resulting from the difference of the original provision by the remaining provision after natural area suppression, Atlantic Forest is the biome with the highest loss, followed by Cerrado, Caatinga, and Pampa, and Amazonia (Kruskal-Wallis chi2 = 685.81, p < 0.001) (Figure 3). The proportion of cells with expressive loss of ES is 62.26% for Atlantic Forest, 50.77% for Cerrado, 44.44% for Pampa, 33.72% for Caatinga, 24.05% for Pantanal, and 16.31% for Amazonia. Considering the cells with ES index equals or greater than the national mean, only 28.3% of essential areas for ecosystem services provision stand in Brazil (Figure 4). This proportion changes among the biomes, being 7.9% for Atlantic Forest, 16.7% for Pampa, 19.6% for Cerrado, 24.5% for Caatinga, 39.8% for Amazonia, and 50.1% for Pantanal.

The negative impacts of natural area conversion in Brazil varied among guilds and regionally. The general pattern for Brazil is a decrescent and significant impact of devegetation on insectivores, frugivores, carnivores, and nectarivores (Kruskal-Wallis chi2 = 391.49, p < 0.001) (Table 4). However, the pattern of impacts on different guilds varied per biome. The frugivores were more impacted in the Amazon, while insectivores and frugivores were more impacted in the Cerrado and Pantanal. In Atlantic Forest, there was no statistical difference between carnivores, frugivores, and insectivores, but these three categories were more impacted than nectarivores (Table 4). The Pampa was the unique biome where the devegetation impact affected equally all guild categories (Table 4).

Considering all feeding guilds, the loss of ES provided by large species was significantly higher than the other two size categories (Yemeni test, p < 0.05), but there was no statistical difference between medium and small size species (Yemeni test, p > 0.05) (Table 5). This result partially corroborated our working hypothesis and was regionally observed in Amazon and Cerrado. We observed, as expected, the loss of bat ecosystem services in the Atlantic Forest, higher in large species than in medium and smaller-size species (Table 5). Furthermore, we did not detect a significant loss of size groups ES for the Pantanal and Pampa biomes (Table 5).

Considering each feeding guild separately, we observed different responses for each biome. For carnivores, the loss of ES per size group was not significantly different for Brazil, Amazon, Pantanal, and Pampa. We observed the expected trend for Caatinga and Atlantic Forest (Table 5). For the Cerrado, we observed the expected trend for all groups except for carnivores and nectarivores (Table 5). For Brazil, Cerrado, and Atlantic Forest, we observed the loss of ES per size group according to expectations for frugivores. The expected trend was observed only for insectivores in the Atlantic Forest and for most other biomes, we found only a partial corroboration. Finally, nectarivores were the only group with a consistent response but presented the opposite pattern than expected (i.e., small species being more impacted than medium-size species).

Our results indicated that reducing area of occurrence and loss of species impacted the provision of ES by bats in Brazil. However, such changes vary among biomes (Table 6). For Atlantic Forest and Cerrado, both factors significantly reduced the ES index values for multiple groups. For instance, for frugivores and nectarivores, the reduction of species richness and occurrence area had significant relationship with the loss of ES in Cerrado (Medium frugivores species loss: t = −4.808, p < 0.001; Medium frugivores occurrence area reduction: t = 7.674, p < 0.001) (Table 6). A similar pattern was found for the Atlantic Forest, where loss of species and occurrence area reduction were relevant, but only for small nectarivores (Small nectarivores species loss: t = 2.494, p < 0.05; Small nectarivores area reduction: t = 7.121, p < 0.05). For the other bat groups only one of the factors significantly reduced the ES index values. For example, small insectivores associated services were affected only by the loss of species, whereas medium insectivores, only by the reduction of occurrence area (Small insectivores: t = −3.535, p < 0.05; Medium insectivores: t = 7.572, p < 0.05). For the biomes of Amazon, Caatinga, and Pampa only the reduction of occurrence area was significantly associated with the loss of ES index values for all guilds and sizes (Table 6).

ES provisioning potential is being lost across Brazil through two main processes. First, by guild occurrence area reduction where land conversion shrinks the available suitable area for species within a guild. This habitat area reduction likely culminates in population declines, negative impacting ES provision numerically (Kremen, 2005; Watling et al., 2020). Second, by the loss of species. Some species may lose all their previously available area to land conversion and be considered locally extinct. This reduces the number of provider species present and adds yet another layer of ES vulnerability, the loss of redundancy within guilds (Elmqvist et al., 2003). One can imagine that ES can be sustained in low provider richness situations if the populations of the few species present are healthy enough to maintain a stable provision (Kremen, 2005). However, habitat loss affects both species population and communities’ composition, likely affecting ES provision and stability synergically (Oliver et al., 2015).

Both processes (i.e., occurrence area reduction and species loss) were not always relevant for all groups within all biomes. Meaning that some groups’ associated services were more impacted by area reduction, species loss or both. This is to be expected since guilds are differently distributed across biomes and biomes have different histories of occupation (Lapola et al., 2014). For example, for biomes where anthropogenic occupation is advanced (e.g., Cerrado and Mata Atlantica), both processes seem relevant, although not for all size groups within guilds. This means that beyond occurrence area reduction, species are being lost and this is also significantly impacting the ES potential of the region. Consequently, these regions are experiencing a strong enough loss of provider richness that outweighs the compensational effects of the remaining species within the assemblage. On the other hand, biomes where anthropogenic occupation is not yet as drastic or follows different spatial patterns (e.g., Amazon and Caatinga), reduction of occurrence area is the main driver of ES index values differences. Thus, area reduction is driving the observed differences in ES values, meaning that the functional impact of species losses is being buffered by the functional effect of remaining species. In the scenario, ES should be better safeguarded through mechanisms of redundancy. Are we losing ecosystem services?

Studies on bat’s provision of ecosystem services have mainly focused on the aggregate function (i.e., the service resulting from the ensemble of providers present, Kremen, 2005). For example, a common practice for the studies on pest suppression is exclosure experiments (Maas et al., 2016). While these studies offer essential insights on local diversity impacts on crop yield and quality (Williams-Guillén et al., 2008; Maas et al., 2013; Maine and Boyles, 2015), little light is shed on the dynamics of provision. Even with the aggregate service measure at hand (e.g., number of fruit sets, plant damage), questions of how much diversity is needed for the desired outcome or which species/traits are associated with the effectiveness of service delivery remain unanswered. As such, these results offer evidence for service existence while making good cases for conservation of bat diversity, accounting for the importance of maintaining multiple species contributions to a service (Elmqvist et al., 2003; Aguiar et al., 2021). In other words, the importance of maintaining high provider richness (Isbell et al., 2011).

Thus, conservation of guild diversity is an essential precautionary approach to ecosystem services persistence in front of uncertainties in the dynamics of provision since it aims to secure function stability (Walker, 1992). As we assumed, function stability implies redundancy, which implies replacement among species (Luck et al., 2009). Through our classification, the potential for redundancy (i.e., the number of species for replacement) decreases as body size increases. This is a pattern observed on mutualistic networks (Zamora, 2000). It means that the local loss of larger species with relatively greater functional importance should come with greater functional consequences due to reduced stability and the trait-limited possibility of compensatory responses from other species. Consequently, this may result in rapid function loss (Kremen, 2005). As reported, we observed such general pattern, but our results cannot ascertain that this is happening because we are not actually measuring ES provision, we calculated a proxy theoretically associated with the potential for ES provision.

Furthermore, in the case of redundancy, the very own capacity of the bat population to respond in a compensatory way should be experimentally investigated. Neotropical bat communities are not consistently structured by competitive interactions (Stevens and Willig, 2000), and bat species are long-lived organisms with low fecundity (Findley, 1993). Thus, bats are likely not responsive to population fluctuations and not capable of the rapid growth such compensatory responses require. Instead, due to their high mobility, function stability might depend on a rescue effect from surrounding populations, which suggests the importance of biodiversity conservation on a landscape scale and the preservation of diversity sources (Tscharntke et al., 2005; Farneda et al., 2018). We argue that areas with above average ES index values should be considered diversity sources for regional ES provision and, as such, should be of conservational focus and restoration actions.

Despite the uncertain role of redundancy in bat function stability, guild diversity can enhance aggregate function through mechanisms of species complementarity (Gaston and Spicer, 2004; Williams-Guillén et al., 2008; Aguiar et al., 2021). Bats exploring similar resources separate their uses in space and time (i.e., niche partitioning) (Arlettaz et al., 1997; Schnitzler and Kalko, 2001; Lopez and Vaughan, 2007; Amaral et al., 2016). Consequently, species performing similar functions at different spatial-temporal scales provide stability to ecosystem services provision (Peterson et al., 1998). This reiterates the importance of conserving the natural levels of diversity of regional bat assemblages. However, in the finer scale of ecological function, species likely differ in their potential contribution to ecosystem services provision, and body size, as we assumed, could be a trait dictating this in bats, as it is on other vertebrates (Luck et al., 2012). If such is proven, it is also a well-established response trait indicating vulnerability to disturbances like habitat fragmentation (Schulze et al., 2000; Farneda et al., 2015) and greater attention directed toward those species.

Identifying potential key provider species and their contributions to aggregate function and the stability of associated ES is crucial and urgent. The current trend of rapid destruction of native habitat may be leading to the disruption of ecosystem functioning and service provision with unknown consequences. The true state of the functional contribution of bat species and assemblages across Brazil to ecosystem service provision is yet to be described. As we stated, our results reflect general patterns, and any attempt to describe the effective distribution of ES provision by bats would be speculative. Many assumptions remain to be tested for our understanding of the dynamic ES provided by bats to be consolidated. For now, we may be losing the benefits provided by these populations without even acknowledging their existence. Although Brazil is one of the countries with the most incredible diversity of bat species (181), this diversity could not be completely represented in our study for a couple of reasons. First, there are not enough records of good enough quality for all species to be modeled. In our study, seventy-eight species (42.6%) could not be modeled due to the lack of records. Secondly, Brazil is large country where only a few areas were well surveyed, remaining many areas with very few surveys. This causes an imbalance on the representation of the environmental space occupied by large ranged species that resulted in badly fit models even for species that had sufficient occurrence records (e.g., Desmodus rotundus, Carollia perspicillata, etc.).

If, on the one hand, the species that we modeled were sufficient to produce a scenario of the potential provision of ecosystem services in Brazil and its biomes, on the other, we detected that there is still a large gap to be filled concerning the occurrence of species in Brazil, as already pointed out by Bernard et al. (2011) and Aguiar et al. (2020). The search for records must be done optimally since the inclusion of new records in little-known regions can substantially alter the distribution modeling results (Aguiar et al., 2015). Furthermore, the study on bat’s contribution should be done at the scale of the service, but country scale descriptions of species distribution must also continue in parallel. This way we can continue to improve our methods for large-scale descriptions of ecosystem service provision and better incorporate this aspect in future conservation plans.