Harvesting of forest products and implications for Afrotemperate bird communities in a montane forest of the Eastern Cape,South Africa

Jessica Leaver ,Johann C.Carstensand Michael I.Cherry

Abstract

Keywords: Human disturbance, Habitat heterogeneity, Habitat modification, Beta-diversity, Elevation, Generalist species, Forest-specialist species

Introduction

Human disturbances and their impact on forest habitats are threatening biodiversity (Bradshaw et al. 2009;Newbold et al. 2014). In developing regions, harvesting of forest products represents the most widespread human disturbance in natural forests (Vermeulen 1996; Luoga et al. 2000; von Maltitz 2003; Kumar and Shahabuddin 2005; Lawes et al. 2007a). While forest management policies in many developing nations aim to balance the socioeconomic benefits of resource use with the conservation of forest biodiversity (Shackleton et al. 2002; Robertson and Lawes 2005), regulation of resource use is often limited, such that de facto open-access systems prevail(Thapa and Weber 1995; Pandit and Thapa 2004; Robertson and Lawes 2005;Sunderlin et al.2005).Several studies have investigated the ecological implications of unregulated forest resource use, revealing significant impacts on forest habitats, from population-level declines of target species (Guedje et al. 2007; Williams et al. 2013) to community-level changes in floristics and structure(Kumar and Shahabuddin 2005; Sassen and Sheil 2013).While these studies show that resource use is a major driver of habitat modification, our understanding of concomitant effects on faunal biodiversity is limited (Laiolo 2004; Shahabuddin and Kumar 2006, 2007; Gardner et al.2016;Asefa et al.2017).Birds are particularly good indicators of environmental change (Gregory and Strien 2010),as well as being essential for the function and regeneration of forest ecosystems(Pimm 1986;Sekercioglu 2006).

Understanding the impact of habitat modification on forest avifaunal communities is challenging given its dynamic nature, and thus requires the consideration of multiple factors.First, avifaunal responses depend on the nature and intensity of habitat change, specifically regarding implications for habitat heterogeneity (Stirnemann et al. 2015; Murray et al. 2017; Schulze et al.2019). For example, disturbances that maintain or enhance habitat heterogeneity may maintain or increase avifaunal abundance and diversity by providing a diversity of resources and niches (Seymour and Dean 2010; Murray et al. 2017; Schulze et al. 2019), while disturbances that reduce habitat complexity are more likely to promote biotic homogenization (Arroyo-Rodríguez et al. 2013;Morante-Filho et al. 2016). Second, responses to habitat change may differ across species based largely on their level of habitat specialisation, with habitat specialists thought to be more sensitive to disturbance than generalists (Devictor et al. 2008b; Clavel et al. 2011). The loss of disturbancesensitive specialist species in response to habitat change may be compensated for by an increase in disturbanceadapted generalist species, thereby masking communitylevel responses to disturbance(Supp and Ernest 2014).Furthermore, given that habitat specialists are generally less wide-spread than other species, they are of greater conservation concern. Third, levels of human-mediated habitat disturbance are often correlated with environmental gradients along which forest environments occur, such as elevation, given that more accessible areas, i.e. at lower elevations, are more likely to be disturbed (Montano-Centellas and Garitano-Zavala 2015). Avifaunal responses to human disturbances may thus be confounded by correlated changes in elevation, given that avian species richness and composition change along elevational gradients due to associated natural changes in habitat conditions (Terborgh 1977; Jankowski et al. 2009, 2013; Monta?o-Centellas and Garitano-Zavala 2015; Hui et al. 2018). Despite the prevalent interaction of elevation and disturbance gradients in forests, few studies have aimed to separate their respective effects on bird communities (Montano-Centellas and Garitano-Zavala 2015). Last, variation in species composition across locations(i.e.beta-diversity)is an important determinant of the number of species that can accumulate at greater scales,and thus a vital component of understanding biodiversity responses to disturbance (Arroyo-Rodríguez et al. 2013). However studies assessing human impacts on avifaunal diversity often overlook the impact that habitat disturbances may have on beta-diversity (Morante-Filho et al.2015).Moreover,few studies consider the impact that habitat disturbances may have on the mechanisms that drive variation in species composition, namely species loss and species turnover (Baselga 2010, 2012). This is important because understanding the mechanisms through which disturbances affect variation in species composition largely informs what conservation actions are necessary.

The Eastern Cape, South Africa harbours 46% of the country’s remaining natural forest cover, and falls within the Maputaland-Pondoland-Albany biodiversity hotspot(Berliner 2009). Thus, while rich in biodiversity, the region is economically one of South Africa’s poorest and least developed provinces, with high levels of unemployment and rural poverty (Statistics South Africa 2018). In this economically impoverished context, forest resources comprise a critical contribution to livelihood strategies for communities close to forests (Shackleton and Shackleton 2004; Stadler 2012). Poorer households in particular may have a high dependence on forest resources for fuelwood, building and fencing material, medicine, food,and increasingly,income earning opportunities through the commercialisation of certain forest products (Paumgarten and Shackleton 2009, 2011; Stadler 2012), particularly medicinal plants (Dold and Cocks 2002; Geldenhuys 2004;Williams et al. 2013). While the National Forest Act(1998) recognises the socio-economic importance of forest resources, and aims to manage natural forests sustainably, the Department of Environment, Forestry,and Fisheries (DEFF), which manages 70% of South Africa’s natural forests, exercises little regulation of resource use (Obiri and Lawes 2002). In the context of increasing commercialisation of medicinal plants, there is increasing concern that unregulated resource use in South Africa’s natural forests is degrading forest habitats (Hoppe-Speer et al. 2015) and compromising forest biodiversity (Castley and Kerley 1996; Krüger and Lawes 1997; Leaver et al. 2019).

In this study, we examine forest disturbance due to resource use,impacts on habitat structure,and concomitant effects on avifaunal communities in Gomo, a representative Afromontane forest occurring along an elevational gradient in the Eastern Cape, South Africa. Specifically,this study comprises of two major components: first, we compare harvest disturbances, resultant harvest-mediated habitat heterogeneity, and avifaunal species richness,composition and beta-diversity across two elevational zones within the forest. Second, we test for associations between harvest-mediated habitat heterogeneity and avifaunal community richness and beta-diversity, while controlling for elevation, to assess whether humanmediated habitat modification influence observed patterns of spatial diversity at the forest-scale; and the mechanisms driving this.

Materials and methods Study site

This study was conducted at Gomo forest (31°0′39.34′′S, 29°20′44.25′′E) within the Alfred Nzo district, in the northern, inland zone of the Eastern Cape Province,South Africa (Fig. 1a and b). The topography of this region is mountainous, with fragmented forest patches within a grassland matrix, together with smaller stands of commercial pine plantations; and scattered rural settlements (von Maltitz et al. 2003). Indigenous forests in this region are classified as Transkei mistbelt forests,and occur along a fragmented band at mid-altitudes(850-1600 m above sea level), confined to fire refugia on south, south-eastern and south-western mountain slopes(Fig. 1b). Mean annual rainfall in the region varies between 600 and 1200 mm, with rainfall occurring predominantly in the summer, between October and March.Heavy mists also occur during these months,contributing to the moist summer conditions in the region. Temperatures are mild, with mean annual temperatures ranging from 14°C to 18°C, although temperatures during the winter months can drop to 2°C, and occasional winter snowfall occurs.

Socio-economically, this district falls entirely within the former homeland of the Transkei, and remains characterised by a weak infrastructure and limited economic opportunities. Consequently, this district has some of the highest unemployment rates in the province (43.5%), with close to 40% of households having monthly incomes below the poverty line, i.e.less than R 800 (～ $56.00) per month (Statistics South Africa 2018). Subsequently, Transkei mistbelt forests represent the inland forest type under the highest resource use pressure in the Eastern Cape(Berliner 2009).

Gomo encompasses a ～500 ha patch of indigenous Transkei mistbelt forest, and represents one of the more limited, larger remnant patches of Transkei mistbelt forest with a larger core area and higher biodiversity value on account of being less affected by edge effects. Environmentally, Gomo is located on a south-easterly slope,with an elevational gradient ranging from 850 to 1500 m above sea level, typical of Transkei mistbelt forests. Furthermore, pockets of commercially managed pine plantation occur along the forest boundary or nested within the forest itself, and a gravel road intersects the length of the forest, features commonly associated with Transkei mistbelt forests (Berliner 2009). Historically, Gomo has endured logging pre-1940, followed by subsistence harvesting in recent times, and is thus representative of the disturbance history of forests in the region. Socioeconomically, a number of rural communities occur less than 3 km from its boundary, characteristic of forests in the region, and is managed by the Department of Environment, Forestry, and Fisheries(DEFF)..

Study design

To investigate the effect of forest product harvesting on the bird community at the forest-scale, 16 circular plots(0.04 ha) were sampled within two distinct elevational zones (Fig. 1c): eight plots were located within the midelevational zone of the forest, close to the road (hereafter mid-zone),and eight plots were located towards the upper forest boundary within the high-elevational zone of the forest, further from the road (hereafter high- zone). Midzone plots occurred at a mean elevation of 1225 m above sea level, and were, on average 82 m away from the road,while high-zone plots had a mean elevation of 1473 m above sea level,and were on average 770 m away from the road. A mean elevational gradient of 250 m existed between mid- and high-zone plots. Within each elevational zone, a transect perpendicular to the forest slope, and representing similar slope conditions within each zone was identified, i.e. respective transects differed in elevation, but were similar with regards to variation in slope.Plots within each elevational zone were then randomly selected within a maximum distance of 50 m either side of each transect, and a minimum distance of 100 m away from the nearest plots. Plots within the mid-zone were more widely spread along the length of the transect given the more consistent slope conditions in this zone. Within the top-zone, a shorter transect of comparable slope conditions was identified, such that plots were located within a narrower band, and placed either side of the transect(Fig. 1c). Within each plot, harvest disturbances, habitat structure and the bird community were recorded.

Fig.1 Location of a)the Eastern Cape Province within Africa,b)Gomo forest within the Eastern Cape,and c)plot location(white circles)within Gomo forest,with the intersecting road shown in brown

Data collection Habitat variables

At each plot, habitat variables were recorded within three nested circular plots: the largest plot was 0.2 ha(radius of 25 m), within which two smaller plots of 0.04 ha (radius of 11.3 m), and 0.01 ha (radius of 5.6 m) were nested. In the 0.2 ha plot, all standing dead trees (henceforth, snags) were recorded by diameter at 1.3 m above the ground (DBH), and cause of death (i.e. natural or due to bark harvesting).In the 0.04 ha plot,the following variables were recorded: DBH of all living stems (＞5 cm DBH),diameter of all harvested trees(henceforth stumps),percentage canopy cover,mean canopy height,percentage coverage of bare ground;leaf litter; grass cover; and herbaceous cover, and foliage density at 0-0.5 m; 0.5-1 m; 1-2 m; 2-5 m; 5-10 m and 10-20 m. Foliage density at each height class was estimated using a telescoping pole eight meters long and marked at each height interval. The pole was sequentially set-up at eight evenly spaced points 11.3 m from the plot centre (i.e. along the 0.04 ha circular plot boundary)and visual estimates of foliage density(as a percentage) at each height class were made from the plot centre. A rangefinder was used to assist with estimates of foliage density beyond the length of the telescoping pole,as well as to estimate mean canopy height at each plot.Lastly,the number of Ocotea bullata stems(＞5 cm DBH)was recorded in each plot as this nationally endangered tree species was shown to be under heavy resource use pressure in Gomo in a previous linked study due to its durable wood harvested for poles, and its medicinal bark which is in high market demand.In the inner-most plot of 0.01 ha (radius of 5.6 m) stem density of saplings was recorded by counting all stems with diameter 1-5 cm.

Bird surveys

Bird surveys were conducted at plots (n=16) during the summer breeding season in the study region (November-December 2017). Non-fixed-radius point-counts(Blondel 1981) were conducted to sample birds . All birds calling over a 10 min period were recorded using a Song Meter SM4 acoustic recorder attached to a tree at a height of 1.5 m near the centre of each plot. Birds seen during the 10 min period were visually identified by JL in the field, and birds recorded on the Song Meter were audibly identified thereafter through playback of recordings by CC. Each site was surveyed three times during the morning period (sunrise +3 h), with repeated surveys conducted on different days. Surveys were conducted in alternating sequence so as to ensure that repeat surveys at each survey site were done at different times within the three-hour morning period. Bird surveys were consistently conducted on dry, still days.Any birds seen or heard which could not be confidently identified were not recorded. Presence/absence data at each plot were pooled, thereby determining species richness as the cumulative number of species recorded within a plot. Following identification, recorded bird species were classified into two groups based on their level of forest dependency, namely, forest-specialist or forest-generalists, based on Oatley (1989) and Hockey et al. (2005) (Additional file 4: Table S1). Forestspecialist species were defined as those that rely on forest resources to survive and reproduce (Oately 1989). Conversely, the forest generalist guild included species that are not, or only partly, dependent on forest resources and thus occur in forests as well as other habitats (Oatley 1989; Neuschulz et al. 2011).

Data analyses Habitat data

Mean DBH (all stems ＞5 cm DBH), stem density and basal area were calculated from tree diameters recorded within plots. Foliage density at each height class interval was calculated as the mean percent density from the eight separate estimates taken.

Comparing harvest rates across elevational zones

A bark harvesting index was assigned to each plot based on the proportion of trees dead due to bark harvesting(i.e. bark harvested snags). This was calculated by dividing the number of bark-harvested snags (standing dead trees ＞10 cm DBH) in a plot by the total number of stems (living and dead ＞10 cm DBH). Given that snags were measured within the 0.2 ha plot, and trees within the 0.04 ha plot, snag and tree abundances were standardized to abundance per hectare, and the intensity of bark harvesting per plot calculated as the overall proportion of bark-harvested snags per hectare:

Number of bark-harvested snags per ha/(total number of living +dead stems per ha).

A pole harvesting index was assigned to each plot based on the proportion of trees (diameter ＞5 cm)harvested per plot. This was calculated based on the accumulated harvestable stems (stumps plus standing stems ＞5 cm diameter) and used as an index of tree harvesting at each plot as follows:

Number stumps per plot/(number stumps + number stems)per plot.

Lastly, a composite Harvest Disturbance Index (HDI)was developed to score each plot according to its overall level of harvest disturbance at the forest-scale. The two measured harvest indices (bark harvest index and pole harvest index) were relativized by their respective maximum values recorded within the forest. Relativized scores were then summed at each plot such that an overall HDI score, ranging from 0 (indicating no harvest disturbance) to a maximum potential score of 2 (indicating the most harvest disturbance), was assigned to each plot. Calculated HDI scores and harvest rates of each forest product were not normally distributed (Shapiro-Wilks test, p ＜0.05 significance threshold). Therefore,Wilcoxon tests were used to compare HDI scores and harvest rates across sample forests.

Effects of harvesting on habitat structure

Harvest disturbance effects on measured structural variables were investigated using linear mixed models(LMMs). The mixed-modelling approach accounted for the nested study design, with elevational zone set as a random effect to account for plots being nested within two distinct, spatially clustered groups. Separate LMMs were used to assess the response of each measured habitat feature to harvesting, with the Harvest Disturbance Index score per plot set as the explanatory variable,with Gaussian errors, using the ‘lme’ function of the nlme package in R version 3.4.3 (Pinheiro et al. 2011; R Core Development Team 2017). Assumptions of normality and homogeneity were assessed using graphical outputs of models. Response variables that were measured as percentages were logittransformed to improve the model assumption of normality. Where variance heterogeneity affected the estimation of harvesting effects, models were run with a constant variance function structure using the ‘varIdent’ function in the nlme package. To quantify the goodness-of-fit for each model, the MuMin package was used to evaluate marginal R2and conditional R2which can be respectively interpreted as the variance explained by the fixed effects only; and by both fixed and random effects (Nakagawa and Schielzeth 2013).

Comparing harvest-mediated habitat dissimilarity across elevational zones

Given that the range of harvest activities, in terms of both nature and extent, was expected to differ across zones due to their varying proximity to the road, it was expected that harvest activities would result in different levels of habitat dissimilarity, i.e. heterogeneity, across elevational zones. This was tested by creating a distance matrix based on scaled Euclidean distances of all habitat variables shown to be significantly affected by harvest activities (derived from the outcome of LMM analyses described above), and calculating the mean pairwise distance of plots within each elevational zone respectively. Each plot was thus assigned a habitat heterogeneity score based on its mean pairwise dissimilarity from other plots within the same zone. Mean harvest-mediated habitat heterogeneity was then compared between mid- and high-zone plots using a non-parametric Wilcoxon test.

Comparing bird species richness and beta-diversity across elevational zones

Bird analyses were conducted from a total of 15 plots. A single high-zone plot was removed from analyses as it had a very small number of observations relative to other plots,on account of the presence of loud calling insects during one of the surveys compromising the quality of the recording.Mean species richness was compared across elevational zones for the forest-specialist and forest-generalist group using T-tests and non-parametric Wilcoxon tests,based on the distribution of the data. To assess variation in taxonomic composition between elevational zones, analysis of similarity (ANOSIM) and non-metric multi-dimensional scaling (nMDS) were used, with each bird group analysed respectively. ANOSIM tested for statistically significant differences in species composition between mid-zone and high-zone plots (S?rensen’s presence/absence index), while nMDS was applied as a visual aid to interpretation of how plots differed between elevational zones. Measurements of variation in species composition across sites (i.e. betadiversity) were quantified by the dissimilarity in species composition using presence/absence species data, and was assessed within each elevational zone. Comparison of within-zone beta-diversity was based on the quantitative S?rensen’s index of pairwise dissimilarity: ?Sor=(b+c)/(2a+b+c), where a is the number of species common to both sites,b is the number of species that occur in the first site but not the second site, and c is the number of species that occur in the second site but not the first site. This measure incorporates change due to species loss, including nestedness of sites where one set of species is a subset of another set, and change due to replacement of one set of species by another, i.e. species turnover (Koleff et al. 2003;Baselga 2010). These processes were distinguished in this study by partitioning total beta-diversity (?Sor) into contributions by turnover (Simpson’s dissimilarity: ?Sim), and nestedness-resultant dissimilarity (?nes), following Baselga(2010). ?Simdescribes turnover without the influence of richness gradients, and ?nesis derived from the difference between ?Sorand ?Sim, accounting for the nestedness component of beta-diversity. Each index varies between 0 and 1, with lower values indicating a greater proportion of shared species richness,and larger values indicating greater dissimilarity between locations. Overall dissimilarity (based on S?rensen’s index),and the proportion of dissimilarity attributed to turnover (based on Simpson’s dissimilarity) and species loss (based on nestedness-resultant dissimilarity),were calculated between every pair of plots within each zone respectively to compare within-zone variability in species composition,and the processes driving this for each bird group separately. Levels of species turnover and species loss were compared across elevation zones for each bird group using t-tests or a Wilcoxon rank sum test,depending on the distribution of the data.

Assessing the effect of habitat modification on bird species richness

Linear mixed models (LMMs) were used to assess the effect of harvest-mediated habitat heterogeneity on bird species richness, with elevational zone again included as a random effect to account for plots being nested within two distinct, spatially clustered groups. Species richness values for each group were included as the response variable in separate models, and harvest-mediated habitat heterogeneity scores per plot were included as the explanatory variable. Response variables were modelled using LMMs with Gaussian errors,using the‘lme’function of the nlme package in R version 3.4.3 (Pinheiro et al. 2011; R Core Development Team 2017). Species richness values were logtransformed to improve the model assumption of normality, and model assumptions were assessed using graphical outputs of models(Zuur et al.2010).

Assessing beta-diversity associations with elevation and habitat modification respectively

Beta-diversity was examined with respect to difference in elevational zone and harvest-mediated habitat heterogeneity. As the distance matrices of elevational zone and harvest-mediated habitat heterogeneity were inter-correlated, two partial Mantel tests were used to test: i) whether beta-diversity was related to differences in elevational zone (i.e. mid-zone vs. top-zone),controlling for differences in harvest-mediated habitat heterogeneity;and ii)whether beta-diversity was related to differences in harvest-mediated habitat heterogeneity,controlling for elevational zone. Beta-diversity was based on i) overall S?rensen’s dissimilarity values, ii) Simpson’s dissimilarity (i.e. dissimilarity attributed to turnover), and iii)nestedness-resultant dissimilarity values(i.e.dissimilarity attributed to species loss), whereas distance matrices for elevational zone and harvest-modified habitat structure where based on Euclidean distances, scaled in the case of the habitat variables. Separate partial Mantel Tests were run to examine correlations between environmental variation and beta-diversity in the forest-specialist guild and forest-generalist guild respectively.

Results

Harvest rates across elevational zones

Harvest rates of medicinal bark, although higher in the mid-zone, did not differ between elevational zones (midzone: 0.06±0.04 vs. high-zone 0.03±0.04, W=18.5, p=0.17; Fig. 2a). Conversely, pole harvest rates were higher in the mid-zone, where pole harvesting was recorded in all but one plot, compared to the high-zone, where pole harvesting was recorded only in a single plot (mid-zone:0.16±0.13 vs. high-zone 0.006±0.02, W=6, p=0.004;Fig. 2b). Consequently, overall harvest disturbance index(HDI) was higher in the mid-zone compared to the high-zone (mid-zone: 0.98±0.51 vs. high-zone 0.27±0.35, W=3, p=0.001; Fig. 3c).

Effects of harvesting on habitat structure

At the forest-scale, harvest disturbances based on overall Harvest Disturbance Index scores (HDI) negatively affected canopy cover, canopy foliage density (5-10 m),abundance of Ocotea bullata stems (＞5 cm DBH), overall abundance of tree stems (＞5 cm DBH) and herb cover. Conversely, woody debris cover, understory foliage density (0-2 m) and snag abundance were positively associated with HDI (Fig. 3; Additional file 5:Table S2).Based on this subset of eight habitat variables significantly affected by harvest disturbances, calculated habitat dissimilarity scores per plot, i.e. harvest-mediated habitat heterogeneity, was higher in the mid-zone compared to the high-zone (mid-zone: 3.38±0.41 vs. highzone 2.25±0.51, t=-4.91, p=0.0003; Fig. 4).

Fig.2 Harvest rates compared across two elevational zones in Gomo forest indicating a)proportion of standing stems(dbh ＞10 cm)dead due to medicinal bark harvesting, b)proportion of stems(dbh ＞5 cm)harvested for poles,and c)overall Harvest Disturbance Index,based on the combined rate of pole and bark harvesting recorded in a plot.Asterisks indicates significant differences across zones

Fig.3 Response of habitat variables to overall Harvest Disturbance Index scores(i.e.combined bark and pole harvest rates per plot).Relationships shown are significant(p ＜0.05),derived from linear mixed models with elevational zone included as a random effect(Table S2)

Bird species richness,occurrence and beta-diversity

Overall, 34 species were recorded during bird surveys conducted at mid- and high-zone plots (18±3 species per survey site, mean±SD, range 13-22). Thirty-three species were recorded in the mid-zone (20±2 species per survey site, mean±SD, range 17-22), whereas 26 species were recorded in the high-zone (16±2) species per survey site; mean±SD, range 13-18). Based on Chao2 estimator of true species richness, the sampling effort yielded 98% and 95% of the “true” species present in the mid- and high-zones respectively. Sample-based species accumulation curves based on Coleman’s method were asymptotic for both zones, further indicating that the sampling effort was sufficient to represent true species richness present in each zone (Additional file 1: Fig. S1). Mean species richness of the forestgeneralist group was higher in the mid- zone, whereas forest-specialist species richness did not differ significantly across elevational zones (Table 1).

Three-quarters of the 34 recorded species were observed in both mid- and high-zones, indicating forestwide distributions of most species. Nonetheless, bird communities of both groups differed significantly across elevational zones in taxonomic composition (ANOSIM,Jaccard’s presence/absence index, 9999 permutations,Forest-specialist group: Global R=0.34, p ＜0.01; Forestgeneralist group: Global R=0.57, p ＜0.01). In support of this, nMDS analyses showed clustering of survey sites according to elevation, based on species presence/absence (Fig. 5). Of the 34 species recorded, only one species, Apalis thoracica (Bar-throated apalis), was absent from the mid-zone, while eight species present in the mid-zone were not recorded at the high-zone, comprising four forest-specialist species: Chrysococcyx cupreus(African Emerald cuckoo), Poicephalus robustus (Cape parrot), Apaloderma narina (Narina trogon) and Bycanistes bucinator (Trumpeter hornbill); and four forestgeneralist species: Dicrurus adsimilis (Fork-taileddrongo), Campethera notata(Knysna woodpecker),Pogoniulus pusillus(Red-fronted tinkerbird)and Apalis flavida(Yellow-breasted apalis) (Additional file 4: Table S1).Interestingly, five of the eight species not recorded in the top-zone were cavity-nesting species (Poicephalus robustus, Campethera notate, Apaloderma narina, Pogoniulus pusillus and Bycanistes bucinator).

Pairwise beta-diversity and the relative contribution of the two mechanisms driving this (i.e. species turnover and species loss)differed across groups and elevational zones(Fig.6).Across both zones, beta-diversity was greater in the forestspecialist community compared to the forest-generalist community. Within the mid-zone, species turnover was the dominant driver of beta-diversity,close to four times that attributed to species loss for both forest-specialist and forestgeneralist communities.In the high-zone,the relative contribution of species loss to beta-diversity increased for both groups. However, species turnover remained the dominant driver of beta-diversity in the forest-specialist community in the high-zone, while species loss became the dominant driver of beta-diversity in the forest-generalist community.Consequently, for the forest-generalist community, the amount of beta-diversity attributed to species turnover was higher in the mid-zone compared to the high-zone (W=402.5; p=0.03), while the relative contribution of species loss was greater in the high-zone compared to the mid-zone(W=144;p=0.002).Conversely,the relative contribution of species turnover (t=1.21; df=44.07; p=0.23) and species loss (W=218.5; p=0.13) did not differ across zones for the forest-specialist group.

Effect of harvest-mediated habitat heterogeneity on bird species richness

Forest-specialist species richness declined in response to harvest-mediated habitat heterogeneity (?=-0.39±0.15;p=0.04), while forest-generalist species richness was positively associated with harvest-mediated habitat heterogeneity (?=0.15±0.06; p=0.03; Fig. 6; Additional file 6:Table S3). Harvest-mediated habitat heterogeneity explained 20% and 32% of the variation in forestspecialist and generalist species richness respectively(Additional file 6: Table S3).

Role of elevational zone and harvest-mediated habitat modification in explaining differences in bird composition

Beta-diversity based on S?rensen’s index was positively correlated with difference in elevational zone when controlling for harvest-mediated habitat heterogeneity in the case of both groups (Table 2). This indicates that there were effects of elevation on avifaunal composition in addition to those caused by variation in harvest disturbances, reflecting results of the ANOSIM and nMDS plots (Fig. 5). Beta-diversity attributed to species turnover (Simpson’s dissimilarity index) was not correlated with difference in elevational zone for both groups,while beta-diversity attributed to species loss (i.e. nestednessresultant index) was positively correlated with elevational zone in the forest-generalist group when controlling for harvest-mediated habitat heterogeneity, i.e. there were effects of elevation on patterns of generalist species loss in addition to those caused by harvest disturbances.There was a positive correlation between S?rensen’s dissimilarity and harvest-mediated habitat heterogeneity when controlling for elevational zone in the forestgeneralist group, indicating that variation in harvest disturbances accounted for compositional variation in this group, but not in the forest-specialist group. However,there was a positive correlation between beta-diversity attributed to species turnover and harvest-mediated habitat heterogeneity in the forest-specialist group when controlling for elevation, indicating that variation in harvestdisturbances accounted for the amount of beta-diversity attributed to species turnover in this group(Table 2).

Table 1 Comparison of species richness across two elevational zones in Gomo forest.Mean± SD species richness is shown for the forest-dependent and forest generalist assemblage in each zone.P-values are Bonferroni-adjusted to account for multiple testing.Text in bold indicates a statistically significant difference across zones

Fig.5 Non-metric multidimensional scaling(nMDS)representing clustering of high-zone(red)and mid-zone(blue)plots by a)forest-specialist species,and b)forest-generalist species,based on species presence/absence(S?rensen’s dissimilarity,95%ellipses)in Gomo forest

Discussion

Habitat modification caused by the harvesting of poles and medicinal bark variably affected forest-generalist and -specialist avifaunal communities in Gomo. Specifically, the variable nature, extent and spatial distribution of harvest activities increased habitat heterogeneity, positively affecting forest-generalist species richness,but negatively affecting forest-specialist richness (Fig. 7). Elevation affected spatial patterns of harvest disturbances and avifaunal diversity in Gomo, despite the relatively small gradient (＜300 m) investigated: harvest rates were higher,and harvest activities more varied at lower elevations,resulting in greater harvest-mediated habitat heterogeneity; and avifaunal species richness, community composition and the mechanisms driving beta-diversity differed across elevational zones. By controlling for elevation, we show that harvest-mediated habitat heterogeneity was positively associated with overall beta-diversity of the forest-generalist community, and the amount of betadiversity attributed to species turnover in the forestspecialist community (Table 2). By controlling for variation in harvest-mediated habitat heterogeneity,we show that natural variation associated with changes in elevation also affected patterns of overall beta-diversity in forest-generalist and -specialist communities. Similarly,beta-diversity attributed to species loss in the forestgeneralist community was positively correlated with changes in elevation. Thus, while overall beta-diversity of the forest-generalist group was affected by elevation and habitat modification, mechanisms driving generalist betadiversity (i.e. species loss) were affected by elevation. On the other hand, overall beta-diversity of the forestspecialist group was affected by elevation only, while the mechanisms driving specialist beta-diversity (i.e. species turnover) were affected by harvest disturbances. Importantly, these findings show that avifaunal communities in Gomo are shaped by natural environmental gradients associated with changes in elevation and human-mediated disturbance gradients associated with harvesting activities,with responses dependent on species’ level of habitat specialization.

Fig.6 Overall pairwise beta-diversity represented as that contributed by species turnover(?sim: dark grey)and species loss(?nes:light grey)within elevational zones in Gomo across the two bird groups:forest-specialist species,and forest-generalist species.Mean dissimilarity values are shown,with higher values indicating greater dissimilarity

Table 2 Partial mantel correlations testing i) the association between avifaunal beta-diversity and difference in elevational zone,while controlling for harvest-mediated habitat heterogeneity, and ii) the association between avifaunal beta-diversity and harvestmediated habitat heterogeneity,while controlling for difference in elevation for plots sampled in Gomo forest, using presence/absence bird data to derive indices of dissimilarity: S?rensen’s index,Simpson’s index and Nestedness-resultant index

Species richness

This study demonstrates that forest-scale habitat heterogeneity is an important predictor of avian species richness,but that the direction of response is dependent on species’level of habitat specialisation. Specifically, the positive response of forest-generalist species to human-meditated habitat heterogeneity reflects the well-established finding that ecological generalists are more likely to benefit from unstable, heterogeneous environments, and thus habitat modification, given their ability to exploit a wide range of habitat conditions (McKinney and Lockwood 1999). Similarly, the decline in forest-specialist species richness in response to harvest-mediated habitat heterogeneity shown in this study can be explained by niche theory which predicts that habitat specialists should benefit from more stable, homogenous environments, and thus be negatively affected by humanmediated habitat modification (Soh et al. 2006; Devictor et al. 2008a; Clavel et al. 2011). The habitat heterogeneity hypothesis (MacArthur and MacArthur 1961) stipulates that resources and niches increase with spatial heterogeneity (Pianka 1972; Bazzaz 1975), which should in turn facilitate the co-occurrence of species (Jeltsch et al. 1999; Palmer 2003) and provide habitat for species with multiple resource requirements (Perkins et al.2000), thereby increasing species richness (Terborgh 1971). In the current study, this applies to forestgeneralists but not to specialist species. Similarly, Stirnemann et al. (2014), showed that increases in habitat heterogeneity did not always result in increased avifaunal species richness in a temperate forest in Australia, and that species’ responses to habitat heterogeneity depended on their ecology. Stirnemann et al.(2010) explains this in terms of increases in niches or resources through increased heterogeneity leading to increased competition for resources between species,resulting in species turnover rather than opportunities for additional species to establish. This may explain our finding that harvest-mediated habitat heterogeneity was positively correlated with species turnover in the forestspecialist group, but not the forest-generalist group: it may provide increased opportunities for forestgeneralist species, thereby increasing species richness,but result in increased competition between forestspecialist species, driving increased species turnover,and an overall loss of species richness.

Fig.7 Response of a)forest-specialist species richness and b) forest-generalist species richness to harvest-mediated habitat heterogeneity, based on linear mixed models(LMMs),controlling for variation in elevational zone

It is also important to consider that increased heterogeneity does not always result in increased niches or resources becoming available (Stirnemann et al. 2014). In the current study, we measured how variable plots within elevation zones were from one another based on a number of harvest-modified variables relating to the amount of cover of different habitat features; and the number of living trees; dead trees; and Ocotea bullata stems. This measure thus combined aspects of spatial habitat heterogeneity and cover, which reflect different needs for birds: cover may relate to amounts of resources while heterogeneity relates to the spatial arrangement of those resources.The combination of these factors may influence whether heterogeneity results in increased species richness: Stirnemann et al. (2014) showed that species richness increased where both cover and spatial heterogeneity were high, but declined where heterogeneity was high but cover was low.Given that the measure of habitat heterogeneity used in the current study combined spatial heterogeneity and cover, it was not possible to investigate how these factors separately affected patterns of species richness. However, while harvest activities did increase spatial habitat heterogeneity across plots, they were also shown to decrease canopy cover, canopy foliage density and herb cover. This may offer an alternative explanation for the respective responses to harvest-mediated habitat heterogeneity by forest-generalist and -specialist species,in that forest-specialists may be unable to benefit from increases in spatial heterogeneity when cover of certain habitat features, such as the canopy and herb layer, have been reduced in the process;while habitat-generalists may be more resilient to open-canopy conditions and thus better able to take advantage of spatial habitat heterogeneity.

Species composition and beta-diversity

While species richness is an important measure to assess community-level responses to disturbance, it provides little information on concomitant changes in species composition. Thus, in addition to species richness, this study assessed beta-diversity and the mechanisms driving this across elevational zones for each avifaunal group separately. Specifically, findings of this study provide insight into the separate effects of harvest-mediated habitat heterogeneity and elevation on patterns of avifaunal community composition at the forest-scale respectively, and how this varied across forest-generalist and forest-specialist species.

Species composition of the forest-generalist and forestspecialist group differed across elevational zones (Fig. 5;Additional file 2: Fig. S2 and Additional file 3: S3). However, the effect of elevation on bird species composition was not fully accounted for by accompanying variation in harvest-mediated habitat heterogeneity (Table 2). Thus,despite the relatively small elevational gradient assessed(＜300 m), variation in environmental conditions other than disturbance,such as moisture availability,soil characteristics and floristic composition, likely affected patterns of beta-diversity across the elevational gradient in Gomo.For example, tree species richness and composition were observed to change across elevational zones,with tree species richness declining with increasing elevation (personal observation).Changes in avifaunal composition across elevational zones may thus be associated with changes in tree species richness,as shown by Jankowski et al.(2013)along an Andean elevational gradient.

While elevation was an important driver of overall forest-scale beta-diversity of forest-specialist and -generalist species, evidence of a particular mechanism driving this was only shown in the case of the forest-generalist community (Table 2). Specifically, the amount of betadiversity attributed to species loss increased with elevation for forest-generalists in Gomo (Table 2). Jankowski et al.(2013) found a similar trend in the avifaunal community along an elevation gradient in the Andes. Our finding indicates that the generalist community in the high-zone was largely a subset of the species occurring in the midzone, supported by the lower species richness of forestgeneralists recorded in the high-zone (Table 1). Thus,natural environmental gradients associated with elevation affected forest-scale beta-diversity of the forest-generalist through a process of species loss, with the high-zone representing an impoverished zone. Conversely, the role of elevation in driving forest-scale beta-diversity of the forest-specialist community was not clearly attributed to either species turnover or species loss when controlling for variation in harvest-mediated habitat heterogeneity.Furthermore, there was considerably higher overlap in taxonomic species composition across zones in the forestspecialist community compared to the forest-generalist community(Fig. 5),reflecting the stronger correlation between elevation and forest-generalist beta-diversity (0.30)compared to that observed between elevation forestspecialist beta-diversity (0.19; Table 2). This suggests that the forest-specialist community was less strongly influenced by natural changes along the elevation gradient at Gomo compared to the forest-generalist community. On the other hand, human-mediated disturbance gradients,namely variation in harvest-mediated habitat heterogeneity, affected the mechanisms driving beta-diversity of the forest-specialist community. Specifically, forest-specialist species turnover was positively correlated with harvestmediated habitat heterogeneity when controlling for elevation (Table 2). Similarly, Palmeirim et al. (2017) showed that the contribution of species turnover to beta-diversity of lizard and amphibian communities increased with human disturbance in the Neotropics. Thus, humanmediated habitat heterogeneity in Gomo was an important driver of forest-specialist species turnover, indicating that harvest disturbances operate as environmental filters for specialist species,but not generalists(Baselga 2009).

Results of this study thus show that avifaunal communities in Gomo are structured by harvest disturbances and elevation, despite the relatively small gradient assessed (＜300 m). Similarly, a study conducted in the Andes showed that bird communities along far greater elevation gradients were shaped by both changes in elevation and disturbance (Montano-Centellas and Garitano-Zavala 2015). Furthermore, Peters et al.(2019) recently showed that variation in species richness and composition of multiple taxa along an elevational gradient of Mount Kilimanjaro was explained by the interaction of land-use intensity and climate, rather than by single drivers. Patterns of Afrotemperate biodiversity along elevational gradients in the Eastern Cape are similarly likely affected by the interaction of multiple factors. Further research is thus needed to assess the drivers of Afromontane bird communities along elevational gradients,in particular,the interactive effect of harvest disturbances and varying environmental conditions along elevational gradients.

Conclusion and conservation implications

The combined effect of harvesting medicinal bark and poles in Gomo affected avifaunal communities at the forest-scale, mediated by harvest-mediated increases in habitat heterogeneity. Specifically, while forest-generalist species richness and overall beta-diversity were positively affected by harvest-mediated increases in habitat heterogeneity, forest-specialist species richness was negatively affected, but species turnover positively affected by harvest-mediated habitat heterogeneity. These results suggest that harvest disturbances and concomitant habitat modifications provided more niches and resources,allowing more opportunities for habitat generalist species, but not for forest-specialist species. Furthermore,harvest-mediated habitat modifications acted as environmental filters for specialist species but not generalists.Findings of this study thus indicate the importance of using different biodiversity metrics when assessing forest biodiversity responses to habitat disturbance. Specifically, the use of a single measure of species richness is cautioned against, as an increase in forest-generalist species in response to disturbance may mask the loss of forest-specialists. However, conclusions drawn from this study are to be considered with caution given that a single spatial scale was considered, thereby limiting insight into disparate patterns that may be revealed at larger scales (Hill and Hamer 2004; Rocha et al. 2015;Morante-Filho et al. 2016). Thus, while findings show that current unregulated rates of harvesting increase habitat heterogeneity at the forest-scale, particularly in accessible areas at lower elevations and/or close to forest roads, with concomitant positive impacts on forest generalist bird communities and negative impacts on forestspecialist species, further research is needed to assess harvest impacts on forest habitats and biodiversity at broader spatial scales. Nonetheless, this study provides previously unexplored, yet important insights into the role of elevation and harvest disturbance in driving spatial patterns of avifaunal diversity in temperate forests of South Africa, specifically the different mechanisms driving beta-diversity at the forest-scale, and how these vary across the relatively small elevation gradient of this Afromontane forest. Results show that elevation has a strong effect on spatial patterns of harvesting patterns,habitat structure and avifaunal communities, despite the small elevational gradient examined.

Supplementary information

Supplementary information accompanies this paper at https://doi.org/10.1186/s40663-019-0207-x.

Additional file 1: Figure S1. Sample-based species accumulation curves based on Coleman’s method from two elevational zones in Gomo forest.

Additional file 2: Figure S2. nMDS ordination of forest-specialist species.

Additional file 3: Figure S3. nMDS ordination of forest-generalist species.

Additional file 4: Table S1.List of species recorded at Gomo over the summer sampling period.

Additional file 5: Table S2.Response of habitat variables to Harvest Disturbance Index scores,derived from linear mixed models.

Additional file 6: Table S3.Response of species richness to harvestmediated habitat heterogeneity derived from linear mixed models.

Acknowledgements

The authors wish to thank Welile Kedama of the Department of Agriculture,Forestry and Fisheries (DAFF)for granting permission to conduct research within state forests of the Eastern Cape, and to all Forestry personnel that assisted with fieldwork.C.J. Geldenhuys is thanked for meaningful discussion that helped guide this study. We are grateful to A. Wannenburgh for providing input on the first draft of this manuscript, and for the study site map, which greatly improved the manuscript.R. Duker is thanked for assistance in developing data collection methodology.

Ethics approval and content to participate

Permission to conduct this study in state forests was approved by the former Department of Agriculture, Forestry and Fisheries (Licence no.: WIFM 04-2016) under section 23(1)(K) of the National Forest Act,1998.

Authors’contribution

JL designed the study, collected and analysed data, and drafted the manuscript. JCC identified bird calls from recordings. MIC conceptualised the study, contributed to data interpretation and wrote the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the National Research Foundation (NRF), South Africa (FBIP 98871). The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

Availability of data and materials

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

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1Department of Botany and Zoology, Stellenbosch University, Private bag X1,Matieland 7602, South Africa.2Wild Bird Trust, Cape Parrot Project, 20 Loch Avenue, Parktown 2193, South Africa.

Received: 19 July 2019 Accepted: 28 October 2019