Composition, richness and nestedness of gallery forest bird assemblages in an Amazonian savanna landscape: lessons for conservation

Study landscape

We carried out our study in a 123,000-ha savanna landscape (0°12′N, 51°6′W) located between the municipalities of Porto Grande and Macapá, Amapá, Brazil (Fig. 1). This savanna landscape is within the Savannas of Amapá, the third largest block of Amazonian savanna (Carvalho & Mustin, 2017). Currently, this region is considered one of the last agricultural frontiers in Brazil, having been changed mainly due to the cultivation of soybeans, corn, and eucalyptus (Mustin et al., 2017; Hilário et al., 2017). The climate is hot (average temperature of 27 °C) and humid (average relative humidity of 81%) (IEPA, 2008). Average annual precipitation from 1961 to 1990 was ~2,700 mm with a well-marked dry season from August to November, when total monthly rainfall is below 50 mm (Silva et al., 1997).

The study region has three main watersheds: Matapi, Curiaú, and Pedreira (Fig. 1). The rivers and streams of the Matapi watershed connect the study region with the upland tropical forests on Precambrian crystalline terrains located in the west. In contrast, the Curuá and Pedreira watersheds’ streams and rivers connect the study region with the seasonally flooded forests of the Amazon Holocene floodplains located in the east (IEPA, 2008).

The landscape matrix is relatively well conserved. It is dominated by upland savannas and flooded grasslands intersected by gallery forest corridors, with human impacts (mostly small-scale agriculture areas and trails) limited to areas that were once occupied by upland savannas. In general, upland savannas have a ground layer dominated by grass species of the genera Rhynchospora, Axonopus, Paspalum, Polygala, Bulbostylis, and Miconia, and a woody layer (3–10 m tall) that include large shrubs and trees such as Byrsonina crassifolia, Salvertia convallariodora, Ouratea hexasperma, Curatella americana, Himatanthus articulatus, Pallicourea rigida, and Hancornia speciosa (Sanaiotti, Bridgewater & Ratter, 1997). Flooded grasslands are at the bottom of some narrow valleys, where soils are shallow and permanently inundated. These grasslands can sometimes include narrow stands of Mauritia and Mauritiella palms (Silva et al., 1997). Gallery forests are narrow (80–500 m) and found only on either cambisols or hydromorphic soils rich in organic matter along the wide valleys of rivers and streams (Silva et al., 1997; Costa-Neto, Miranda & Rocha, 2017). They are evergreen with a well-defined canopy composed of 20–30 m tall trees and a humid understory with many ferns, epiphytes, and palms. The most common plant species were Mauritia flexuosa, Euterpe oleracea, Mauritiella aculeata, Desmoncus sp., Annona paludosa, Coccoloba sp., Ficus sp., Symphonia globulifera, Virola sp., Lecythis sp., and Hymenaea parvifolia (Costa-Neto, Miranda & Rocha, 2017).

Site selection

Our study sites were 26 gallery forests associated with low-order streams distributed across the study landscape (Fig. 1). We used two criteria to choose these sites. First, we placed 8–9 sites in each of the three major watersheds in the study region to capture the regional environmental heterogeneity (Fig. 1). Second, we selected the gallery forests that were at least 1.5 km apart within each watershed to enhance spatial independence. In each site, we placed a 500 m long transect to measure habitat structure and count birds, including their abundance and richness. Perpendicularly to these forest transects, we set 400 m linear transects to measure how far species recorded in the gallery forests moved into the savanna. Each transect was marked at 50 m intervals (henceforth, “distance classes”) with colored ribbons nailed to tree trunks.

Habitat attributes

We quantified habitat attributes of gallery forests by taking measurements of both landscape and vegetation structures. To analyze the landscape structure around each site, we first mapped the landscape’s land cover/use at a spatial resolution of 10 m using Sentinel-2 images (2015–2017) and ArcGis 10.3. Using ArcGIS 10.3, we measured the gallery forest width at the transect centroid as well as estimated the percentage of forests, savannas and anthropogenic ecosystems within 500 m radius around the transect centroid. We used 500 m radius to avoid overlap with the area of influence of other transects.

To measure vegetation structure, we set three plots measuring 75 m × 2 m (150 m²) located in the beginning, middle, and end of each transect. In these plots, we took three measurements: (1) canopy cover, (2) understory foliage density, and (3) tree height.

Canopy cover was measured using a leaf coverage index ranging from 0% to 100%. This index is estimated by analyzing the canopy photographs using the program “Gap Light Analysis Mobile App” (GLAMA) (Tichý, 2016). We took four photographs at points located in each plot’s corners, with the camera vertically positioned at the height of 1.6 m above ground. The average leaf coverage of the 12 photographs (four photographs × three plots) was used to measure the site’s canopy cover.

Understory foliage density was measured as the number of foliage contacts between 0.5 and 2 m above ground with an aluminum pole 2 m tall and 5 mm in diameter. The pole was walked in front of the body in the middle of the plot (2 m wide), and all the touches of the three plots were summed and used as a measure of the site’s understory foliage density (Mantovani & Martins, 1990).

Tree height was measured using a 4.5 m graduated ruler and a Bushnell Hypsometer (Scout 1,000 Arc). The laser was focused on the highest branches or leaves, and the hypsometer data was recorded. Then, the formula

$$\text{sin}(\text{Angle}\times \pi \text{/180})\text{*Distance}\text{from}\text{the}\text{object+1}.\text{59}(\text{eye}\text{height})$$was used to obtain the height of each tree. We used the site’s average tree height to characterize each site.

Bird sampling

In the forest transects, birds were counted using 10-min, unlimited radius point counts, one of the most commonly used methods for sampling birds in tropical regions (Bibby et al., 2000; Vielliard et al., 2010). With this method, the observer records all birds seen or heard within a pre-established period. Species flying above the canopy or flying through the sample area were not recorded. A total of 78 points (three points 200-m apart within each transect) was sampled, which was replicated four times. In the savanna transects, birds were counted using the fixed distance transect method, which consists of counting all birds detected visually or aurally within 50 m on each transect side (Bibby et al., 2000).

Because we counted birds twice during the rainy season (April and June 2018) and twice during the dry season (July and September 2018), we are confident that we covered all critical periods of the region’s annual bird cycle. Birds were counted between 06:00 and 10:30, the peak period of bird activity, to maximize detection. The sequence by which transects were sampled was determined randomly one day before the sampling. J.P. surveyed all forest transects while J.C.S sampled all savanna transects. Both surveyors have extensive experience studying birds in the study region.

In each survey, the observer recorded the point location, the start and end time of the count, the recorded species’ identity (observed or heard), and the number of individuals. In the savanna transects, forest’s distance interval was also recorded. In the forest point count, the surveyor tried to avoid double counting vocal species by noting their position in relation to the points. Species flying over the transect were noted but not counted, and therefore were not included in the analysis. Olympus Binoculars (7 × 32) were used in all counts. When necessary to identify a species, vocalizations were recorded with a Zoom H1n recorder and directional microphone Yoga-Ht81.

Species functional traits

We collected data for each species’ five functional traits included in the analysis: body mass, dispersal ability, diet, foraging stratum, and habitat specificity. These five traits were chosen because previous studies have indicated that they are related to species abundance or range size across different scales, from local to global (e.g., Pigot et al., 2018; Valente & Betts, 2019; Spake et al., 2020). Body mass for all species was gathered from the literature (Wilman et al., 2014). We used the Kipp’s Index of wing morphology (KI) as a proxy for dispersal ability because this information is available for all bird species in a standard format (Sheard et al., 2020). Based on field observations and literature (Wilman et al., 2014), we classified species into four dietary groups (nectarivores, herbivores that eat mostly fruits and/or seeds, insectivores that eat primarily insects and other invertebrates, and omnivores that combine herbivore and insectivore diets) and five foraging strata (ground, understory, mid-level, canopy, and edges). As an indicator of a species’ habitat specificity, we used the weighted-average distance in which a forest species was recorded in the different category distances from the forests along the savanna transects.

Species selection

We gathered 3,770 records from 143 species in all our study sites (Table S1). However, we restricted our analyses to 99 species (representing 2,411 records) associated with forest habitats within savanna landscapes.Thus, we excluded all species considered as forest independent by Silva (1995). Because the methods that we used do not provide a reliable estimate of abundance and occupancy for all groups of birds (Bibby et al., 2000), we also excluded from our analyses all species of Psittacidae (parrots and macaws), Ramphastidae (toucans and toucanets), obligate waterbirds, raptors, nocturnal species, and aerial insectivores.

Species-level statistical analyses

We estimated abundance and occupancy for each species. Species abundance was estimated at site and at landscape level by dividing the number of points which the species was detected by the total number of points sampled (Hutto, Pletschet & Hendricks, 1986). The occupancy of a species was estimated by dividing the number of sites where the species was present by the total number of sampled sites. We used Pearson’s correlation test to assess the hypotheses that species abundance and species occupancy are associated.

To evaluate if functional traits explain the abundance and occupancy of gallery forest bird species within savanna landscapes, we used Phylogenetic Generalized Least Squares (PGLS) models to avoid problems associated with the statistical nonindependence of related species (Martins & Hansen, 1997). Foraging stratum, a categorical variable, was added to the models using a dummy coding (Mundry, 2014). Phylogenetic distances among species were estimated based on an updated version (available in http://vertlife.org/phylosubsets) of the Jetz et al’s supertree (Jetz et al., 2012) supertree based on the Hackett et al. (2008) bone. Before proceeding with PGLS, we first examined the variance inflation factors (VIF) to ensure that the predictor variables were independent. All variables presented VIF < 3 and were used in the model selection (Dormann et al., 2013). Models were generated and ranked considering Akaike’s information criterion corrected for a small sample size (AICc, Burnham & Anderson, 1998). Models with delta values (Δ_i) < 2, and high values of Akaike weights (w_i) (i.e., closest to 1), were considered to be those with the most robust support. We computed the best set of models, based on AICc (corrected Akaike Information Criterion), using the “MuMIn” package (Barton, 2020) and the model-averaging procedure. To average models, we computed mean values of estimates assuming (full averages) and not assuming (conditional average) zero values for predictors in models where they did not occur. PGLS model generation and selection were carried out in R using the PGLS function in the package “caper” (Orme et al., 2018). We have followed the recommendation of Revell (2010) and estimated the phylogenetic signal simultaneously using Pagel’s λ (Pagel, 1999) with the regression model.

Assemblage-level statistical analysis

We measured the three types of alpha diversity: taxonomic, functional, and phylogenetic. To describe alpha diversities, we used the framework described by Chiu & Chao (2014), which is based on Hill numbers. Hill numbers are defined by parameter q, which considers the relative abundance of species in determining the estimation of diversity, which facilitates the comparison of data (Hill, 1973; Chiu & Chao, 2014; Roswell, Dushoff & Winfree, 2021). In our case, we only used q values that represent taxonomic, functional, and phylogenetic richness (q = 0), where the abundance of species is ignored (Hill, 1973; Chiu & Chao, 2014). All Hill numbers were estimated with the R package “hillR” (Li, 2018). For functional richness, the Hill numbers incorporate an array of functional distances constructed from the functional traits of the species (see Chiu & Chao, 2014). For phylogenetic richness, the Hill numbers incorporate a phylogenetic tree (Li, 2018).

We used hierarchical partitioning (Chevan & Sutherland, 1991; Cisneros, Fagan & Willig, 2015) to identify the habitat attributes that best accounted for variation in each of the three dimensions of biodiversity. Statistical significance of the independent contribution of each explanatory variable was determined using a randomization approach with 1,000 iterations and an alfa-level of 0.05 (MacNaly & Walsh, 2020). Hierarchical partitioning and associated randomization tests were executed using the R package ‘hier.part’ (MacNaly & Walsh, 2020).

To test the hypothesis that less diverse bird assemblages are a nested subset of more diverse assemblages, we carried out taxonomic, phylogenetic, and functional nestedness analyses. Presence-absence matrices were first constructed where species were in the columns and sites were in the rows. Taxonomic nestedness was then estimated using the NODF index (Nestedness Metric Based on Overlap and Decreasing Fill). We chose NODF because it has more robust statistical properties than other indices and quantifies the degree to which each site is nested in each of the other sites (Almeida-Neto et al., 2008). We evaluated the significance of the taxonomic nestedness using the fixed–fixed null model (999 permutations) based on the “quasiswap” algorithm (Miklós & Podani, 2004). Both NODF estimation and the significance test were conducted using the R package ‘vegan’ (Oksanen et al., 2019).

To estimate functional (traitNODF) and phylogenetic (phyloNODF) nestedness, we used an extension of the NODF index called treeNODF index (Melo, Cianciaruso & Almeida-Neto, 2014), the same phylogeny used for the Phylogenetic Generalized Least Squares (PGLS) models, and a functional dendrogram created by using five functional traits (body mass, wing morphology, dispersal potential, diet, and foraging stratum). The functional dendrogram represents species similarity for the five functional traits and was generated from the function gawdis and UPGMA clustering algorithm. We used the gawdis function because there are problems in combining quantitative and categorical traits into multi-trait dissimilarities using Gower distance (Pavoine et al., 2009). Function gawdis balances the different traits when computing multi-trait dissimilarities, finding weights that minimize the differences in the correlation between the dissimilarity of each trait and the multi-trait (Bello et al., 2020). In general, the treeNODF index assesses the proportion of functional/phylogenetic diversity present in functionally/phylogenetically impoverished assemblages that are present in functionally/phylogenetically rich assemblages (Melo, Cianciaruso & Almeida-Neto, 2014). In addition, we partitioned the traitNODF and phyloNODF into their two components: S.fraction and topoNODF. The S.fraction represents the degree to which assemblages are or are not nested due to having assemblages composed of the same or different species. In contrast, topoNODF represents the degree to which assemblages are nested or not within the functional dendrogram or phylogenetic tree (Melo, Cianciaruso & Almeida-Neto, 2014). The treeNODF index was estimated using the R package ‘CommEcol’ (Melo, 2019). The significance of the observed traitNODF and phyloNODF and their component values (S.fraction and topoNODE) were determined using a permutation null model (999 permutations).

Species abundance and occupancy

Species abundance and species occupancy are positively correlated (Pearson’s correlation coefficient, r = 0.91, df = 97, p < 0.001). Among the 99 species included in our analyses (Table S1, Fig. S1), most of them have a low abundance index (range = 0.95–125.9, median = 12.5) and low occupancy index (range = 0.038–1, median = 0.26). Four species had the highest abundance indices (Fig. S1). Three of them were recorded in all 26 sites: a hummingbird (Phaethornis ruber), a small insectivore flycatcher (Lophotriccus galeatus), and an omnivorous thrush (Turdus leucomelas). However, the most abundant species (a small insectivore flycatcher, Tolmomyias flaviventris) was recorded in 25 sites. Among the rarest species, 26 species were recorded in less than two sites (Fig. S1). All of them had low abundance indices (range = 0.96–2.88, median = 1.92). Most of the rare species are interior forest species, including, for instance, a large tinamou (Tinamus major), six woodcreepers (Dendrexetastes rufigula, Dendrocolaptes certhia, Dendrocolaptes picumnus, Lepidocolaptes albolineatus, Nasica longirostris and Xiphorhynchus obsoletus) and two antbirds (Myrmoderus ferrugineus, Mymophylax atrothorax).

The best models (i.e., the ones with the lowest AICc) predicting species abundance (Table S2) and species occupancy (Table S3) from functional traits included four out of five functional traits examined: habitat specificity, Kipp’s Index, foraging stratum groups, and body mass. However, habitat specificity is the only one of these traits that has a positive and significant correlation with both abundance and occupancy (Table 1, Fig. 2).

Species assemblages

We detected different patterns of relationships between indicators of landscape and vegetation structure with diversity estimates (Fig. 3; Table S4). The proportion of anthropogenic area around the sites was negatively correlated with taxonomic, functional and phylogenetic diversity. In contrast, forest and savanna areas are associated positively with functional and phylogenetic diversity. Understory foliage density is negatively associated with phylogenetic diversity.

The nestedness analysis indicated that the less diverse assemblages are nested subsets of the most diverse assemblages (Fig. 4). This pattern holds when analyzing taxonomic (NODF = 53.5, p < 0.01; Figure 4), phylogenetic (phyloNODF = 64.6, p < 0.01) and functional (treeNODF = 67.1, p < 0.01) nestedness. In addition, we found that phylogenetic and functional nestedness is driven mostly by changes in taxonomic species composition (S.fraction = 52.9 and 54.6, respectively) rather than by functional or phylogenetic tree topology (topoNODF = 14.1 and 10.1, respectively).