Emergent properties of species-habitat networks in an insular forest landscape.

Deforestation and fragmentation are pervasive drivers of biodiversity loss, but how they scale up to entire landscapes remains poorly understood. Here, we apply species-habitat networks based on species co-occurrences to test the effects of insular fragmentation on multiple taxa-medium-large mammals, small nonvolant mammals, lizards, understory birds, frogs, dung beetles, orchid bees, and trees-across 22 forest islands and three continuous forest sites within a river-damming quasi-experimental landscape in Central Amazonia. Widespread, nonrandom local species extinctions were translated into highly nested networks of low connectance and modularity. Networks' robustness considering the sequential removal of large-to-small sites was generally low; between 5% (dung beetles) and 50% (orchid bees) of species persisted when retaining only <10 ha of islands. In turn, larger sites and body size were the main attributes structuring the networks. Our results raise the prospects that insular forest fragmentation results in simplified species-habitat networks, with distinct taxa persistence to habitat loss.

INTRODUCTION

The ongoing biodiversity crisis has been primarily driven by the synergistic effects of habitat loss, fragmentation, and consequent degradation (–). Within fragmented landscapes, habitat patches can be functionally connected by interpatch species movements (). Considering habitat patches and species as connected units at the whole-landscape scale could provide innovative perceptions about ecosystem responses to human-induced disturbances, maximizing efforts on biodiversity conservation and the delivery of ecosystem services (–).

The effects of habitat fragmentation sensu lato at the entire landscape scale can be examined using species-habitat networks, which comprise a bipartite network tool recently adapted to landscape ecology by Marini et al. (). A species-habitat network represents a landscape in which each species is a node that is linked to another set of nodes matching the habitat patches where that species occurs (Fig. 1, A and B) (). The interaction links within the network are established by the occurrence of a given species at a given patch. Emerging properties of entire species-habitat networks can then be quantified using traditional metrics largely used in network studies, such as connectance, modularity, nestedness, and robustness (), further elucidating the functioning and resilience of ecosystems [see () for a recent review]. To date, this approach has yet to be tested with empirical data [but see ()].

Fig. 1.

Conceptual diagram illustrating the fragmented landscape, co-occurring species, and the resulting species-habitat networks and their properties.

The diagram represents both all taxa combined [multitaxa, (A)] and each individual taxon [single-taxon (B)] species-habitat networks based on species co-occurrences. In (A), the fragmented landscape is represented in gray, including six hypothetic habitat fragments, each of which is represented by a different color. The color code attributed to each of these fragments is maintained in (B). Network properties were examined at the landscape level (C) and node level (D), which includes both site and species traits. Landscape-level properties included network connectance, modularity, nestedness, and robustness (C), and node-level properties included normalized degree and nestedness contribution (D). At the site level, predictors included landscape-, patch-, and habitat-related metrics such as forest patch size, degree of isolation, and habitat quality. At the species level, predictors included species functional traits hypothesized to affect species persistence across the landscape such as the degree of dispersal capacity and habitat specialization.

Fragmented landscapes are typically characterized by species-area relationships (SARs) (, ), implying higher local extinction rates at smaller habitat patches [() but see ()]. Species responses to insular habitat fragmentation as induced by river damming, in which habitat patches become embedded within a uniformly hostile open-water matrix, are particularly severe (). In the aftermath of fragmentation, tropical forest specialists are especially prone to undergoing local extinctions (), while a set of habitat generalists that can cope with different forms of habitat transformation tends to persist (). Such nonrandom species reassembly is expected to disrupt habitat-species networks in the following ways (Fig. 1, A and B). First, connectance (here, the rate between actual and potential species occurrences across the focal sites within the landscape), a theoretical property related to the stability of complex systems (), is expected to be low. Second, modularity—a common descriptor of system resistance to collapse (), here interpreted as the strength of module formations comprised by subgroups of patches occupied by the same species—is expected to be low as the remaining generalist species would still connect most of the network (, , ). Third, nestedness, indicating the degree to which the network is formed by sites and species with fewer interactions (i.e., sites harboring few species and species occurring in a few sites) that are connected to species-richer sites and widely distributed species (), is expected to be high. Fourth, robustness, indicating system tolerance to secondary extinctions, corresponding to the sequential loss of species following distinct scenarios of area-driven habitat-patch removal from the landscape (, ), is expected to be low (Fig. 1C). Large tracts of tropical forests often become subdivided and isolated as small habitat patches, while large patches become increasingly rarer across the fragmented landscape (). The removal of larger forest sites is expected to exert a greater impact on network robustness (), likely inducing secondary extinctions of exclusive species (i.e., those that occur at only a single site). The loss of larger sites is likely to affect mostly species groups characterized by overall low dispersal capacity (e.g., dung beetles and frogs) and larger spatial requirements (e.g., large-medium mammals). Conversely, small forest patches proportionally harbor more species than one or a few larger patches of equal total area (). Thus, removing smaller sites from the landscape is also expected to cause secondary extinctions but to a lesser extent compared to larger sites.

Furthermore, both specific habitat features and species traits may disproportionately contribute to the emergent structure of species-habitat networks (). From either a habitat or species perspective, that is expected to be the case for species-rich habitat patches and widespread species occurring in many patches, respectively. Under the lens of network theory, those habitat patches or species would be characterized by higher degrees of centrality (, ). Likewise, habitat or species characteristics that translate into an environmental filter acting across the gradient of fragmentation, such as variation in island size or species extinction proneness (e.g., large-bodied, forest-dependent species), are expected to more strongly contribute to the overall network structure, here mediated by nestedness contribution (Fig. 1D) ().

We examined the emerging properties of species-habitat networks applied to an insular fragmented forest landscape composed of eight forest vertebrate, invertebrate, and plant taxa, including medium-large mammals, nonvolant small mammals, lizards, understory birds, frogs, dung beetles, orchid bees, and trees (Fig. 1, A and B). These taxa are critical components of tropical forest dynamics and ecosystem functioning (), ultimately providing a comprehensive synopsis of biodiversity trajectories following habitat loss and insular fragmentation (). Data on species incidence of each taxon were obtained from 22 widely distributed forest islands within the Balbina Hydroelectric Reservoir and three adjacent continuous forest sites, located in Central Brazilian Amazonia (Fig. 2). Given the expected nonrandom species extirpations across the gradient of habitat area, we hypothesized species-habitat networks at the landscape level to be characterized by low connectance, low modularity, and high nestedness (Fig. 1C). We expected network robustness to be low, particularly when larger forest sites are initially lost from the landscape, compared to either the simulated removal of smaller sites in the first place, or a random site removal irrespective of area. Statistical significance of landscape-level metrics was evaluated by applying null models that hold the number of species constant at each site but vary species composition (). We then tested whether site- and species-level attributes predict their role within the network structure. We hypothesized larger, less isolated islands and continuous forest sites harboring higher habitat quality to have higher numbers of species, as denoted by higher degree, thereby showing a stronger effect on the ordered patterns of species extirpations, as denoted by higher nestedness contributions (). In turn, species occurring at a larger number of sites, denoted by higher degree, which more strongly contribute to a nested pattern, were expected to exhibit both high dispersal capacity and generalist habitat and dietary requirements (Fig. 1D) ().

Fig. 2.

Sampling sites within the Balbina Hydroelectric Dam Landscape in Central Brazilian Amazonia.

Sampling sites include 22 forest islands (in red and highlighted by a 500-m buffer) and three continuous forest sites in the mainland (indicated by red squares). The inset map shows an aerial view of the Balbina insular landscape (photo credit: E. M. Venticinque), including the reservoir open-water matrix.

RESULTS

The multitaxa species-habitat network across the Balbina landscape included 608 species occurring at all 25 survey sites. In total, this comprised 4098 interactions, 58.3% of which were established by tree species that dominated the network (Fig. 3A). Central network positions are depicted by sites harboring more species and/or ubiquitous species occurring at a larger number of sites (i.e., higher degree) and/or by those yielding a higher contribution to network nestedness (Fig. 3). Likewise, peripheral positions are occupied by sites that contributed less to the overall structure because they were either species poor (i.e., lower degree) or did not contribute to network persistence (i.e., low nestedness contribution). Network properties at the site or species level (i.e., degree and nestedness contribution) were all highly intercorrelated (σ > 0.85).

Fig. 3.

Graphic representation of the species-habitat networks.

We considered networks both (A) combining all taxa and each taxon separately: (B) medium-large mammals, (C) small nonvolant mammals, (D) lizards, (E) understory birds, (F) frogs, (G) dung beetles, (H) orchid bees, and (I) trees. Nodes correspond to survey sites (squares, dark green colored) and species (circles, color-coded for each taxonomic group). In (A), nodes are sized according to degree (log10 x) of sites and species, and links between sites and species are color-coded according to the taxonomic group to which each species belongs. In (B) to (I), nodes are sized according to the normalized degree and color-coded for each taxonomic group. In (G), the smaller module is not scaled but amplified for visualization purposes.

Landscape-level network structure

The Balbina insular landscape was characterized by significantly low connectance considering the networks of birds, trees, and all taxa combined (Fig. 4A and table S4). Given that our null model controls for the effect of habitat area, this means that these groups were less distributed throughout the landscape than expected if their occurrence were determined solely by SARs. In contrast, connectance was not significant for the other study taxa despite the relatively wide range of observed values. For instance, dung beetles were poorly connected (C = 0.18, z = −1.85), while medium-large mammals realized almost 50% of their potential species-habitat interactions throughout the landscape (C = 0.47, z = 0.34).

Fig. 4.

Landscape-level properties characterizing the network structure of each taxonomic group and all groups combined.

Landscape-level properties include (A) connectance, (B) modularity, (C) nestedness, and (D) robustness, given the sequential loss of forest sites ordered by decreasing area. Taxa include medium-large mammals, small nonvolant mammals, lizards, understory birds, frogs, dung beetles, orchid bees, and trees and all taxa combined. Solid dots represent the observed values of each metric. Boxplots indicate the median, first and third quartiles, and minimum and maximum expected values of each metric obtained using the null model r1 (see main text for details). Red asterisks indicate statistically significant values given by −2 > z > 2 (table S4). Each taxonomic group is represented by the same icon, and the corresponding boxplot is colored as in Fig. 3. Taxa are ordered by observed values.

Medium-large mammal, small mammal, orchid bee, understory bird, and dung beetle species and all taxa combined showed modularity values significantly higher than that of the null model (Fig. 4B). Dung beetles showed the most significant modular structure (M = 0.39, z = 3.34), while medium-large mammals showed the opposite trend with low modularity (M = 0.17, z = 5.99). Conversely, trees, lizards, and frogs were more evenly distributed across the landscape, precluding the formation of distinct modules (table S4).

Assemblages of each taxon and all taxa combined, except for dung beetles and birds, were significantly and highly nested (Fig. 4C and table S4), meaning that these groups are composed of some generalist species that occur at a large number of sites, while few sites harbor a large number of species. Medium-large mammals [nestedness, based on the non-overlap and decreasing filling (NODF) = 87.5, z = 8.16] and lizards (NODF = 78.4, z = 4.66) accounted for the most nested structures (Fig. 4C).

As hypothesized, the simulated sequential removal of larger forest sites from the landscape resulted in a significant extirpation of many species of all taxa except for lizards, dung beetles, and orchid bees; yet, all significant robustness values were above 0.5 with a negative z score (Fig. 4D and table S4). This indicates that the large-to-small sequential loss of sites resulted in secondary extinctions of fewer than one species per site removed, but these networks remain less robust than expected on the basis of our null models. Robustness patterns, however, diverged markedly when sites were removed sequentially from small to large (table S4). In this case, trees and all taxa combined were more robust than expected from our null models, indicating that removing small sites would be buffered by large sites, which still contained most species, given the overall nested assemblage structure. Species declined much faster when sites were removed from the largest to the smallest (Fig. 5A), compared to the reverse sequence (Fig. 5B). For instance, if forest loss occurred at all islands larger than 10 ha, which represent nearly half of all islands in the Balbina archipelago, species persistence would be markedly reduced across the entire landscape, particularly for birds (~20% persisting species) and dung beetles (~5%) (Fig. 5A). No significant robustness values were observed when forest sites were sequentially lost at random (table S4).

Fig. 5.

Proportion (%) of species persisting across the Balbina landscape.

Two alternative scenarios are provided given the sequential loss of forest sites ordered by decreasing (A) and increasing forest area (B). The number of secondary species extinctions following habitat removal was obtained using the second.extinct function () from the bipartite R package (), which considers the total number of species recorded for each taxon. Results are indicated for each taxon: medium-large mammals, small nonvolant mammals, lizards, understory birds, frogs, dung beetles, orchid bees, and trees. Each taxon is represented by the same icon and corresponding color-coded circles, as in Fig. 3. Lines correspond to the adjusted trend using the geom_smooth function considering span = 1.5, based on the ggplot2 R package ().

Predictors at the site level

Forest area was consistently a critical predictor of site degree across all taxa (Fig. 6A and fig. S1), with medium-large mammals showing the steepest relationship (Fig. 7A). Site degree was further affected by the interaction between forest area and distance to the mainland, with positive and negative relationships for lizards and small mammals, respectively. Proportion of closed-canopy cover was positively related for lizards, birds, and dung beetles. Burn severity was also positively related for birds, as was the interaction of burn severity and forest area for medium-large mammals, beetles, and trees (Fig. 6A). Determinants of species richness for all taxa combined further reflected those observed for trees (table S5). Nonetheless, forest area was also the most important driver structuring the nested species-habitat networks, affecting all taxa, except for lizard and orchid bee assemblages (Fig. 6B). While the effect of forest area was positive for most taxonomic groups, the nestedness contribution decreased with forest area for medium-large mammals (Fig. 7B and fig. S2). Trees were additionally affected by the interaction between forest area and burn severity, and all taxa combined were affected by the proportion of closed-canopy cover (table S6).

Fig. 6.

Predictors of site-level network properties for each taxon and all taxa combined.

Network properties included (A) site normalized degree and (B) site nestedness contribution. Predictors included forest area (log10 x; Area), distance to mainland continuous forest (Dist), proportion of closed-canopy forest (%Closed-canopy), and burn severity (Burn). The interaction terms between Area and Dist and Area and Burn were retained in the models only in case this ensured a significant explanation of the response variable, represented by the asterisks. Circles are sized according to the estimate obtained from averaged models (see further details on model results in tables S5 and S6). Blue and red circles denote positive and negative estimates, respectively. Only statistically significant variables (P ≤ 0.05) in average models are shown. Each taxon is represented by the same icon as in Fig. 3.

Fig. 7.

Relationship between species-level network properties and forest area (log10 x; hectares).

Network properties included (A) species normalized degree and (B) species nestedness contribution. Lines represent the adjusted linear model between the network property and forest area. Each taxon is represented by the same icon, and corresponding circles and lines are color-coded as in Fig. 3.

Predictors at the species level

Both species degree and its contribution to nestedness were related to specific species traits (Table 1). In particular, both wide-ranging medium-large mammals and small-bodied birds and frogs occurred at fewer sites (lower degree). In contrast, large-bodied orchid bees occurred at a large number of sites (higher degree) (figs. S3 and S4A). Tree species typical of the subcanopy and canopy forest strata occupied fewer sites. Bird species that most contributed to nestedness were characterized by small body size and higher forest habitat specificity (figs. S4B and S5). Tree species most contributing to nestedness were primarily large seeded and characterized by an old-growth strategy in the shaded understory (Table 1).

Table 1.

Results of model averaging explaining species normalized degree and nestedness contribution at the species level for each taxon according to their traits.

For each taxon and each network property, we provide the estimate, SE, and P value of each functional trait. Species normalized degree of understory birds, dung beetles, and orchid bees was log-transformed. All models were fitted with a normal distribution, except for those on the nestedness contribution of both medium-large mammals and frogs, which were fitted with a log-link function. Results on the model intercept are provided only when any categorical trait was included in the corresponding model. Significant results (P ≤ 0.05) are shown in bold.

Taxon	Species traits	Normalized degree			Nestedness contribution
		Estimate	SE	P value	Estimate	SE	P value
Medium-large mammals
	Home range size (log10 x)	−1.066	0.048	0.035	−0.026	0.048	0.610
	Trophic level	−0.076	0.051	0.148	−0.045	0.044	0.333
	Body mass (log10 x)	0.020	0.054	0.716	−0.032	0.045	0.496
	Group size (log10 x)	−0.038	0.046	0.432	0.038	0.041	0.376
Small mammals
	Body mass (log10 x)	0.078	0.052	0.167	−0.200	0.322	0.560
	Trophic level	0.048	0.055	0.418	0.085	0.311	0.800
	Matrix tolerance (log10 x)	0.003	0.059	0.959	−0.514	0.285	0.092
	Locomotion habitat	−0.002	0.059	0.964	−0.235	0.308	0.477
Lizards
	Prey size	0.135	0.078	0.120	0.093	0.058	0.150
	Thermoregulation mode	−0.085	0.087	0.375	−0.068	0.062	0.321
	Body length (log10 x)	−0.002	0.089	0.984	−0.028	0.064	0.690
Understory birds
	Body mass (log10 x)	−0.080	0.036	0.028	−0.314	0.112	0.006
	Trophic level	−0.036	0.038	0.350	0.215	0.118	0.070
	Forest dependency	0.001	0.037	0.988	0.401	0.116	0.001
Frogs
	Body size (log10 x)	−0.114	0.041	0.007	−0.084	0.052	0.121
	Reproduction mode	0.059	0.042	0.173	0.055	0.051	0.299
	Habitat diversity	−0.031	0.041	0.478	−0.050	0.051	0.343
Dung beetles
	Body size (log10 x)	−0.046	0.061	0.471	0.029	0.175	0.872
	Diet: generalist	−0.006	0.127	0.964	−0.012	0.363	0.974
	Relocation: roller	0.155	0.196	0.449	−0.413	0.551	0.473
	Relocation: tunneler	0.157	0.173	0.385	−0.581	0.487	0.253
Orchid bees
	Body length	0.033	0.016	0.060	0.064	0.049	0.205
Trees
	Seed mass	0.051	0.029	0.077	0.211	0.086	0.015
	Vertical strata: emergent	0.038	0.094	0.690	−0.025	0.295	0.932
	Vertical strata: subcanopy	−0.184	0.056	0.001	−0.182	0.176	0.301
	Regeneration: old growth	0.087	0.066	0.189	0.410	0.207	0.049
	Regeneration: pioneer	0.133	0.102	0.195	0.192	0.319	0.548
	Wood density	−0.022	0.027	0.413	−0.106	0.086	0.219
	Dispersal mode: biotic	0.125	0.080	0.117	0.407	0.246	0.100
	Dispersal mode: both	0.100	0.118	0.397	0.121	0.363	0.739

DISCUSSION

Our study revealed the emergent properties of species-habitat networks considering the spatially correlated trajectories of multiple vertebrate, invertebrate, and plant taxa across a unique quasi-experimental fragmented tropical forest landscape created by a large dam. Overall, widespread, nonrandom local species extinctions were depicted by the low connectance and high nestedness of the networks, indicating weak tolerance to further habitat loss and fragmentation as shown by low robustness values when forest loss initially penalized continuous forest sites and the largest islands. Our innovative approach amounts to a bipartite network in which species and forest sites are potentially connected, given patch-scale assemblages across the gradient of fragmentation. This ensured that we could characterize a complex landscape as a unit. From a theoretical perspective, we used a unique empirical dataset, including a large number of phylogenetically independent taxa surveyed on an exceptional number of sites, to reveal the properties of species-habitat networks in a man-made archipelagic landscape. This allowed us to unveil previously unknown patterns, such as the network simplification that characterizes species-habitat networks (, ). By considering all eight vertebrate, invertebrate, and plant taxa, we further encompass several ecological functions implicit within the system, such as seed dispersal, pollination, arthropod pest control, and herbivory, providing us with a holistic assessment of ecosystem responses to human-induced disturbances.

The low network connectance in Balbina indicates that a considerable number of all possible site-species interactions have been discontinued (), illustrating that islands within human-induced insular landscapes succumb to systemic local extinctions (), as broadly observed in other fragmented landscapes (). Our site-level analyses further demonstrate that these local extinctions were consistently driven by forest area, whereas the interaction between island size and isolation, coupled with indicators of habitat quality (forest canopy continuity and burn severity), was poorly informative for each taxon. Consistent with the Island Biogeography Theory, we would expect a larger role of island distance from the mainland (), but our findings are in agreement with previous expectations for forest patches surrounded by an inhospitable aquatic matrix (, ). Notwithstanding, connectance was lower than expected by chance only for birds, trees, and all taxa combined. Birds and trees contained the largest number of species (N birds: 119; N trees: 335), whereas all taxa combined were represented by 154 species. The lack of significance in connectance for all other taxa is likely due to the choice of a conservative null model, which constrained the number of species per site (), thereby leaving fewer possibilities for randomization in small networks. From a species perspective, functional traits also contributed to species persistence across the fragmented landscape. For instance, medium-large mammals with limited spatial requirements occurred at a higher number of sites, suggesting that wide-ranging species were primarily restricted to large forest patches (). Similarly, most widespread bird and frog species were small bodied (). Conversely, large-bodied orchid bees likely benefit from prolonged flight autonomy and high dispersal capacity ensuring that they could traverse greater matrix distances in search of resources (, ), facilitating their persistence in a larger number of sites. In addition, (sub)canopy tree species were restricted to a small number of sites, likely due to stand replacement on more degraded small islands, which are prevalent in the study area (95% of all islands are <100 ha) ().

Nestedness was the most preeminent property of the Balbina species-habitat networks. Here, highly nested networks indicate that the landscape is composed by a core of species-rich sites and widespread distributed species to which species-poor sites and species occurring at few sites are associated with (). The generalist species–habitat interactions (i.e., between species-rich sites and ubiquitous species) play an important functional role in this system as they can potentially (re)connect the network, maintaining or even reestablishing lost ecological functions performed by interacting species. These widespread distributions can maintain even small, highly isolated islands connected to the entire network, thus contributing to the overall functional connectivity across the landscape. This also indicates that habitat-specialist interactions, which would form more modular structures, are rare, as occurred to some extent for five of the eight taxa. In general, modularity showed relatively low but significant values for these five groups. As an exception, the dung beetle network was more modular, given that a unique dung beetle species occurred only at a single island. Moreover, nestedness was highest for medium-large mammal networks, followed by lizards and small mammals. This supports the notion that spatial nestedness tends to be higher for species foraging at higher trophic levels (), which occur at lower densities and require larger areas of contiguous habitat ().

Site characteristics and species traits were also important in contributing to network nestedness. For instance, forest area played a major role in driving nestedness for most taxa. At the species level, large-bodied, forest-dependent bird species favored communities holding a nested subset of the links in better connected communities (). Likewise, old-growth, large-seeded tree species, which most contributed to functional nestedness, tend to interact with a smaller coterie of highly specialized dispersal agents, such as mid-sized to large-bodied frugivores. These potential dispersers are more extinction prone on small islands (), limiting forest regeneration therein, as observed in other semi-defaunated tropical forest remnants (). Furthermore, when all taxa were combined, network properties were very similar to those observed for trees, likely because of the high species richness of this group. Similarly, tree assemblages reflected the SAR applied to multiple taxa (). One interesting exception was, however, the importance of habitat quality. Our results showed a positive effect of the proportion of closed-canopy forest on network nestedness, which only became evident when all taxa were combined. This reinforces the value of high-quality habitat in maintaining tropical forest biodiversity.

Our simulated large-to-small removal of forest sites from the network showed a consistent decline in the proportion of remaining species richness for all taxa, which was particularly severe for birds and dung beetles. Both of these groups showed higher levels of modularity and nonsignificant nested structures, shedding light on the open network-stability debate (, ). Our findings point to higher robustness of nested structures following the loss of the most connected nodes within the overall network. Tropical forest birds have already been pinpointed as one of the most sensitive groups to habitat change (). At the site level, forest area had the strongest effect on the degree for both birds and dung beetles, which, unlike other groups, were affected by variables related to habitat quality (canopy openness, burn severity, and its interaction with island area). It is possible that such habitat dependency of birds and dung beetles reduces the stability of their species-habitat networks that are most sensitive to changes under those habitat conditions. Species-habitat networks in a fragmented landscape also highlighted the importance of habitat quality for avian assemblages (). On the other hand, orchid bees comprised the most robust networks. Balbina orchid bees showed a flat slope in their SAR (), being relatively unaffected by forest area, which explained only 10% of their rarefied local species richness (). This is likely related to the higher dispersal capacity of this group, especially the most abundant large-bodied species (e.g., Eulaema bombiformis and Eulaema meriana), which continues to potentially connect the overall landscape (i.e., higher species degree). Over and above differences between taxa, the overall low robustness of the species-habitat networks most evident when habitat loss was ordered from the largest to the smallest site highlights their sensitivity to further habitat loss and fragmentation (). In tropical forest landscapes, large tracts of intact forests are often converted into multiple small forest patches, whereas very small patches tend to vanish from the landscape (). Given that small islands rarely harbor exclusive species, our results forecast scenarios in which removal of relatively small islands (<10 ha), which are widespread across the Balbina reservoir, would exert only a minor impact on the proportion of species persisting in the network. However, if forest loss initially occurred in relatively large islands (100 to 1000 ha), which are few and far between, then a collapse in the species-habitat networks would occur for multiple taxa with profound implications for ecosystem functioning across the entire archipelagic landscape.

Our network approach is, however, limited in precluding inferences on the actual functional connectivity of the landscape. Connections between any two sites within the network did not necessarily correspond to “real” linkages (movement of individuals) across the landscape (). Ensuring that links connecting sites within the networks illustrate the functional connectivity of the landscape would require considering both species dispersal abilities and space. Although there have been attempts to measure species dispersal capacity [e.g., observations of mammals swimming between land masses () and orchid bee records within the aquatic matrix ()], data on dispersal are severely limited for most taxa at Balbina and other tropical landscape (). Inferences on meta-community dynamics are also limited, given by the fact that islands surveyed in this study were not contiguous to one another. Hence, we cannot rule out the possibility that unsurveyed islands could facilitate functional connectivity.

In conclusion, we show simplified species-habitat networks and a high prevalence of nonrandom local species extinctions for the different vertebrate, invertebrate, and plant taxa across a vast archipelagic landscape. However, species persistence was highly variable across taxa. Results at the node level are consistent with more traditional approaches (), reinforcing the dominant role of forest patch size in determining the persistence of species diversity in fragmented landscapes. Given the growing trends in habitat insularization in the aftermath of river damming, future hydropower development should avoid flooding large expanses of forest and creating myriad small islands. Incorporating such recommendations would minimize biodiversity loss and consequently buffer disruptions in ecosystem functioning at the scale of entire landscapes.

MATERIALS AND METHODS

Study area

This study was carried within and around the insular landscape of the Balbina Hydroelectric Reservoir in Central Brazilian Amazonia (1°48′S, 59°29′W; fig. S7). The reservoir was created in 1987 following the damming of the Uatumã River, a left-bank tributary of the Amazon River. In total, 312,900 ha of primary forest was flooded within the 443,772-ha hydroelectric reservoir, leaving 3546 islands corresponding to the former hilltops. Most islands and surrounding continuous forests consist of dense closed-canopy terra firme forest. Although part of the reservoir is legally protected by biological reserve, many islands experienced ephemeral understory fires during the 1997–1998 El Niño drought (). The local mean annual temperature is 28°C, and the mean annual rainfall is ~2376 mm (). Data on eight taxa were obtained at 22 islands. Island selection aimed to achieve a wide range of spatial configurations in terms of island size (means ± SD: 223 ± 364 ha; range = 1.52 to 1460 ha) and isolation distances to the nearest continuous forest in the mainland (4351 ± 3386 m; range = 44 to 11,972 m; table S1). As a baseline, we also surveyed three sites that were widely spaced in the surrounding mainland continuous forest (Fig. 2). Including continuous forest sites in the mainland in our analysis further allowed us to represent the extreme end of the gradient of habitat loss and fragmentation.

Dataset

Species incidence data for eight taxonomic groups were collected between 2011 and 2016, across the 25 surveyed sites within the Balbina landscape. Surveyed islands were not adjacent to each other, and their variation in size reflects the overall island size distribution in Balbina (95% of all islands are <100 ha). Sampling effort varied with the forest area of each site for medium-large mammals, small mammals, lizards, frogs, orchid bees, and trees, but not for birds and dung beetles, which followed a standardized sampling effort (see details on taxon-specific sampling effort in table S1). We provide brief descriptions on survey data acquisition below, but the related references should be consulted for full details including those on the authority granting ethics approval and permits obtained. Medium- and large-sized diurnal and nocturnal mammals (hereafter, medium-large mammals) correspond to those of >100 g and detectable by line transect censuses, indirect sign surveys, and/or camera trapping. A total of 28 medium-large mammal species were recorded including primates, carnivores, xenarthrans, ungulates, and rodents (). Small-mammal surveys were conducted along transects using live traps set at ground, understory, and subcanopy heights and 100-liter pitfall units connected by a drift fence, resulting in a total of 20 species recorded, including rodents and marsupials (). Terrestrial diurnal lizard assemblages were sampled using 100-liter pitfall traps connected by a drift fence (), and a total of 14 species were recorded. Understory birds were surveyed using mist nets that allowed us to record 119 species (). Frog assemblages were surveyed using autonomous recording units, and their calls allowed the identification of 35 species (). Dung beetles were represented by 33 species surveyed using pitfall traps baited with a mixture of pig and human dung (). Euglossine or orchid bees (24 species) were surveyed using scent trap arrays baited with cineole, methyl salicylate, methyl cinnamate, and vanillin (). Large trees (including arborescent palms) were surveyed using floristic inventories carried out within quarter-hectare vegetation plots. All live trees ≥10 cm in diameter at breast height within each plot were measured, number-tagged, and identified to species level by an expert botanist, resulting in 335 species ().

Landscape, patch, and habitat quality metrics

All metrics, except for burn severity, were derived from multispectral RapidEye imagery (5-m resolution with five-band color imagery) from 2011 and 2012, using ArcMap 10.1 (). To do so, we conducted a semisupervised classification to obtain four land cover classes (closed-canopy forest, open-canopy forest, bare ground, and water) [for further details on image acquisition and processing, see ()]. At the landscape scale, we measured the proportion of open water within a 500-m buffer threshold surrounding each focal island (%Water). We considered this buffer size to maximize spatial independence between neighboring sites. At the patch scale, we measured the following: forest area (Area), defined as the island area excluding bare ground; island shape, defined as the Perimeter:Area ratio (); and distance to continuous forest (Dist), given by the Euclidian distance from the edge of each focal island to the nearest lake margin. Regarding habitat quality, we obtained a measure of ground burn severity (Burn) and the proportion of closed-canopy forest (%Closed-canopy) within each site. Burn was scored in situ on the basis of the extent to which each forest site had been affected by previous surface fires and the number of charred trees and height of char marks on each tree (). Because several analogous patch variables could not be obtained for the three continuous forest sites, we assigned their metrics to closely approximate “real-world” values. Thus, we assumed Area to be equal to one order of magnitude (10-fold) larger than our largest island, Dist to equal zero, and %Water to be the minimum value of 0% (table S2).

Species traits

To examine the effects of species traits on species-level network properties, we considered a set of uncorrelated traits that are primarily related to the inherent species capacity to disperse and persist in the insular landscape. Traits included body size (all taxa), trophic level (all but lizards, frogs, and orchid bees), home range size and group size (medium-large mammals), degree of forest dependency (small mammals and birds), thermoregulation mode (lizards), and preferred habitat type (small mammals, lizards, and frogs). For trees, we considered wood density, seed mass, forest vertical strata, seed dispersal mode, and seedling regeneration strategy. Part of the species traits was measured in situ using individuals captured during surveys, whereas other traits were extracted from available databases and the literature (for a description and source of each species trait, see table S3).

For each taxon and all taxa combined, we compiled an a adjacency matrix in which i corresponds to a sampling site, j corresponds to a given species, and the ij element is filled with either the presence or absence of species j at site i. We first obtained a graphic representation of the species-habitat networks, using the igraph R package (). We then characterized the network structure by estimating the (i) connectance, as the realized proportion of species-habitat interactions given all potential interactions () obtained using the networklevel function; (ii) modularity, which tests whether species-habitat interactions form distinct and cohesive subgroups that strongly interact among themselves compared to nodes in other regions of the structural network configuration (), was obtained by applying the computeModules function, which uses the algorithm provided in (); (iii) nestedness, based on the non-overlap and decreasing filling (NODF) of the presence-absence of data in the observed adjacency matrix (), was also obtained using the networklevel function; and (iv) robustness, which provides a measure of resilience of the species-habitat network to disturbances mediated by both the random loss of forest sites and the sequential loss of sites (ordered by both decreasing and increasing forest patch area). Network robustness is interpreted as the proportion of all species remaining after the simulated removal of any given site from the network. We obtained the area below the extinction curve and the measure of robustness using the second.extinct () and the robustness functions (), respectively.

Given that we aimed to examine habitat fragmentation effects other than the canonical SAR (), we contrasted the observed values for each of the network metrics against random matrices generated using the “r1” null model (). This model (i) fixes row marginal totals, maintaining the number of species constant at each site i; (ii) fixes the number of columns and thus the number of j species in the matrix; and (iii) allows the column marginal totals to vary, meaning that the number of sites in which species j occurs can vary, but the number of species occupying each site i remains constant. In summary, this null model randomizes species composition at each site, so that contrasting with our observed values ensured an understanding of whether local extinctions were not a random process. We ran 100 iterations for each taxonomic group and all taxa combined, using the function oecosimu in the vegan R package (). Statistical differences between observed and random expectations were assessed using z scores.

Node-level network structure

For both site- and species-level (i.e., node-level) analyses, we calculated two centrality properties: (i) normalized degree (hereafter, degree), corresponding to, at the site level, the number of species occupying each site divided by the total number of species and, at the species level, the number of sites occupied by each species divided by the total number of sites (). Using the normalized version of degree ensures minimal disparity between intrinsic differences in the numbers of species in each taxon, and (ii) nestedness contribution, at the site level, denotes the extent to which the interaction links established by each site i affect network nestedness and, at the species level, the extent to which the interaction links established by species j affect nestedness contribution (). In particular, for each row and column, this function compares observed nestedness to an ensemble of nestedness values generated by randomizing the interactions of a given focal site or species. Nestedness contributions are the z scores from this comparison (). Calculations are based on the NODF index and use a probabilistic null model to randomize interactions (). These properties were estimated using the specieslevel and nestedness.contribution functions, respectively. All functions used to obtain landscape- and species-level network properties can be found in the bipartite R package ().

Statistical analyses

Predictors at the site level

Before statistical analyses, we checked the distribution of each network property for each and all taxonomic groups combined. We log-transform network properties whenever this ensured improved model fitting. To control for high levels of variable interdependence, we performed a Pearson correlation matrix; %Water and Shape were highly correlated with Area (r > 0.75). Given that Area is a strong predictor of species diversity for all taxa (), we further excluded %Water and Shape from subsequent analyses. Explanatory variables were previously standardized (mean = 0, SD = 1) to place coefficient estimates onto the same scale.

To examine the effects of habitat quality, patch, and landscape variables on both normalized degree and nestedness contribution, we first accounted for potential spatial autocorrelation, applying generalized least squares (GLS) models using the REML method in the nlme R package (). For each network property of each taxonomic group and all groups combined, we applied six GLS models, each of which included one of five fixed variance structures (i.e., spherical, linear, rational quadratic, Gaussian, and exponential correlations) both with and without the model spatial structure represented by the geographic coordinates of each site (latitude and longitude) (). GLS models were compared on the basis of their Akaike Information Criteria corrected for small sample size (AICc). The model with the lowest AICc was selected, except if ΔAICc ≤ 2 (ΔAICc = AICc − AICcmin in which i = ith model), in which case, we performed an analysis of variance (ANOVA) with all models with ΔAICc ≤ 2. If no differences were observed between those models, then we selected the one with the simplest structure. If the best model excluded the fixed variance structure, then we used a linear model (LM) or a generalized linear model (GLM) according to the best fitted distribution. We then constructed a candidate model set using all additive combinations of the four explanatory variables retained (Area, Dist, %Closed-canopy, and Burn), and models were ranked on the basis of their AICc, using the MuMIn R package (). As observed for some taxa in Balbina, the effect of island size (Area) on species diversity can be further boosted by their isolation distance (Dist) [e.g., orchid bees: (); dung beetles: ()] or their degree of burning severity (Burn) [trees: ()]. Hence, we additionally ran the model selection using the interaction term between Area and Dist and Area and Burn, which were retained only in case this ensured a significant explanation of the response variable. To account for model uncertainty in multimodel inference, a model-averaging approach was performed using all alternative models.

Predictors at the species-level

We related the species-level normalized degree and nestedness contribution with the species traits for each taxonomic group. We inspected data distributions and high levels of variable interdependence as performed at the site level. Lizard thermoregulation mode was highly correlated with habitat type (r > 0.70), and body length of bees was highly correlated with body width and wingspan, so we included only one of those traits (body length) in subsequent analyses. We then applied LM or GLM according to the data distribution fitted into each of those models. Model selection followed the same structure described above for predictors at the site level. All data analyses were performed in R ().