Why a landscape view is important: nearby urban and agricultural land affects bird abundances in protected areas

Study area

Our study area consisted of a heterogeneous landscape (a square with coordinates at the NW corner: 25S 27E and SE corner: 27S 29E) that included the greater Gauteng Province of South Africa. It consisted of a rich mix of urban and other heavily human-modified land-use types, as well as protected areas, and the study area included the cities of Pretoria and Johannesburg (Fig. 1), which are two of the most densely populated cities in South Africa (Statistics South Africa, 2012).

The study area was approximately 35,000 km², and comprised of eight land-use types (Fig. 1): mines (0.80% of total land use); plantations (0.32%); waterbodies (2.80%); degraded (2.54%); protected area (6.40%); urban (8.13%); agriculture (28.71%); and natural land (50.30%). Here, natural land refers to land that is not primarily used for any of the other aforementioned land uses. Therefore, in addition to representing the naturally occurring vegetation (grass and trees), it can also represent small holdings, open plots alongside roads or between agricultural area, and recreational land uses (such as sports-fields, parks and lawns). Land-use data were provided by the South African National Biodiversity Institute (South African National Biodiversity Institute (SANBI), 2009) at a 30 m × 30 m resolution.

Urban, agricultural and natural land made up approximately 87% of the land-use cover over the study area, and are known to be influential in affecting bird distributions (Brandon, Redford & Sanderson, 1998; Knapp et al., 2008; Santos et al., 2008; Zhang et al., 2011). We therefore exclusively examined these land-use types in addition to protected area in this analysis.

Species detection/non-detection data

We used bird detection data from the second Southern African Bird Atlas Project (SABAP 2) which started in June 2007 (Harebottle et al., 2007) and was on-going in 2018, when we performed the analysis. SABAP 2 is a citizen science project whereby registered volunteers submit checklists of birds they observed during a fixed time period within a pre-defined area called a pentad, which is 5′ × 5′ in dimension (unit is arcminutes; approximately 61 km² in area). Volunteers must have spent at least 2 hours, but not more than five days searching for birds within each pentad. Only the presence of bird species is recorded per pentad, not the number of birds seen. Observers were asked to sample all habitats within the pentad. Unusual records were scrutinized by a vetting committee, who either accepted or rejected the record based on supporting information (Harebottle et al., 2007). Every checklist constitutes a detection (if the species was recorded) or non-detection (if the species was not recorded) for all species that occur in the area. We treat each checklist as an independent survey to a particular grid cell.

We considered only common, resident bird species within the study area, and omitted any nomadic, alien, and migratory species, totalling 198 species. We included bird atlas data that were collected and submitted to the project between the beginning of January 2014 and end of December 2015. The years 2014 and 2015 had enough data to support robust data analyses, whilst being short enough a time period to assume abundance did not change markedly in this period (see “Discussion”). Because pentads within the study area were not surveyed the same number of times, like Broms et al. (2014), we randomly selected 100 checklists for pentads that had more than 100 checklists (Fig. A1, Appendix). The study area covers 576 pentads (a 24 pentad by 24 pentad grid), for which 10,400 checklists were submitted at an average of approximately 18 checklists per pentad (min. 1 and max. 468)

Each species was assigned to a guild based on information in Hockey, Dean & Ryan (2005). The definition of a guild is based on a species’ primary food source, and its primary foraging mode. We distinguished between seven guilds, namely: frugivores (species that primarily consume fleshy fruit, totalling nine species); gleaners (species that primarily consume insects and other invertebrates caught off plants, totalling 30 species); granivores (species that primarily consume seeds and grains, totalling 48 species); ground-feeders (species that primarily consume insects and invertebrates caught off the ground, totalling 62 species); hawkers (species that primarily consume insects and other invertebrates caught in the air, totalling 11 species); predators (birds of prey, species that primarily consume the flesh of vertebrates, totalling 19 species), and vegivores (vegetative herbivores; species that primarily consume vegetative parts of plants, totalling 19 species). Table A2 (Appendix) shows the species that constitute each guild. In total, we considered 198 species.

Analyses

To model the abundance of each species per pentad, we used an extension of traditional occupancy models, known as the Royle-Nichols model of abundance (Royle & Nichols, 2003). Briefly, occupancy models are a class of models which use detection/non-detection data to estimate the probability that a species occurs within a specified area (a pentad in this case). The Royle-Nichols model infers abundance based on detection/non-detection data. These models account for the fact that most species are not observed perfectly in each habitat in which they occur (MacKenzie & Kendall, 2002; Pellet & Schmidt, 2005). Failure to account for non-detection may bias parameter estimates (Boulinier et al., 1998; MacKenzie et al., 2006; Nichols et al., 1998).

We tested for collinearity among the covariates using the Variance Inflation Factor. The Variance Inflation Factors for all our covariates were below 5 (Table A1, Appendix), suggesting multicollinearity wasn’t a problem in our dataset (Kock & Lynn, 2012; Montgomery & Peck, 1992).

Interpretations of model beta coefficients

In the absence of urban and agricultural area, the relationship between the proportion of protected area and bird abundance in a pentad is represented by ${\beta }_1$ (protection–abundance relationship). Our primary aim was to examine how this relationship changed with increasing proportions of urban and agricultural area in the same pentad. Therefore, the model parameters ${\beta }_4$ and ${\beta }_5$ (Eq. (5)), which estimate the effects of the interactions between the protected area and urban ( ${\beta }_4$) or agricultural ( ${\beta }_5$) area within the same pentad were of most interest. They indicate the degree to which the slope for the linear protection–abundance relationship changed when the proportion of urban ( ${\beta }_4$) or agricultural ( ${\beta }_5$) area within the same pentad changes. Species with a positive ${\beta }_4$ and ${\beta }_5$ value indicate that the slope of the protection–abundance relationship increases (i.e. becomes more positive) as the amount of urban or agricultural area increases within the pentad, meaning that the effect of protected area on the abundance of birds becomes more positive when urban or agricultural area neighbours protected area. The opposite is true for negative ${\beta }_4$ and ${\beta }_5$ values. We further examined variation in ${\beta }_4$ and ${\beta }_5$ through guilds using simple data aggregation, and a hierarchical Bayesian analysis (see sections below).

For what percentage of species does the protection–abundance relationship increase or decrease with increasing proportions of urban and agricultural land in the same pentad?

In response to increasing proportions of urban or agricultural area nearby protected area, the slope of the protection–abundance relationship may increase (a significantly steeper slope describing the protection–abundance relationship), decrease (a significantly less steep slope), remain the same (a slope with no significant change), or may even change sign completely (change from a positive slope to a negative slope, or vice versa). The type of change in the protection–abundance relationship for each species is indicated by the estimates for the interaction coefficients ( ${\beta }_4$ and ${\beta }_5$), specified in Eq. (5). Interpreting the estimates for these coefficients provides a good understanding of how bird abundances are predicted to change within protected area with increasing proportions of agricultural and urban area near protected area. However, a more thorough understanding of how this occurs is gained from interpreting these interaction coefficients with the estimates for the other land-use type covariates that pertain to land-use types ( ${\beta }_0-$ ${\beta }_3$, main effects in Eq. (5)). Thus, we used the main and interaction effects that pertain to land-use types ( ${\beta }_0-$ ${\beta }_5$), and examined eight scenarios of the way in which urban and agricultural area near protected area could potentially influence the relationship between bird abundance and the proportion of protected area (Figs. 2A–2H). There are, of course, many more hypothetical scenarios that could be considered, but the ones considered in this study are those that best align with the ecological hypotheses presented earlier. Distinguishing between these scenarios allows for a better understanding of the overall fitted relationships of the main and interaction effects in Eq. (5) as estimated by the model.

Figure 2 is divided into two columns which show positive interactions (left-hand column) and negative interactions (right-hand column). With positive interactions, the protection–abundance relationship becomes more positive as urban or agricultural area within the same pentad increase, as illustrated by the dotted line having a more positive slope than the solid one. With negative interactions, the protection–abundance relationship becomes more negative as urban or agricultural area within the same pentad increase, as illustrated by the dotted line having a more negative slope than the solid line.

In scenarios where pentads are without protected area, local abundance is a function of the main effects of urban or agricultural area, illustrated by examining values of the y-intercepts, and relative to natural land. The intercept for the solid line indicates the average abundance for a land-use type scenario of 0% protected area, and 100% natural land. The intercept for the dotted line indicates the average abundance for a land-use type scenario of 0% protected area, and 50% either urban or agricultural area, and 50% natural land. Scenarios A, C, E and G show situations where the species is more abundant in urban/agriculture area (dotted line y-intercept) than in natural land (solid line y-intercept—positive main effect); and in scenarios B, D, F and H, species are less abundant (negative main effect).

We fitted an abundance model to each species using Eqs. (1)–(5). Then, we categorised each species into scenarios A–H (Fig. 2), based on the way the species reacted to increases in urban area and agricultural area that surrounded protected areas. As a single guild was assigned to each species, we counted the frequency of species within each guild that fell into each interaction scenario (A–H). This indicated the manner in which the protection–abundance relationship for each guild is expected to change given increases in urban/agricultural land near protected area (assuming all other factors in that relationship remain constant). This analysis counts only the frequency with which species’ protection–abundance relationship is predicted to increase or decrease, rather than the magnitude of this change; this is considered in the next section below.

What is the average change in magnitude of the protection–abundance relationship with increasing proportions of urban and agricultural land in the same pentad?

We used a hierarchical Bayesian analysis to estimate the average change in the magnitude of the average protection–abundance relationship for each guild, along with associated credible intervals. This analysis used the species-specific mean and standard error estimates for the interaction terms $\hat{\beta }$₄ (urban) and $\hat{\beta }$₅ (agriculture) from Eq. (5) in the Royle-Nichols model of abundance to estimate a mean $\hat{\beta }$₄ and $\hat{\beta }$₅, and associated credible intervals, for each guild. The basic structure of this model was similar to a linear mixed-effects model, with guild as a random factor and normally distributed errors. However, instead of treating the mean interaction estimates as if they were observed values, we modelled them as coming from a normal distribution using the means and standard errors as estimated by the Royle-Nichols model of abundance. This approach is similar to the analysis described in McCarthy & Masters (2005) (see also Lloyd et al., 2014). We used non-informative priors for the mean interaction response per guild. We implemented this in the program to do still JAGS with 50,000 iterations and 25,000 burn in and 3 MCMC chains. The Gelman-Rubin diagnostic indicated that this model converged, and all R-hat values were <1.01. The JAGS code for this model is provided in the Appendix (Model A1).

Royle-Nichols abundance models

For what percentage of species does the protection–abundance relationship increase or decrease with increasing proportions of urban and agricultural land in the same pentad?

The protection–abundance relationship became more positive, and hence the overall abundance per pentad was expected to increase for 58% of all species as urban area near protected areas increased in proportion (Table 2A, row 8). The protection–abundance relationship became more positive for 49% of all species as agricultural area near protected area increased in proportion (Table 2B, row 8).

Of the positive interaction scenarios (A, B, C, D; Fig. 2) C and D were observed most frequently in total, for increases in proportion of both urban (Table 2A, row 8) and agricultural area (Table 2B, row 8). In total, scenarios C and D together comprise 81% of all the positive interaction scenarios as urban area increased, 86% of all the positive interaction scenarios as agricultural area increased. Scenario D, in particular, contributed 63% and 86% of the positive interaction scenarios of species for which the protected area by urban land use interaction or the protected area by agricultural land use interaction was significantly different to zero.

Of the negative interaction scenarios (E, F, G, H; Fig. 2) E and G were observed most frequently in total, for increases in proportion of both the urban (Table 2A, row 8) and agricultural area (Table 2B, row 8). In total, scenarios E and G made up 75% of all negative interaction scenarios as both urban and agricultural area increased. In particular, scenario E contributed 64% and 69% respectively to all the negative interaction scenarios of species for which the protected area by urban land use interaction or protected area by agricultural land use interaction was significantly different to zero (Table A3, panel 2, Appendix).

What is the average change in magnitude of the protection–abundance relationship with increasing proportions of urban and agricultural land in the same pentad?

For the granivores and ground-feeder guilds, the average protection–abundance relationship became significantly more positive with a higher proportion of urban area in the same pentad, whilst this relationship became more negative for frugivores, and neither more positive nor more negative for the remaining four guilds (gleaners, hawkers, predators, and vegivores), Fig. 3A.

As the proportion of agricultural area increased within pentads, the average protection–abundance relationship became significantly more positive for six of the seven guilds: frugivores, vegivores, predators, gleaners, ground feeders and hawkers. For granivores, the average change in the protection–abundance relationship was close to zero, and non-significant (Fig. 3B).

Figure A3 (Appendix) presents a visual summary of the parameter estimates across guilds and Table A2 (Appendix) shows model results for each species within each guild.

The effect of agricultural areas near protected areas on the protection–abundance relationship

The average protection–abundance relationship became more positive in response to an increasing proportion of agricultural area within the same pentad for six of the seven guilds (frugivores, vegivores, predators, gleaners, ground-feeders and hawkers), and none decreased significantly (Fig. 3). We initially hypothesised that the protection–abundance relationship for ground-feeders, hawkers and predators would increase (become more positive) as agricultural and urban areas near protected areas increase in proportion; our results thus concur with this hypothesis. For these guilds, interaction cases C and D (Fig. 2) were the most common ways in which the average protection–abundance relationship increased (Table 2). C and D represent cases in which a negative protection–abundance relationship becomes positive as agricultural area near protected area increases in proportion. These results suggest that in the absence of agricultural land use (and for a scenario where protected areas are neighboured by only natural land-use type), guilds ground-feeders, hawkers, and predators are more abundant outside protected areas, and thus, more abundant within the natural land, compared to inside protected areas. However, as the proportion of agricultural land use surrounding protected areas increases, guilds ground-feeders, hawkers, and predators become more abundant inside protected areas, compared to outside. This result shows a synergistic relationship between protected areas and agricultural land, and on average, these guilds prefer protected areas to agricultural land. This result could possibly be due to the activities associated with farming practices which decrease the quality of the soil, air, and water within entire landscapes, rendering agricultural area less suitable for a broad range of bird guilds (Billeter et al., 2008; Stoate et al., 2001).

The effect of urban areas on the protection–abundance relationship

We predicted the average protection–abundance relationship would increase for the insect eating guilds (ground-feeders, hawkers and gleaners) as the proportion of urban areas near protected areas increased, and decrease for granivores. However, our results only partially supported our hypothesis; on average, guilds ground-feeder and granivores became significantly more positive in their average protection–abundance relationship, whilst this relationship decreased significantly for frugivorous species (Fig. 3A).

The most frequent scenario by which the average protection–abundance relationship increased for the ground-feeder guild was via interaction cases C and D, and via cases D for the granivorous guild. C and D describe a scenario where a negative protection–abundance relationship becomes more positive as urban area increases near protected area. These results suggest that in the absence of urban land near protected areas, ground-feeder and granivorous guilds are more abundant inside protected areas, compared to outside. However, as the proportion of urban area near protected areas increases, ground-feeder and granivorous guilds become more abundant inside protected areas. These results indicate a synergistic relationship between protected areas and urban land for the ground-feeder and granivorous guilds. Our results could possibly be due to the negative effects of dense urban areas on wildlife, such as persecution, predation, or pollution (Blair, 1996), or even the design and use of the urban areas by people (Paker et al., 2014).