Interplay between local and landscape-scale effects on the taxonomic, functional and phylogenetic diversity of aerial insectivorous neotropical bats

Human-modified landscapes are globally ubiquitous. It is critical to understand how habitat loss and fragmentation impact biodiversity from both a local habitat context and landscape-scale perspective to inform land management and conservation strategies. We used an experimentally fragmented landscape in the Brazilian Amazon to investigate variation in aerial insectivorous bat diversity in response to local habitat and wider landscape characteristics, applying a multiscale approach. We conducted bat acoustic surveys at 33 sites, comprising old secondary forests and fragments of primary forest. Taxonomic, functional and phylogenetic diversity facets were calculated within a Hill numbers framework. We analysed responses to fragment size, interior-edge-matrix gradients, as well as local vegetation structure, continuous forest cover, edge density and patch density across five spatial scales (0.5−3 km) surrounding detector locations. Compared with continuous forest, secondary forest matrix around the smallest fragments harboured lower diversity. The overall negative effect of the matrix became less pronounced with increasing fragment size. In contrast, forest edges generally contained higher taxonomic, functional and phylogenetic diversity. We found subtle scale-sensitive associations for functional diversity, responding positively to forest cover (at the 1 km scale) and negatively to edge (1 km scale) and patch density (2.5 km scale). Despite a low-contrast matrix of tall secondary forest surrounding fragments after ~ 30 years of forest recovery, aerial insectivorous bat diversity is not comparable to continuous primary forest. Assemblage functional diversity responds to compositional and configurational landscape characteristics at scales deserving further evaluation at guild and species level.


Introduction
Globally, over the past 300 years, there has been a net forest loss of ~ 7-11 million km 2 , primarily due to logging and conversion to agricultural use (Foley et al. 2005). Loss of suitable habitat to sustain species populations leads to an overall increased risk of extinction (Ceballos et al. 2015(Ceballos et al. , 2017Powers and Jetz 2015). Forest patches of varying sizes and shapes, embedded in various land cover types, are conspicuous features of present-day human-modified landscapes (Melo et al. 2013). The link between remaining native habitat and the species persisting in the aftermath of habitat fragmentation was first described through the lens of Island Biogeography Theory, which sought to predict species loss on islands surrounded by a 'hostile' matrix of water (MacArthur and Wilson 1967). Yet, mounting evidence now supports that the structure and composition of the intervening matrix significantly affect species' persistence in adjoining forest fragments, a dynamic that is better conceptualised under alternative theoretical frameworks, such as Countryside Biogeography (Daily 1997;Mendenhall et al. 2013). Through this framework, community changes in modified landscapes are forecast based on the interaction between the species' spatial requirements and their tolerance towards matrix habitats, thus offering a more realistic portrait of biodiversity persistence in landscapes with matrix habitats more salubrious than water (Mendenhall et al. 2014;Wolfe et al. 2015;Farneda et al. 2020).
In the Amazon, deforestation rates dropped from 30,000 km 2 /year in the 1980s to 5843 km 2 /year in 2013 (Davidson et al. 2012;Nepstad et al. 2014). However, since 2013, alongside the main driver of deforestation in the Amazon, pasture expansion for cattle production (Skidmore et al. 2021), compounding threats such as oil palm plantations (Butler and Laurance 2009), expanding soy agriculture (Rosa et al. 2017) and dam development (Lees et al. 2016) have contributed to a sharp increase. Much of this deforested land has been abandoned, and a recent study estimates a total of 262,791 km 2 of recovered secondary forests in Brazil between 1986 and 2018 (Silva Junior et al. 2020). Secondary forests make up a significant proportion of fragmented tropical landscapes (Chazdon et al. 2009;Chazdon 2014). A growing body of literature supports that these regenerating forests can alleviate fragmentation impacts and support diverse assemblages and overall ecosystem functioning (Spake et al. 2015;Farneda et al. 2018;Rocha et al. 2018;Rozendaal et al. 2019).
The Neotropics are a major hotspot of chiropteran diversity and the Amazon basin, with over 160 species, is especially rich . Throughout the region, bats play essential roles in countless ecological networks and provide valuable ecosystem services such as seed dispersal, pollination and arthropod population regulation (Aguiar et al. 2021;Ramírez-Fráncel et al. 2022). They are demonstrably sensitive to habitat loss, fragmentation, and habitat degradation (reviewed in Meyer et al. 2016) and have become a popular indicator group of environmental disturbance (Cunto and Bernard 2012). Yet, while the consequences of anthropogenic forest fragmentation for phyllostomid bats have received substantial research attention (e.g., Klingbeil and Willig 2009;Avila-Cabadilla et al. 2014;García-García et al. 2014;Arroyo-Rodríguez et al. 2016;Muylaert et al. 2016;Farneda et al. 2021;, aerial insectivores-non-phyllostomid counterpartshave largely been neglected and their responses to habitat fragmentation remain understudied (but see e.g., Estrada-Villegas et al. 2010;Rodríguez-San Pedro and Simonetti 2015;Núñez et al. 2019).
Measures of diversity based solely on traditional species counts fail to encapsulate the complexities associated with the distinct functional roles of different species or the evolutionary history contained within a given assemblage (Cadotte et al. 2013). A multifaceted approach considering the complementary taxonomic, functional and phylogenetic dimensions of diversity can provide a more detailed and comprehensive understanding of the drivers of biodiversity change across human-modified landscapes (Swenson 2011). However, studies simultaneously assessing how multiple dimensions of Neotropical bat diversity are affected by habitat gradients are still scarce (Cisneros et al. 2014;Frank et al. 2017;Carrasco-Rueda and Loiselle 2020;Carvalho et al. 2021), particularly for aerial insectivores (but see e.g., Pereira et al. 2018).
Amongst the research questions commonly posed by fragmentation studies, the role of fragment size, edge effects, and compositional vs. configurational aspects of the landscape has rarely been addressed for aerial insectivorous bats. Here, we set out to help fill this gap by investigating how local vegetation structure and landscape composition and configuration affect this ensemble along a disturbance gradient formed by continuous primary forest (CF) and primary forest fragment interiors (I), forest edges (E) and secondary forest matrix (M) habitats (hereinafter IEM gradients; sensu Rocha et al. 2017a). Specifically, we aimed to address two objectives: (1) Quantify between-habitat differences in the taxonomic, functional and phylogenetic diversity of aerial insectivorous bats along IEM and fragment-size gradients. We anticipated (i) that assemblages in CF are taxonomically, functionally and phylogenetically most diverse, with diversity being eroded through fragment interiors (< 10 ha) and further still in the matrix, (ii) similar levels of diversity in CF and larger fragments (100 ha), with a reduction in the three biodiversity dimensions with decreasing fragment size (10 and 1 ha), (iii) a positive response of all three diversity facets at the fragment edges as more species are able to take advantage of foraging opportunities along the ecotone.
(2) Assess the importance of the following variables as determinants of assemblage diversity, using a multiscale approach; local vegetation structure, landscape composition (forest cover), and landscape configuration (edge density and patch density). We predicted that (i) taxonomic, functional and phylogenetic diversity would exhibit only subtle responses to the local and landscape characteristics, (ii) local vegetation structure would generally have a weaker effect than landscape metrics and, when present, would manifest at the smallest scale (0.5 km), possibly turning neutral as the scale increases, (iii) forest cover would be the most important predictor of all three diversity facets at the landscape scale, similar to findings for the species richness and abundance of phyllostomids (Rocha et al. 2017a).

Study area
Fieldwork was conducted at the Biological Dynamics of Forest Fragments Project (BDFFP), located ~ 80 km north of Manaus, Central Amazon, Brazil (see Fig. 1). Established in 1979, the BDFFP is the world's largest and longest running experimental study focused on habitat fragmentation, landscape dynamics, forest regeneration, and regional and global changes affecting plant and animal communities (Laurance et al. 2018). The topography of the ~ 1000 km 2 study area is relatively flat, with many small streams cutting through the nutrient-poor soil. Annual rainfall varies from 1900 to 3500 mm with a moderately strong dry season from June to October (Laurance et al. 2011). The area is characterised by a mosaic of primary forest fragments embedded in a second-growth forest matrix surrounded by large areas of continuous terra firme forest. The fragments were first isolated through logging and burning for cattle pasture in the early 1980s. Fragments of different sizes (1 ha, 10 ha and 100 ha) were isolated from CF by distances ranging from 80 to 650 m. Following the abandonment of the cattle ranches, the second-growth forest matured around the fragments over ~ 30 years (Laurance et al. 2018), creating a landscape of low structural fragment-matrix contrast (Fig. 1). A 100 m-wide strip of regrowth vegetation has been cleared on at least five occasions since fragments' creation to maintain their integrity. The last re-isolation event prior to data collection for this study occurred between 1999 and 2001 ).

Bat surveys
Acoustic data were collected between 2011 and 2013 in both the interior and at the edges of eight forest fragments (three of 1 ha, three of 10 ha, two of 100 ha), eight secondary forest sites (located 100 m into the matrix from the nearest fragment edge) and nine control sites in three continuous forest areas.
Recordings were obtained with SM2Bat + detectors with omnidirectional microphones SMX-US (Wildlife Acoustics Inc., Massachusetts, USA). At each site, the detector was positioned at ~ 1.5 m height and programmed to record from 18:00 to 06:00. Recordings were captured in real time with a full spectrum resolution of 16-bit, with a high pass filter set at fs/32 (12 kHz) and a trigger level of 18SNR. Sites were sampled for one night each survey during the first year (2011)  López-Baucells, unpublished data). Calls were identified to species level when possible or assigned to groups of taxa with similar calls (sonotypes). A total of 21 species/sonotypes were identified (Table S1 in Online Supplementary Material). Since it is not possible to estimate abundance with acoustic data, we used activity as a proxy of abundance based on the number of bat passes (Rowse et al. 2016). The sample unit, a bat pass, was defined as any call sequence with a maximum duration of 5 s, which contained at least two distinguishable echolocation pulses López-Baucells et al. 2021). Bat activity was quantified by the total number of bat passes per night per species/sonotype.

Calculation of diversity response metrics
We calculated taxonomic (TD), functional (FD) and phylogenetic (PD) alpha diversity using a unified framework based on Hill numbers. Hill (1973) integrated species richness, the converted Shannon entropy and Gini-Simpson index measures into a family of diversity measures by order q or the effective number of species. The sensitivity of the measure to the relative frequency of species is determined by the parameter q and is expressed in units of species.
Hill numbers offer advantages over standalone diversity measures, including satisfying the mathematical replication principle, allowing for direct comparison across orders q due to values being expressed as units of effective number of species. They can be partitioned into independent group components and by doing so can be generalised to taxonomic, functional, and phylogenetic diversities providing a unified framework for measuring biodiversity ). Using the R package 'hillR' (Li 2018), we calculated each diversity facet based on total bat activity per site and per night using the Hill numbers framework. Diversity values become more sensitive to common species as q increases. When q = 0, species/sonotype abundance is ignored (species richness); q = 1, all species/sonotypes are weighted by their abundance equally (Shannon diversity (the exponential of entropy)); q = 2, greater weight is placed on common than rare species/sonotypes (Simpson diversity (inverse of the Simpson index)). We calculated all three q values (0, 1 and 2), representing a full diversity profile illustrating the species abundance differentiation of the assemblage. For FD, we selected functional traits that have been shown to indicate potential vulnerability to habitat fragmentation in aerial insectivorous bats (Núñez et al. 2019). Trait information on echolocation call structure, frequency and alternation, body mass, aspect ratio and relative wing loading, as well as vertical stratification were considered, encompassing both continuous and categorical data (Table S2 in Online Supplementary Material). To adhere with the Hill number framework, each 'species/sonotype' must be a distinct entity with its own (attribute) diversity data . Complete trait data was not available for Saccopteryx gymnura, thus, calls of this species were removed from the analysis (1817 calls in total). Also, for molossid sonotypes, body mass, aspect ratio and relative wing loading were calculated using the mean of values for individual species within each sonotype. To quantify PD, phylogenetic information was extracted from a species-level supertree for bats (Shi and Rabosky 2015) and pruned to include only the species of aerial insectivorous bats known to occur at the BDFFP ( Fig. S1 in Online Supplementary Material). Again, the 'hillR' methodology required a single entity species to be selected from the supertree, therefore to deal with sonotype data, we selected a single representative species, based on the likelihood of occurrence at the BDFFP, for the PD analysis; Molossidae II-Molossus rufus, Molossidae III-Eumops auripendulus, Molossops I-Molossops neglectus and Promops I-Promops centralis (See Table S3 in the Online Supplementary Material).
Pteronotus alitonus was not present in the supertree and so was replaced by its closest congener, Pteronotus parnellii (Pavan et al. 2018). It has been argued that this approach of substitution by close congeners does not cause serious changes in the distance matrices (Cisneros et al. 2014).
Local and landscape predictor variables  Table S4 in Online Supplementary Material). This first component was positively correlated with average DBH of trees > 10 cm, canopy height and percentage canopy cover. Although these features describe more structurally complex habitats, the metric represented a gradient, including negative values, characteristic of secondary regrowth (matrix), with an increased number of woody stems and liana density.

Landscape structure
Landscape composition and configuration were quantified using 2011 LandSat Thematic Map-per™ satellite images (30 m spatial resolution) with continuous forest and secondary forest land cover classes identified. Collinearity is a common problem with landscape predictor variables. Therefore, we selected the same three landscape metrics as Rocha et al. (2017a), which were acceptable based on variance inflation factor calculations (VIF). Metrics of landscape composition (primary forest cover) and configuration (patch density, edge density) were calculated using the R package 'landscapemetrics' (Hesselbarth et al. 2019). Circular buffers were defined with radii of 0.5, 1, 1.5, 2, 2.5, and 3 km around the 33 sampling sites. In selecting these buffer sizes, we took into consideration the observed scale of effect for bats (Jackson and Fahrig 2015) and  (Meyer and Kalko 2008). At each spatial scale, we first fitted a linear model between the landscape composition variable (primary forest cover) and the landscape configuration variables (edge and patch density) (Trzcinski et al. 1999). The residuals of the simple linear regression were then extracted and new configuration variables were created for use in subsequent modelling (Bélisle et al. 2001;Klingbeil and Willig 2010).

Modelling taxonomic, functional and phylogenetic diversity
The diversity metrics (TD, FD, PD/q = 0, q = 1, q = 2) were modelled (1) with the fragment size/IEM variables (see below) and (2) the local-and landscapescale variables, using Bayesian generalised linear mixed-effects models (GLMMs) as implemented in package MCMCglmm (Hadfield 2010). A measure of significance of the difference between effects (p MCMC ) can be produced via the proportional overlap of the distribution estimates produced, accompanied by an estimate of the mean and 95% credible intervals without post-hoc tests (Sweeny et al. 2021). For fragment size/IEM variables we fit a set of nine models, three for each diversity facet (TD, FD, PD/q = 0, q = 1, q = 2), specifying a single categorical fixed effect with combined information on the fragment size and IEM gradient (10 categories: continuous forest interior, 100 ha interior, 100 ha edge, 100 ha matrix, 10 ha interior, 10 ha edge, 10 ha matrix, 1 ha interior, 1 ha edge, 1 ha matrix; Rocha et al. 2017a) and incorporated research camp location as a random effect.
To model local-and landscape-scale variables, we fit a set of models using each by site alpha diversity metric (TD, FD, PD/q = 0, q = 1, q = 2) in turn, with four fixed effect local and landscape-scale variables (local vegetation structure, continuous forest cover, edge density, patch density) for each buffer size (0.5, 1, 1.5, 2, 2.5 and 3 km) and research camp location as a random effect. Each model set contained six full models and nine sets were run in total. All models were fitted with a Gaussian error distribution and a "non-informative" prior, which is weakly informative and is equivalent to an inverse-gamma prior with shape and scale equal to 0.001 (Gelman 2006;Wilson et al. 2010).

Parameter sampling
Using the MCMCglmm package, each model chain was run for 50,000 iterations. As the chain begins to run the early samples may show a strong dependence on the starting parametrisation. As such, we allowed 5000 iterations to pass before the samples were stored (burn-in period) and estimates were retained every 10 iterations (thinning interval) following burn-in. We then evaluated convergence through (a) visual check of parameter time series representations, i.e. trace plots, (b) calculation of the lag k autocorrelation statistic to check lag progress and independence of posterior distribution samples, and (c) calculation of the Gelman-Rubin diagnostic statistic (comparison of four chains). All point estimates of potential scale reduction factor were < 1.1, indicating good convergence (Gelman and Rubin 1992). All models achieved convergence. Posterior distributions for the predictor variables were obtained. As model output, we report posterior means, 95% credible intervals and Bayesian p-values (p MCMC ) indicating the significance of variables with a threshold of *p MCMC < 0.05, **p MCMC < 0.01 ***p MCMC < 0.001.

Results
Across the 33 sites, 281,425 bat passes were analysed belonging to 20 different species/sonotypes from five different families (Table S1 in Online Supplementary Material). Given that there were a number of species/ sonotypes with low call numbers we present here the results for Hill numbers of order q = 2, interpreted as the effective number of dominant or very abundant species in the assemblage. Results for q = 0 and q = 1 are included in Online Supplementary Material for comparison and completeness.

Fragment size and IEM gradients
The secondary forest matrix around the smallest fragments (1 ha) harboured significantly lower TD and PD compared to CF interiors (Fig. 2). The same was true for PD in the matrix of the 10 ha fragments, with a similar pattern of erosion exhibited for FD in the 1 ha and 10 ha fragments and TD in the 10 ha fragments, albeit not significantly so. The negative effect of the matrix gradually became less pronounced with increasing fragment size, with minimal differences evident between the 100 ha matrix sites compared with CF interiors. In contrast, the edges of the 1 ha and 10 ha fragments contained significantly higher TD, FD and PD compared with CF interiors (Fig. 2). This positive edge effect extended to even the largest fragments, with the 100 ha edge sites harbouring significantly greater PD than CF interiors (see Table S5 in Online Supplementary Material for q = 2 (inverse Simpson) modelling results). Modelling results for q = 0 (Table S6 in Online Supplementary Material) and q = 1 (Table S7 in Online Supplementary Material) showed a larger number of significant effects. In particular FD was significantly reduced along with TD and PD in the matrix of the 1 ha fragments, contributing to the general pattern of significant erosion of all three diversity facets in the 1 ha fragments (Fig.  S3 in Online Supplementary Material).

Influence of local and landscape-scale predictors
Overall, the assemblage response to local-and landscape-scale variables was subtle, with scalesensitive associations for FD, whereas no significant relationships were found for TD or PD (Fig. 3). Vegetation structure does not appear to be a particularly relevant predictor of any diversity dimension. Vegetation structure had a negative relationship with FD across all diversity metrics (see Table S6; Table S7; Table S8 in Online Supplementary Material for modelling output for q = 0; q = 1; q = 2), but this manifested as a significant result only for the 2.5 and 3 km scale for q = 1 (Fig. S5 in Online Supplementary Material). FD responded significantly and positively to forest cover (compositional metric), and significantly and negatively to edge density (configurational metric) at the 1 km scale. These responses weakened with increasing buffer size. Finally, at the 2.5 km scale, patch density had a significant negative effect on FD (Fig. 3).

Discussion
Our study indicates that despite the advanced-stage secondary forest around the BDFFP fragments (~ 30 years), the diversity recovery of the aerial insectivorous bat assemblage in forest fragments and matrix habitats is still incomplete. Whilst our results indicated variation in FD, we found no clear evidence of the local-or landscape-scale variables analysed to accurately predict the diversity responses in our study area.

Responses to fragment size and IEM gradients
Despite the advanced stage of maturation of the secondary regrowth, its aerial insectivorous bat assemblage shows evidence of degradation compared to levels observed in CF interiors. These findings are consistent with similar phyllostomid bat studies at the BDFFP (Farneda et al. 2015Rocha et al. 2017b;Aninta et al. 2019). The loss of TD and FD in the secondary forest is likely to reduce ecosystem services provided by aerial insectivorous bats such as arthropod suppression (Kunz et al. 2011;Puig-Montserrat et al. 2015). Most importantly, our results support that specialist forest-dwelling bats might be strongly affected by deforestation even after ~ 30 years of forest recovery. The fact that PD is also lower in the secondary forest reflects depletion of evolutionary richness due to the loss of the overall genetic diversity of the assemblage (Struebig et al. 2011;Rivera-Ortíz et al. 2015;Edwards et al. 2021). The close relationship between the decrease in PD and the low structural complexity of secondary forests compared to CF suggests that habitat fragmentation reduces total evolutionary history by eliminating distantly related species in less complex habitats. Forest fragment interiors were either not significantly different or slightly more diverse than CF, as was the case for the 100 ha fragments. This aligns with our prediction that diversity levels in CF and the larger fragments (100 ha) would be homogeneous, contrasting previous research on phyllostomids at the BDFFP, which found significantly lower diversity in the fragments than CF (Rocha et al. 2017a;Farneda et al. 2018;Aninta et al. 2019;. As fragment size decreased, we found increasing erosion of all three facets of diversity in the intervening secondary matrix. Larger fragments (10 and 100 ha) and adjoining secondary forests tended to harbour higher TD, FD and PD than the smaller fragments (1 ha). The fact that TD, FD and PD in secondary forests next to the largest fragments (10 and 100 ha) were higher than in those adjoining smaller fragments (1 ha) suggests that, in this landscape, the largest fragments act as important reservoirs of aerial insectivorous bat diversity for the nearby regenerating areas.
In line with our prediction, all three diversity metrics peaked at the primary-secondary forest interface. Fragment edges were generally more diverse in TD, FD and PD than the CF and fragment interiors. Rodríguez-San Pedro et al. (2019) reported a similar positive response of aerial insectivorous bats to forest edges. The tall secondary forest at the BDFFP provides a low-contrast matrix that could facilitate connectivity, buffering the isolation effects of the smallest fragments and function as commuting corridors (Van Houtan et al. 2007;Jantzen and Fenton 2013;Kalcounis-Rueppell et al. 2013). Coleoptera, Hymenoptera, Lepidoptera and Araneae are among the most species-rich canopy invertebrates in tropical forests (Basset 2001). Coleoptera and Hymenoptera are more abundant in secondary forest compared to the other habitats at the BDFFP (De Aquino et al. 2021). This increased abundance of prey items could also provide increased foraging opportunities simultaneously, for edge-, open-and interior-specialists in and around the edge habitat (Ingala et al. 2021).

Influence of local-and landscape-scale predictors
We anticipated subtle responses across the three diversity facets, but we only found statistically significant responses for FD. Uncovering significant responses at the community level is often difficult due to diversity metrics amalgamating species-specific responses that may cancel each other out (Klingbeil and Willig 2009;López-Bosch et al. 2021).
Our results concur with some other studies on the effects of fragmentation on aerial insectivorous bats, in which community responses were muted. For instance, Estrada-Villegas et al. (2010) found that sonotype abundance levels were indistinguishable in a land-bridge island system which comprised interior and edge mainland sites and island sites (near and far/large and small) regardless of the sonotype group analysed and the level of island isolation or size. Falcão et al. (2021) compared activity levels and sonotype composition across two landscapes (forested vs. deforested/pasture-dominated) in the Brazilian Atlantic Forest and found that aerial insectivorous bat community richness was not related to any landscape descriptors.
We predicted that vegetation structure would have an effect at the smallest spatial scale, possibly turning neutral at large scales. However, we found no such effect on any of the diversity facets at any scale. Our findings contradict that of Blakey et al. (2017) who utilised LiDAR technology to comprehensively scan and characterise vegetation structure and found evidence that understorey forest structure was related to overall taxonomic and functional bat diversity at the community level. In a study in Panama, density of obstacles or vegetation clutter, which restrict the flight manoeuvrability of aerial insectivorous bats, were the main factor explaining both species richness and total abundance (Estrada-Villegas et al. 2012). Aerial insectivorous bats are known to utilise vertical space in a number of ways; foraging below, at and above canopy level, with much less activity within the forest interior (Marques et al. 2015). Perhaps vegetation structure is a limiting factor, albeit the effect may not be of the same magnitude as for phyllostomids at the BDFFP, where vegetation structure was a relevant predictor of total abundance (Rocha et al. 2017a).
We only observed responses to landscape composition and configuration at the assemblage level from a FD perspective. We expected that the amount of forest cover would be an important predictor of all three diversity facets at the landscape scale, but it was only important at the 1 km scale. Some recent studies analysed the influence of landscape composition and forest cover on aerial insectivorous bat assemblages in agricultural systems (Azofeifa et al. 2019;Rodríguez-San Pedro et al. 2019;Put et al. 2019), illustrating the nuanced responses to forest cover embedded in "hard" matrix types (i.e., rice fields, vineyards, grazing lands). For instance, Azofeifa et al. (2019) found forest cover to have no effect on overall species richness or composition in a study comparing two rice field sites in Venezuela. However, forest cover had a differential effect on functional groups, with background-cluttered space species responding more sensitively to changes in forest cover than uncluttered/ open space specialists. A study of vineyards in central Chile found that preservation of native vegetation cover positively affected bat diversity and species richness and that areas of landscape with more irregular and smaller patches and higher edge density facilitated a more diverse assemblage (Rodríguez-San Pedro et al. 2019).
We also found that FD responded to landscape configuration. At the 1 km scale, edge density had a negative effect, indicating that an increase in the total perimeter of edge habitat resulted in a functionally less diverse assemblage. A study quantifying edge effects across the interface of primary and secondary forest at the BDFFP suggested that the consequences of edge effects on some aerial insectivorous bat species could potentially extend for more than 2 km (Yoh et al. 2022). In contrast, at the 2.5 km scale, patch density had a positive effect, suggesting that the spatial configuration of patches at this landscape scale facilitates a more diverse assemblage. These two results for edge density and patch density agree with Chambers et al. (2016), who investigated scale dependence of habitat associations and scaling patterns of landscape metrics about bat occurrence in forests of southwestern Nicaragua and found that edge density and patch density may be as important as compositional metrics in predicting bat capture rates across multiple scales.
We have limited evidence to suggest a specific scale of effect for this community. It is reasonable to theorise that to gain access to both foraging and roost sites (landscape complementation), the species within the BDFFP community are responding, not at a single scale but across multiple scales, to both local and landscape features to secure the necessary resources (Ethier and Fahrig 2011). We suggest that further species-specific analysis might show the complexity of responses from species within an assemblage (Rodríguez-San Pedro et al. 2019). We are currently exploring species and guild level responses, which might uncover interesting and possibly conflicting patterns that may be masking an effect in this study (Gomes et al. 2020). Appel et al. (2021) investigated the interaction between habitat and moonlight at the BDFFP and analysed aerial insectivores on a speciesspecific basis. They found that Saccopteryx bilineata, Saccopteryx leptura, Centronycteris maximiliani, Cormura brevirostris, Eptesicus brasiliensis and Furipterus horrens had lower activity in secondary forest and Pteronotus alitonus and Furipterus horrens showed lower activity in fragments compared to continuous forest. It is understandable to see how these contrasting species-specific preferences can make community-level interpretation challenging.

Conclusions
Our study shows that at the BDFFP aerial insectivorous bat diversity in secondary regrowth is still not comparable with that of undisturbed forest even after ~ 30 years of recovery. The low-contrast matrix at the BDFFP does, however, appear to create opportunities for aerial-hawking bats to take advantage of forest edge habitat created as a result of fragmentation. For aerial insectivorous bat diversity, specifically, this may help to buffer some of the negative isolation effects of the smaller remnants. It also reinforces the importance of including > 10 ha forest patches in land management and conservation strategies to maximize bat diversity in human-modified landscapes. Whilst the effects of fragmentation manifest as different responses in aerial insectivores and phyllostomids, our findings reinforce the irreplaceable value of oldgrowth forest in tropical landscapes for both groups. We found significant variation in functional diversity which might be overlooked with a traditional taxonomic focus, and we observed responses to both local and landscape-scale variables at the 1 and 2.5 km scale. We therefore recommend that future studies of tropical bats follow a multidimensional biodiversity approach integrated with a multiscale analysis when assessing responses to fragmentation in human-modified landscapes.
Author contributions All authors contributed to the study conception and design. AL-B and RR collected data at the BDFFP. Data analysis was performed by AL-B and SR. The first draft of the manuscript was written by SR and all authors commented on previous versions of the manuscript and gave final approval for publication.
Funding This work was supported by the Portuguese Foundation for Science and Technology under grants PTDC/ BIABIC/111184/2009 (CM), SFRH/BD/80488/2011 and 2020.01129.CEECIND/CP1601/CT0004 (RR), PD/ BD/52597/2014 (ALB) and by the Foundation for Research Support of the State of Amazonas FAPEAM 062.01173/2015 (PEDB). Additional funding was provided by a Bat Conservation International student research fellowship to ALB and RR. This research was conducted under ICMBio (Instituto Chico Mendes de Conservação da Biodiversidade) permit (26877-2 and 26877-3) and constitutes publication 851 in the BDFFP Technical Series.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.