Erosion of phylogenetic diversity in Neotropical bat assemblages: findings from a whole-ecosystem fragmentation experiment

The traditional focus on taxonomic diversity metrics for investigating species responses to habitat loss and fragmentation has limited our understanding of how biodiversity is impacted by habitat modification. This is particularly true for taxonomic groups such as bats which exhibit species-specific responses. Here, we investigate phylogenetic alpha and beta diversity of Neotropical bat assemblages across two environmental gradients, one in habitat quality and one in habitat amount. We surveyed bats in 39 sites located across a whole-ecosystem fragmentation experiment in the Brazilian Amazon, representing a gradient of habitat quality (interior-edge-matrix, hereafter IEM) in both continuous forest and forest fragments of different sizes (1, 10, and 100 ha; forest size gradient). For each habitat category, we quantified alpha and beta phylogenetic diversity, then used linear mixed-effects models and cluster analysis to explore how forest area and IEM gradient affect phylogenetic diversity. We found that the secondary forest matrix harboured significantly lower total evolutionary history compared to the fragment interiors, especially the matrix near the 1 ha fragments, containing bat assemblages with more closely related species. Forest fragments ≥ 10 ha had levels of phylogenetic richness similar to continuous forest, suggesting that large fragments retain considerable levels of evolutionary history. The edge and matrix adjacent to large fragments tend to have closely related lineages nonetheless, suggesting phylogenetic homogenization in these IEM gradient categories. Thus, despite the high mobility of bats, fragmentation still induces considerable levels of erosion of phylogenetic diversity, suggesting that the full amount of evolutionary history might not be able to persist in present-day human-modified landscapes.

Bat Sampling. Bats were sampled with ground-level mist nets in eight primary forest 133 fragments -three of 1 ha, three of 10 ha and two of 100 ha-and in nine control sites in three 134 areas of continuous forest. Fragments and controls were sampled in the interior, at the edges, 135 and in the secondary forest matrix, resulting in a total of 39 sites (for continuous forest, only 136 three edge and matrix sites were sampled; cf. Rocha et al. 2017 To assess the effect of phylogenetic information per se on PD, we eliminated any potential 166 effect of species richness by calculating SESPD (standardized effect size of PD) under a 167 'richness' null model (Swenson 2014). Similarly, we quantified phylogenetic structure per se 168 by calculating SESMPD (standardized effect size of MPD) using the tip-shuffling null model 169 (Swenson 2014). However, as the local phylogeny consisted of closely related phyllostomids 170 with quite a balanced topology, MPD could underestimate phylogenetic clustering in the 171 terminal part of the phylogeny. We therefore also calculated the mean nearest taxon distance 172 (MNTD) and computed SESMNTD using the tip-shuffling null model as it is more sensitive 173 than MPD for detecting phylogenetic clustering in the terminal part of the local phylogeny 174

Phylogenetic alpha and beta diversity across the IEM and forest size gradient 188
Alpha Diversity. We modelled phylogenetic diversity in relation to the IEM and forest size 189 gradient to understand how the two affect the evolutionary dimension of biodiversity. While  Table S1). We therefore dropped the block 195 effect, leaving us with a simple linear model to assess significance of the fixed effects. The 196 linear models were fitted using the R package 'stats' (R Core Team 2017), with IEM gradient 197 (interior, edge, or matrix) and forest size (CF, 100 ha, 10 ha, or 1 ha) as explanatory variables, Results 216

Comparison within assemblages: phylogenetic diversity is lowest in the smallest fragments 217
Rarefied PD showed a decreasing trend from continuous forest and fragment interiors to 218 forest edges and matrix, with the matrix sites adjacent to the 1 ha fragments harbouring the 219 lowest PD overall ( Fig. 1). In agreement with this, the model including IEM gradient and 220 forest size as additive predictors of PD received overwhelming support (wi = 0.99) compared 221 to the other candidate models (Table 1). The decrease in PD from interior towards edge and 222 matrix was significant (P <0.001 for both), with edges experiencing larger decreases in PD 223 compared to the matrix (Table S2). However, multiple comparison tests did not support any 224 significant differences in PD between edges and matrix, irrespective of forest size (Table S3). 225 Also, there was no significant loss in PD from larger towards smaller fragments except for 1 226 ha fragments (P <0.001, Table S2). Multiple comparison tests confirmed the significantly 227 lower phylogenetic richness in 1 ha fragments relative to larger fragments and continuous 228 forest across the entire IEM gradient (Table S3). According to the local phylogeny, most of 229 the old species lineages, i.e. those with long tree branches, do not appear in the 1 ha fragments 230 (Fig. S1). 231 When the effect of species richness on PD was accounted for (SESPD), only the IEM gradient 232 emerged as the best-supported predictor in the model (wi = 0.78, Table 1). SESPD for forest 233 fragments showed lower phylogenetic richness than expected given the number of species 234 present, particularly for fragment edges and the matrix (Fig 2A). SESPD was also significantly 235 lower in the matrix than in forest interiors, particularly for the 1 ha fragments (Fig 2A).
In line with PD, observed MPD of phyllostomid assemblages in most habitat types was lower 237 than expected under the simulated null communities, particularly for the interiors of 1 ha 238 fragments, fragment edges, and the matrix around the larger (10 and 100 ha) fragments (Fig.  239   2B). Only the assemblages in the matrix surrounding 1 ha fragments, although characterized 240 as relatively phylogenetically over-dispersed compared to the other assemblages based on 241 SESMPD, showed terminal clustering based on its significantly lower SESMNTD (Fig. 2C). 242 Patterns of MPD were best explained by a model that included only the IEM gradient as 243 explanatory variable (wi = 0.33, Table 1). There was, however, substantial model selection 244 uncertainty, with the null, size, and IEM + size models also being supported as plausible 245 (ΔAICc <2, Table 1). According to the single variable models IEM gradient and size (Table  246 S2), MPD was significantly reduced in the matrix (P = 0.03) and 1 ha fragments (P = 0.04), 247 respectively, suggesting that the matrix around 1 ha fragments contains more closely related 248 lineages as indicated by its significantly lower SESMNTD. For SESMPD, the effect of the IEM 249 gradient was even more tenuous (wi = 0.36) as the null model received the strongest support 250 (wi = 0.57, Table 1). 251

Comparison between assemblages: no clear pattern of lineage replacement 252
Dendrograms based on the total evolutionary history shared between assemblages (Pβtotal, 253 Pβrich, and Pβrepl) were substantially different from the one based on relatedness between 254 lineages within assemblages (COMDIST). UPGMA clustering based on total phylogenetic 255 beta diversity (Pβtotal) suggested that the interiors of continuous forest and forest fragments 256 harbour similar amounts of phylogenetic richness, except for the 1 ha fragments (Fig 3A), 257 further confirming the phylogenetic erosion of 1 ha fragments. The similarity in total 258 phylogenetic richness between the interiors of continuous forest and larger forest fragments 259 (10 and 100 ha) was maintained after Pβtotal was partitioned into Pβrich and Pβrepl. For Pβrich, 260 however, the interior sites of the larger fragments clustered together, unlike Pβtotal which 261 grouped the interior of continuous forest closer together with that of 100 ha fragments (Fig  262   3B). For Pβrepl, the interior sites were completely dispersed (Fig 3C), suggesting that the 263 difference in Pβtotal compared to Pβrich is caused by different lineages contained within the 264 interiors. The position of the IEM habitat categories in the Pβrepl dendrogram are unlikely due 265 to their geographic proximity as there was no significant relationship between Pβrepl and 266 geographic distance based on the Mantel test (Pearson's Mantel statistic r = 0.031, P = 0.299, 267 Table S4). 268 UPGMA clustering of COMDIST revealed that the investigated assemblages were closely 269 related and were not clustered according to either the same IEM gradient or the same forest 270 size categories (Fig. 4) and fragments. This is also unlikely due to spatial effects as there was no significant 280 relationship between COMDIST and geographic distance (Pearson's Mantel statistic r = 281 0.055, P = 0.230, Table S4). 282

Discussion 283
Using metrics of phylogenetic alpha and beta diversity, we showed that there was not only a 284 disproportionate reduction of phylogenetic richness in the smallest forest fragments (1 ha) but 285 also phylogenetic homogenization of assemblages at forest edges and in the matrix as 286 assemblages in these categories of the IEM gradient contained closely related bat lineages. 287 Whereas the two environmental gradients investigated in this study explained quite well the 288 observed variation in total phylogenetic richness (PD), they could not adequately account for 289 how closely related lineages clustered in small forest fragments (MPD). Additionally, the 290 results of the clustering analysis based on phylogenetic beta diversity metrics were counter to 291 our predictions since the resulting dendrograms did not clearly cluster assemblages from the 292 same category of habitat quality and forest size together. 293

Environmental filtering is strongest in the matrix surrounding small forest fragments 294
Across the full disturbance gradient of interior-edge-matrix, phylogenetic richness was 295 curtailed in both edge and matrix habitats, especially in the case of the 1 ha fragments. Thus Micronycteris species (Fig S1), and contain lineages that are phylogenetically clustered in the 319 terminal branches. Nevertheless, it is important to note that matrix sites adjacent to the larger 320 fragments (10 ha and 100 ha) did not differ much from their edges in terms of phylogenetic 321 richness and structure, indicating that the habitat quality gradient investigated in this study 322 represents a strong environmental filter in both large (≥10 ha) and small (1 ha) forest 323 fragments. 324 Our model selection results, however, did not provide strong support for the environmental 325 gradient investigated in this study to explain the observed phylogenetic clustering, which was 326 weak as the values of SESMPD and and SESMNTD did not deviate much from the 5 th and 95 th 327 quantiles of the null model. This is potentially due to the influence of landscape-level 328 attributes which encompass wider environmental gradients (Rocha et  incidence-based metrics, however, was limited to assess among-species variation and could 420 not account for within-species variation that is better captured by abundance-based metrics.