Do protection gradients explain patterns in herbivore densities? An example with ungulates in Zambia’s Luangwa Valley

Summary

Ungulate populations face declines across the globe, and populations are commonly conserved by using protected areas. However, assessing the effectiveness of protected areas in conserving ungulate populations has remained difficult. Using herd size data from four years of line transect surveys and distance sampling models, we modeled population densities of four important herbivore species across a gradient of protection on the edge of Zambia’s South Luangwa National Park (SLNP) while accounting for the role of various ecological and anthropogenic variables. Our goal was to test whether protection was responsible for density dynamics in this protection gradient, and whether a hunting moratorium impacted herbivore densities during the studies. For all four species, we estimated lower densities in partially protected buffer areas adjacent to SLNP (ranging from 4.5-fold to 13.2-fold lower) compared to protected parklands. Density trends through the study period were species-specific, with some species increasing, decreasing, or remaining stable in all or some regions of the protection gradient. Surprisingly, when controlling for other covariates, we found that these observed differences were not always detectably related to the level of protection or year. Our findings highlight the importance of accounting for variables beyond strata of interest in evaluating the effectiveness of a protected area. This study highlights the importance of comprehensively modeling ungulate population density across protection gradients, identifies lands within an important protection gradient for targeted conservation and monitoring, documents prey depletion and expands our understanding on the drivers in a critical buffer area in Zambia.

Keywords:

Grasses – Conservation science – Population density – Wildlife – Surface water – Zebras – Lagoons – Herbivory

1 Introduction

In sub-Saharan Africa, ungulates have ecological and economic value through their top-down effects on plant communities and their bottom-up effects as prey for large carnivores [1]. Despite their importance, many protected ungulate populations have recently declined at a rate comparable to populations with less protection [2,3] and face rapid human encroachment [4]. Protected areas (PAs) are an important tool to protect wildlife from human activities and reduce human-wildlife conflict [5]. Strictly protected areas with no permanent settlements and no consumptive use are often bordered by areas with some lower level of protection, creating a gradient to buffer edge effects [6] and source-sink dynamics [7]. However, it remains unclear how effective these protection gradients are in protecting ungulate species of ecological importance or conservation concern [5].

Despite the intuitive benefits of PAs, assessing the effectiveness of PAs in protecting ungulate populations is difficult. First, to provide a valid test of the effect of protection gradients, studies must control for ecological differences between PAs and adjacent buffer zones that could affect ungulate density and distribution. Protected areas are not randomly distributed [8] but instead are generally placed where wildlife densities are high, while buffer zones are often designated in areas with lower wildlife densities. Thus, differences in animal density between PAs and buffer zones can exist due to natural ecological differences between locations unrelated to the effectiveness of their protection. Second, ungulates cannot be surveyed with perfect detection in most habitats [9]. Methods exist to account for imperfect detection when estimating animal densities [10,11], but many studies of protection gradients assume perfect detection or use an index to convert counts that rely on untested assumptions about detection probability. Finally, even when methods that account for detection are used, a common approach is to model the density of groups with distance sampling [10] and then convert group density to individual density using a mean group size [10,12]. This conversion uses either mean group size across all observations or across focal categorical variables, such as vegetation types [13]. Because ungulate group size is typically influenced by the same variables that affect the distribution of herds, this approach may not be entirely accurate. Important broad studies have aided our understanding of the effectiveness of PAs to conserve ungulates [14–17], but their inferences have been constrained by one or more of these limitations.

Zambia contains several PAs important for regional conservation of ungulates and the large carnivores that depend on ungulate prey [18], and most of these PAs face human encroachment that is approaching rapidly or already has reached the PA itself [19]. Most Zambian PAs are buffered by Game Management Areas (GMAs) that allow some human settlement and support consumptive uses of wildlife and resources, such as legal trophy hunting [20]. Rapid human population growth in GMAs has brought increased pressures of illegal bushmeat harvest and habitat conversion, challenging the effectiveness of Zambia’s GMAs as buffers for adjacent PAs [21]. South Luangwa National Park (SLNP) is Zambia’s premier PA for wildlife-based photo-tourism and conserves regionally-important populations of several threatened and endangered species yet faces rapid human pressure from adjacent GMAs [19]. Recent studies of large carnivore demography and dynamics [20,22–25] and studies of bushmeat poaching patterns [21,26,27] suggest that the depletion of ungulate populations from bushmeat poaching may be affecting the ecological integrity of this protection gradient. Aerial surveys also suggest overall ungulate declines in the Luangwa Valley and lower ungulate density in GMAs compared to PAs [21], but these studies are generally inconclusive due to low precision of density estimates and no correction for variability in detection. Finally, it has been suggested that a temporary moratorium in all trophy hunting from 2013 to 2014 allowed increased poaching activity in the absence of the primary wildlife-based tourism activity in GMAs [28]. Despite this concern and reported increased poaching in communities in GMAs adjacent to SLNP [29], there has not been sufficient data to test the relationship between the hunting moratorium and herbivore population trends. Therefore, to better inform management actions in this important PA, there is a clear need for unbiased and precise estimates of ungulate abundance in SLNP and adjacent GMAs, for tests of the relationship between ungulate density and the level of protection, and for tests of trends in density over time and across management actions, while addressing the difficulties introduced above.

Here we use data from repeated, stratified line transect sampling to fit distance sampling and group size models to estimate population densities of impala (Aepyceros melampus), puku (Kobus vardonii), plains zebra (Equus quagga), and warthog (Phacochoerus africanus) across the South Luangwa Protection Gradient from 2012–2015. These species are abundant in the area, important prey for carnivore species of concern, primary targets for illegal bushmeat trade [30] and are expected to be negatively impacted by the absence of trophy hunting during the 2013–2014 moratorium. For each species we model both group density and group size as functions of top-down, bottom-up, abiotic, and anthropogenic covariates, and then estimate differences in population density across space and time while controlling for these effects. Our approach addresses the above difficulties in assessing the efficacy of protected areas as it accounts for imperfect detection, robustly integrates ecological and anthropogenic covariates that impact both the distribution and size of ungulate groups, compares ungulate densities within a PA and GMA that are similar ecologically but differ in human usage, and assess population trends with the temporary cessation of trophy hunting activities. The results improve our understanding of the efficacy of protection gradients in buffering the impacts of human encroachment and bushmeat poaching and provide baseline estimates of density for ecologically important ungulate populations that face increasing anthropogenic pressures.

2. Materials and methods

2.1 Study area

Our study area (hereafter the South Luangwa Protection Gradient, or SLPG) was within a 1 200 km² complex of grasslands, scrublands, and forests situated along the perennial Luangwa River, the primary eastern boundary of SLNP and the adjacent Lupande GMA (hereafter the GMA; Fig 1) [31,32]. The Luangwa River is only a seasonal barrier for wildlife and attracts the highest densities of wildlife in the area as it is the largest and most reliable perennial water source in the region, particularly during the dry season (May-November).

The SLPG includes three regions of interest for this study. West of the Luangwa River are fully protected SLNP lands (818 km²), where the primary human activity is photographic safari tourism and Park management activities. East of the Luangwa River there is a section of SLNP (151 km²) that is also popular for photographic tourism but is thought to be exposed to more illegal bushmeat poaching due to open borders and a public road that bisects the area. While this is a small, unsurfaced dirt track, it supports considerable local foot and bicycle traffic, and some vehicular traffic. Finally, the GMA lands (231 km²), also east of the Luangwa River, contain a growing human population (annual growth of 3.8%) [33] that has raised an array of conservation concerns (see Section 1). Trophy hunting operates in leased, unfenced concessions in the GMA (save for a moratorium from 2013–2014) [20], targeting an array of wildlife species including those studied here [31]. In short, our study area encompasses areas of high conservation importance across a protection gradient with associated variation in human influence and other potentially limiting factors.

2.2 Study design

2.2.1 Survey design

We used line transects and distance sampling [10] to estimate herd density (herds/km²) and herd size (individuals/herd) for ungulates across the SLPG. We surveyed 15 transects by vehicle across our study area (120.5 km in total), with sufficient spacing between transects to minimize the risk of double-counting groups [10] (Fig 1). We established most transects in a generally east-west orientation (i.e., perpendicular to the Luangwa River) to sample across the range of environmental and anthropogenic variables of interest. Topography and vegetation impassible by vehicle forced us to follow roads for sections of transects, which we accounted for in our modeling of detection (See Section 2.2.2). Each transect was split into segments defined by observed changes in vegetation cover type or a maximum length of 2km, resulting in a total of 97 segments across the three regions of interest in the SLPG (West of the Luangwa and National Park, n = 62; East of the Luangwa and National Park, n = 15; GMA, n = 20). These segments allowed us to estimate the role of covariates on a finer scale than for entire transects (e.g. distance to the Luangwa River; See Section 2.2.2). We were not able to achieve balanced sampling across the SLPG due to human settlement and agriculture activity in the GMA. GMA transects were primarily on primitive roads or off-road, to minimize any sensitivity of animals to vehicle noise. As illegal harvests of wildlife are usually done on foot, these species did not appear to be sensitive to our driven transects. We recorded herbivore group observations and covariate values on the segment level as required by our analysis described below.

We conducted 10 surveys of all 97 segments during daylight hours (0900–1700) at the beginning (June), middle (August), and end (October) of the dry season (May-November) from 2012–2015 (n = 970 segment-surveys). June 2012 and August 2015 were not surveyed due to logistic constraints. We could not conduct any surveys during the rainy season (December-April) as much of the study area is inaccessible. During each survey, the vehicle was driven at a maximum speed of 15 km/h, with two observers seated on the rooftop scanning for animal herds. When a herd or single animal was seen, the observers recorded the species composition and size of the herd, and the distance (aided by a laser rangefinder) and azimuth to the herd following standard distance sampling protocol [10]. We reduced our risk of double-counting animal herds during each survey by surveying all segments over a period of 3–4 days.

2.2.2 Factors affecting group density and size

We considered a suite of bottom-up, top-down, anthropogenic, and abiotic variables that could affect herd density and herd size of the focal species on the segment level (Table 1). We measured ten covariates that characterize bottom-up effects of vegetation composition and structure on herd density and herd size. While surveying segments we recorded the presence or absence of green grass and of grassy lagoon patches, and the predominant grass height category (short = <10cm, intermediate: 10-100cm, long: >100cm). We suspected that vegetation composition would influence both herd density and herd size, and our ability to detect herds. To model detection, we assigned each segment to one of six vegetation structure categories, as well as one of three simplified vegetation categories (Table 1): because transects were split into segments at irregular intervals that were defined by changes in vegetation structure, each segment was relatively homogenous in structure. We considered these two structure covariates separately during our analyses as they were two descriptions of the same underlying variable. To model group density and herd size, we classified vegetation composition and heterogeneity within a 240m-wide buffer around each segment (hereafter buffer) to characterize vegetation beyond the view of the survey team. Using Landsat 7 Global Land Survey 2010 imagery we partitioned pixels into 10 vegetation classes of differing vegetation densities using the clara() function in the cluster package [34] in R [35]. We then combined vegetation classes into four cover types based on vegetation density and validated these cover type classifications during our line transect surveys. These cover types included closed scrubland (right-skewed, and therefore log-transformed), closed woodland, open woodland, and open grassland. We estimated the proportion of buffer composed of each cover type for each segment using the sp package in R [36,37] and estimated the density of edges between cover types within each buffer (km edge/km²) using the perimeter tool in QGIS 2.0.1 (www.qgis.org).

Differences in predation risk could affect herd density or herd size across the SLPG, thus we quantified predation risk across our study area by measuring the utilization of the area by the lion (Panthera leo) population (the most abundant large carnivore within the study area) [20]. Intensive lion population monitoring was ongoing during this study, with prides and coalitions monitored across the 1200 km study area detailed in this manuscript. We fit a single kernel utilization distribution to 7 785 lion locations collected over five years (2008–2012) from lions equipped with GPS radiocollars and direct observation of 18 lion prides and 14 male coalitions [20]. Though this spatial information overlaps only one year of this study, ongoing lion studies indicate little change in lion distribution [25]. We used sp, rgdal [38], and plyr [39] packages in R to estimate the distribution of daily distance moved from six lionesses equipped with GPS radiocollars (range: 0m - 17 300m) and used the 90th percentile (6 974m) as a smoothing parameter for the utilization distribution (UD) [40]. We fit the UD using the bivariate normal kernel function in the adehabitatHR package [41], with an output grid of 300m x 300m. We converted the UD to a raster using the raster [42] and sp packages and extracted and standardized lion UD values for each segment’s midpoint, thereby quantifying the risk of herbivores encountering the dominant large carnivore species in each segment.

The impact of anthropogenic activities was characterized for each segment by classifying whether it lay within the PA or partially-protected GMA (hereafter protection), on the east or west side of the Luangwa River (hereafter side), and the distance from the segment midpoint to the nearest road (also right-skewed, and log-transformed). We could not test for an interaction between protection and side as there are no GMA areas within the study area on the western side of the Luangwa River. We also could not include distance to park boundary in our analysis as this boundary is mostly defined by the Luangwa River. Despite the Luangwa River defining a boundary in the SLPG, any effect of this boundary in our study is confounded by the abiotic influence of the Luangwa River as the SLPG’s perennial water source. Therefore, we chose to treat the Luangwa River as an abiotic variable, and not an anthropogenic variable (see below). Finally, as some segments followed roads, we classified each segment’s path type to account for any biases in detection for road-based segments (Table 1) [43].

We measured six abiotic covariates to characterize the environmental conditions potentially influencing herbivore density and distribution in the Luangwa Valley (Table 1). Availability of water is a limiting factor for wildlife in the Luangwa Valley, with water sources diminishing and disappearing with the progression of the dry season. Therefore, we measured distances from each segment’s midpoint to the perennial Luangwa River and the nearest seasonal stream and recorded the dry season stage and whether standing water was present during surveys of each segment. We also recorded the year of each survey to determine annual trends in group density and size. Finally, we recorded whether there was evidence of a fire that had burned through the segment for each survey, as post-fire “green flush” provides herbivores with high-quality forage and potentially reduces predation risk [44].

2.3 Analytical methods

2.3.1 Herd density analysis

We used a multinomial generalized distance sampling model to estimate the ‘super-population’ of herds (λ) for each species within each segment, while accounting for imperfect detection (p) and varying availability for detection at the time of the survey (φ)[45]. We used the gdistsamp() function to fit candidate models for each species in the unmarked package [45] in R. We truncated herd sighting data for each species to exclude outlier distances that would compromise estimation of detection probability.

To focus on accurate estimation of the parameter of interest (group density), we evaluated a set of models in three steps using Akaike’s Information Criteria (AIC; Table 2). We first identified the best-supported model for detection across all combinations of covariates thought to influence detection using hazard, half-normal, and uniform detection functions, while estimating only a mean for availability and super-population (Tables 1 and 2). Using the best-supported detection model and a uniform super-population model, we next identified the best supported model for availability across all combinations of covariates thought to influence availability (Tables 1 and 2). In the third step, while modeling detection and availability using their best-supported models, we used AIC model selection to identify the best supported models (≤ 1 AIC) for each of the four types of effect on the super-population of herds: bottom-up, top-down, anthropogenic, and abiotic (Table 2). Finally, we took the top model for each of these four types and used AIC scores to identify the final density models (≤ 1 AIC), selecting from a set with all additive combinations of effects in the top models for each type (Table 2). This multi-stage approach to limit the number of models compared was necessary to reduce computation time. We used deviance goodness-of-fit tests for each species’ top herd density model(s) to examine model fit. We used these final models and model averaging (when more than one model received comparable support from the data) to estimate herd super-population (λ) and availability (φ) using the predict() function, and derive herd density estimates (D^) as the product of these two parameter estimates.

2.3.2 Herd size analysis

We fit zero-truncated Poisson (ZTP) regression models using the vglm() function in the VGAM package [46,47] in R to estimate segment-specific mean herd sizes for the four focal species while controlling for the effects of covariates described above (Table 1). First, we used a ZTP model to test whether herd size was affected by distance from the transect, as is expected if herd size affects detection, and only used herd size observations at 0m from the transect if this there was evidence of this effect (p < 0.15) [10]. Next, we dropped covariates with badly imbalanced sampling or that were highly correlated with other covariates (Pearson’s r > 0.6) from consideration. We then used reverse step-wise likelihood-ratio (LR) tests to select a herd size model for each species, and then confirmed this model selection using forward step-wise LR tests. We evaluated model fit by fitting a linear regression of Pearson’s residuals on predicted herd sizes and assumed adequate fit if the estimated intercept and slope were not detectably different from 0 [46]. Using the final ZTP model and the predict() function we estimated differences in herd size across covariate ranges and variation in herd size while accounting for the non-random distribution and correlation of covariates across the study area. To avoid extrapolation, we only estimated mean herd size for segments with covariate values within the range documented during herd observations.

2.3.3 Population density analysis

To estimates population densities for each species, we multiplied mean herd density by estimates of mean herd size for each segment and used non-parametric bootstrapping to estimate mean population density and its variance across SLPG regions and years. We also predicted herd density and herd size with all covariates other than protection, side, and year fixed at their mean value, estimated mean population density and variance across those strata of interest, and calculated differences across SLPG regions and across years. In summary, our comprehensive approach allows a test for the effect of protection and year on population density that does not ignore the ecological and abiotic differences across the protection gradient.

3. Results

3.1 Herd density and herd size

During the 10 surveys of 97 segments, we detected 890 impala herds, 478 puku herds, 175 zebra herds, and 169 warthog herds. After truncating datasets to maintain suitable detection probabilities (400m for puku, 300m for the other three species), our final sample sizes for our herd density analysis were 836 impala herds, 452 puku herds, 155 zebra herds, and 163 warthog herds (S1 Fig). Distance had a positive association with impala and puku herd size (p<0.0001 and p = 0.005, respectively), but not with zebra and warthog herd size (p = 0.43 and p = 1, respectively). These results indicated that large impala and puku herds were more detectable at large distances, so to avoid bias we only used herds that were directly on the transect to estimate herd size for these two species. Our final sample for our herd size analysis included 122 impala herds (mean = 7.1 individuals, range: 1–75), 56 puku herds (mean = 6.6 individuals, range: 1–111), 155 zebra herds (mean = 5.1 individuals, range: 1–24), and 163 warthog herds (mean = 2.6 individuals, range: 1–9).

The best supported herd density model(s) for each parameter varied between species (Table 3), and coefficient estimates varied in magnitude and sign between models (Table 4). Despite this variation two variables of primary interest—side of Luangwa River and protection status—were included in top herd density models for all species. Herd density was higher in the National Park for all species. Segments west of the Luangwa River were estimated to have higher impala and zebra herd densities, whereas puku and warthog densities were lower or not detectibly different than in segments east of the Luangwa River, respectively. Survey year influenced availability for zebra and warthog herd density, with herd density detectably higher for zebra in 2015, but detectably lower for warthog in 2015.

Like herd density models, herd size models and coefficient estimates varied between species (Table 5). Segments on the western side of the Luangwa River were associated with smaller herds for impala (21% smaller) and puku (16% smaller), while neither side of the Luangwa River nor protection status were detectably correlated with zebra or warthog herd sizes. Impala and puku herds were smallest in 2012, with detectably larger herds during the remainder of the study period (except 2014 for puku).

3.2 Population density

Impala had the highest population density, with an average density of 27.69 animals/km², and showed great spatiotemporal variation in density (range among segments: 0.24–601.18 animals/km²). Puku also occurred at variable densities, an average density of 8.11 animals/km² (range: 0.003–172.60 animals/km²). Zebra and warthog occurred at lower and less variable densities, with average densities of 2.41 animals/km² (range: 0.03–28.09 animals/km²) and 1.76 animals/km² (range: 0.05–13.82 animals/km²), respectively.

3.2.1 Protection effects

Herbivore densities varied widely across the SLPG, and the highest densities for all species were within the fully-protected SLNP (Fig 2). Within SLNP, density varied between portions of the park that were west and east of Luangwa River. Impala and zebra densities did not detectably differ within these two areas of the national park, while both puku and warthog occurred at their highest densities in eastern parklands. Thus, there were no consistent differences between areas east and west of the river that had the same legal protection. All species occurred at lower densities in the GMA, ranging from 4.5-fold to 13.2-fold lower than densities within SLNP, with no overlap in 95% CIs.

3.2.2 Temporal effects

Density estimates varied among years, but annual differences from regional mean densities were species-specific, with no consistent pattern for across species (Fig 2). Impala density did not detectably change across the study period in any section of the protection gradient. Puku density increased from 2012 to 2013 in western parklands but did not differ from 2014 and 2015 estimates or the overall mean density for that region. Zebra densities did not change in GMA and eastern parklands but increased in western parklands. Finally, we detected a decline in warthog densities across all regions in the protection gradient.

3.2.3 Predicted population density across the SLPG and the study period

After controlling for all covariates, we found that the observed differences in population density just described were not always detectably related to the level of protection (Fig 3) or year (Fig 4). Across the SLPG, predicted mean species densities were consistently higher in protected areas, but some of these estimated differences had low precision (Fig 3). For impala, populations in protected parklands east and west of the Luangwa River were estimated to have 1.3- and 0.66-fold greater densities compared to GMA densities, but all CIs showed that predicted mean density overlapped with GMA estimates. Puku densities were 9.5-fold greater in the eastern parklands (with >80% confidence that a difference existed due to protection), while densities in western parklands were not detectably different from the GMA. Predicted zebra densities were 5-fold and 12-fold greater in the eastern parklands and western parklands, respectively, with high confidence that a protection effect existed for both comparisons. Predicted warthog densities were higher in protected regions compared to the GMA, but CIs overlapped predicted GMA densities. Throughout the study period we could not detect any density trends for any species after controlling for other covariates (Fig 4), including during the trophy hunting moratorium (2013–2014).

In summary: (1) Ungulate densities were consistently higher in the better protected SLNP. (2) There was evidence that higher protection status yields higher density, after controlling for ecological and abiotic differences between NP and GMA, but this evidence was mixed and sometimes not strong. (3) After controlling for other effects, there was no consistent pattern of change in ungulate densities over the four years of observation.

4. Discussion

Our study of the SLPG supports previous findings that lands with better protection generally hold higher densities of herbivores [14–17], supporting continued protection and monitoring in the SLPG. In particular, the segment of SLNP east of the Luangwa River still supports densities of ungulates comparable to or greater than in SLNP west of the Luangwa River, despite having no physical boundary to prevent human incursions and a heavily-used public road passing through it. We estimated the lowest densities for all focal species in the SLPG within the GMA, which supports conclusions from previous studies in the SLPG that suggested depleted herbivore populations within Lupande GMA [23,26]. Given that our GMA transects surveyed areas within 6 km of the Luangwa River, where wildlife densities are thought to be highest and away from dense human settlements, the densities we estimated within this buffer area are likely high compared to the rest of the GMA. These inferences are restricted to dry season conditions, as we were not able to survey transects in wet season conditions.

However, after controlling for a suite of bottom-up, top-down, anthropogenic, and abiotic covariates, there is no clear evidence that these dynamics were driven by protection status or year alone. The non-anthropogenic variables included in our analysis incorporate ecosystem alteration by humans, and thus our interpretation of the role of the SLPG regions and time are proxies for anthropogenic mortality and associated risk effects [48]. Isolating these forces from ecological variables is not entirely clear; for example, poaching efforts are non-random and are correlated with both anthropogenic and ecological variables [27]. We did not allow trends to vary across the SLPG, so any differences in trend between regions for a species are due to covariates other than year (Fig 2), supported by lacking trends when all covariates other than year were held constant. With these considerations, our findings indicate that an array of ecological and anthropogenic variables influence herbivores in protected area networks characterized by national parks (or other strictly protected areas) and adjoining buffer zones.

If we had not applied our rigorous modeling of both group density and size, and instead followed previous approaches (see Section 1), our inferences would have been confined to the variables of interest (protection and year) and would not acknowledge the effects of other important variables. With our approach, ecological variables could play an important role in our modeling process, highlighting the importance of considering ecological conditions in the SLPG along with the protections provided by SLNP. Herbivore density can vary across gradients of protection, but our findings indicate that the protection status of an area is generally insufficient to capture the complexities of factors influencing herbivores across this protection gradient.

4.1 Differences in density across the South Luangwa Protection Gradient

Illegal wire snare poaching has been a well-recognized threat to wildlife communities in the Luangwa Valley [30] and has been identified as a major threat to the persistence of large carnivores [26]. In this study, GMA transects represented areas that were prone to high risk for snare occurrence relative to risk in SLNP due to proximity to human activity and greater law enforcement efforts within SLNP [27]. Therefore, we would expect reduced herbivore densities in GMA relative to SLNP regions if bushmeat poaching was impacting herbivore densities. After controlling for all other variables, our density estimates were higher in PA regions of the SLPG, but only zebra densities reflected this expectation with 95% confidence (puku reflected this expectation with 80% confidence; Fig 3). The detected role of the SLPG in zebra density dynamics indicates evidence of increased impacts of bushmeat poaching and other human activities, particularly as larger ungulate species are vulnerable to overharvest [3]. While we cannot identify bushmeat poaching as the sole driver behind the role of the SLPG in zebra density dynamics (survival data are unavailable), other sources of mortality are unlikely drivers in this study. Trophy hunting occurred with limited quotas and small harvest rates relative to our population density estimates, and lion utilization was accounted for in our modeling of group density and size (Table 4). Therefore, the predictable impact of wire-snare poaching on zebra (and likely puku) is a signal for concern but demonstrates that the protection gradient is providing protection from human exploitation.

Despite concerns of increased poaching in the absence of wildlife-based economies during the trophy hunting moratorium in GMAs, there is no evidence of any coinciding large-scale herbivore decline in our study area. Anti-poaching efforts in the study area were conducted jointly by the Department of National Parks and Wildlife and Conservation South Luangwa during this period, and the level of anti-poaching effort and support increased during the moratorium [29,49–51]. We acknowledge that this is not necessarily the case in all concessions but does indicate the potential for greater anti-poaching investment by operators in these areas. We did not test for interactions between year and region of the protection gradient, so our estimates of population trends could hide local declines in the GMA. However, these results indicate no evidence that the 2013–2014 trophy moratorium was detrimental to the species in this study.

Estimated densities across the SLPG and throughout the study period are clearly influenced by ecological variables, indicating that habitat alteration likely plays an important role in decreasing available resources and reducing herbivore densities in the GMA. Colleagues [19] documented rapid rates of human encroachment in Lupande and other GMAs around key Zambian protected areas. Our GMA transects did not pass directly through any majorly altered areas (e,g, agriculture or villages), but such areas were present across the GMA region of the SLPG. In addition to human encroachment, altered composition of the wildlife community in the GMA, particularly the poaching of elephants and rhinos (extirpated by 1995) [52], likely influences vegetation structure in this region [53]. The role of habitat conversion and changing wildlife community composition should be further investigated in the SLPG, alongside ongoing attention to combat illegal bushmeat poaching.

4.2 The future role of ground transects in the South Luangwa Protection Gradient

Future monitoring is critical to track population densities for these and other herbivore species, particularly in the face of a growing human population. While multiple approaches are implemented to achieve this monitoring, we advocate continued monitoring by stratified, ground-based transects. Aerial survey data collected by the Zambia Department of National Parks and Wildlife supports the notion that the SLPG supports notably high densities for all four species [54], and therefore remains an important area to protect and monitor. While ground-based distance sampling surveys cannot match the spatial scale achievable by aircraft, ground-based surveys provide more accurate density estimates to identify population trends and aid in our understanding of the function of protected areas. The difference in costs between the two methods is well-documented [55], and population monitoring could be supplemented with data collected by law enforcement patrols [56] and citizen science initiatives [57] to offset costs and to achieve a large study area. Ground-based distance sampling methods should be integrated into long-term monitoring of the SLPG and in other protected areas, focused on critical areas for ungulate populations and accounting for the dynamics of ecological covariates in both group density and size, as demonstrated in this study.

Supporting information

S1 Appendix [xlsx]
Detection and covariate data from line transect surveys.

S1 Fig [tif]
Histogram of ungulate group detections.