Do protection gradients explain patterns in herbivore densities? An example with ungulates in Zambia’s Luangwa Valley
Elias Rosenblatt aff001; Scott Creel aff001; Paul Schuette aff001; Matthew S. Becker aff001; David Christianson aff001; Egil Dröge aff005; Thandiwe Mweetwa aff001; Henry Mwape aff001; Johnathan Merkle aff001; Jassiel M’soka aff001; Jones Masonde aff006; Twakundine Simpamba aff006
Authors place of work:
Zambian Carnivore Programme, Mfuwe, Eastern Province, Zambia
aff001; Department of Ecology, Montana State University, Bozeman, Montana, United States of America
aff002; Alaska Center for Conservation Science, University of Alaska Anchorage, Anchorage, Alaska, United States of America
aff003; School of Natural Resources and the Environment, University of Arizona, Tucson, Arizona, United States of America
aff004; Wildlife Conservation Research Unit, University of Oxford, Tubney, United Kingdom
aff005; Department of National Parks and Wildlife, Chilanga, Zambia
Published in the journal:
PLoS ONE 14(10)
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.
Conservation science – Grasses – Herbivory – Lagoons – Population density – Surface water – Wildlife – Zebras
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 . 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 . Protected areas (PAs) are an important tool to protect wildlife from human activities and reduce human-wildlife conflict . 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  and source-sink dynamics . However, it remains unclear how effective these protection gradients are in protecting ungulate species of ecological importance or conservation concern .
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  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 . 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  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 . 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 , and most of these PAs face human encroachment that is approaching rapidly or already has reached the PA itself . 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 . 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 . 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 . 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 , 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 . Despite this concern and reported increased poaching in communities in GMAs adjacent to SLNP , 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  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 km2 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 km2), 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 km2) 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 km2), also east of the Luangwa River, contain a growing human population (annual growth of 3.8%)  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) , targeting an array of wildlife species including those studied here . 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  to estimate herd density (herds/km2) 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  (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 . 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  in R . 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/km2) 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) . 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 . Though this spatial information overlaps only one year of this study, ongoing lion studies indicate little change in lion distribution . We used sp, rgdal , and plyr  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) . We fit the UD using the bivariate normal kernel function in the adehabitatHR package , with an output grid of 300m x 300m. We converted the UD to a raster using the raster  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) .
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 .
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 (φ). We used the gdistsamp() function to fit candidate models for each species in the unmarked package  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) . 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 . 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.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/km2, and showed great spatiotemporal variation in density (range among segments: 0.24–601.18 animals/km2). Puku also occurred at variable densities, an average density of 8.11 animals/km2 (range: 0.003–172.60 animals/km2). Zebra and warthog occurred at lower and less variable densities, with average densities of 2.41 animals/km2 (range: 0.03–28.09 animals/km2) and 1.76 animals/km2 (range: 0.05–13.82 animals/km2), 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.
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 . 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 . 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  and has been identified as a major threat to the persistence of large carnivores . 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 . 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 . 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  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) , likely influences vegetation structure in this region . 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 , 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 , and population monitoring could be supplemented with data collected by law enforcement patrols  and citizen science initiatives  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.
1. Gordon IJ, Hester AJ, Festa‐Bianchet M. Review: The management of wild large herbivores to meet economic, conservation and environmental objectives. J Appl Ecol 2004;41: 1021–1031.
2. Western D, Russell S, Cuthill I. The status of wildlife in protected areas compared to non-protected areas of Kenya. PloS one 2009; 4:e6140. doi: 10.1371/journal.pone.0006140 19584912
3. Ripple WJ, Newsome TM, Wolf C, Dirzo R, Everatt KT, Galetti M, et al. Collapse of the world’s largest herbivores. Sci Adv 2015;1:e1400103. doi: 10.1126/sciadv.1400103 26601172
4. Wittemyer G, Elsen P, Bean WT, Burton ACO, Brashares JS. Accelerated human growth at protected area edges. Science 2008;321: 123–126. doi: 10.1126/science.1158900 18599788
5. Geldmann J, Barnes M, Coad L, Craigie ID, Hockings M, Burgess ND. Effectiveness of terrestrial protected areas in reducing habitat loss and population declines. Biol Conserv 2013;161: 230–238.
6. Woodroffe R, Ginsberg J. Edge effects and the extinction of populations inside protected areas. Science 1998;280: 2126–2128. doi: 10.1126/science.280.5372.2126 9641920
7. Pulliam HR. Sources, sinks, and population regulation. Am Nat 1988;132: 652–661.
8. Joppa L, Pfaff A. Reassessing the forest impacts of protection: the challenge of nonrandom location and a corrective method. Ann. N.Y. Acad. Sci. 2010;1185: 135–149. doi: 10.1111/j.1749-6632.2009.05162.x 20146766
9. M'soka J, Creel S, Becker MS, Murdoch JD. Ecological and anthropogenic effects on the density of migratory and resident ungulates in a human-inhabited protected area. Afr J Ecol 2017;55: 618–631.
10. Buckland ST. Introduction to Distance Sampling: Estimating Abundance of Biological Populations. Oxford: Oxford University Press; 2001.
11. Royle JA, Dawson DK, Bates S. Modeling abundance effects in distance sampling. Ecol 2004;85: 1591–1597.
12. Ferguson MC, Barlow J, Fiedler P, Reilly SB, Gerrodette T. Spatial models of delphinid (family Delphinidae) encounter rate and group size in the eastern tropical Pacific Ocean. Ecol Model 2006;193: 645–662.
13. Rduch V. Population characteristics and coexistence of puku (Kobus vardonii) and impala (Aepyceros melampus) in and around Kafue National Park, Zambia. Mamm Biol 2016;81: 350–360.
14. Caro TM, Pelkey N, Borner M, Campbell KLI, Woodworth BL, Farm BP, et al. Consequences of different forms of conservation for large mammals in Tanzania: preliminary analyses. Afr J Ecol 1998;36: 303–320.
15. Stoner C, Caro T, Mduma S, Mlingwa C, Sabuni G, Borner M. Assessment of effectiveness of protection strategies in Tanzania based on a decade of survey data for large herbivores. Conserv Biol 2007;21: 635–646. doi: 10.1111/j.1523-1739.2007.00705.x 17531042
16. Ogutu JO, Piepho HP, Dublin HT, Bhola N, Reid RS. Dynamics of Mara–Serengeti ungulates in relation to land use changes. J Zool 2009;278: 1–14.
17. Schuette P, Creel S, Christianson D. Ungulate distributions in a rangeland with competitors, predators and pastoralists. J Appl Ecol 2016;53: 1066–1077.
18. Purchase G, Mateke C, Purchase D. A review of the status and distribution of carnivores, and levels of human- carnivore conflict, in the protected areas and surrounds of the Zambezi Basin. Unpublished report. The Zambezi Society, Bulawayo. 2007. pp 1–77.
19. Watson FGR, Becker MS, Milanzi J, Nyirenda M. Human encroachment into protected area networks in Zambia: implications for large carnivore conservation. Reg Environ Change 2014;15: 415–429.
20. Rosenblatt E, Becker MS, Creel S, Droge E, Mweetwa T, Schuette PA et al. Detecting declines of apex carnivores and evaluating their causes: An example with Zambian lions. Biol Conserv 2014;180: 176–186.
21. Lindsey PA, Nyirenda VR, Barnes JI, Becker MS, Mcrobb R, Tambling CJ, et al. Underperformance of African protected area networks and the case for new conservation models: Insights from Zambia. PLoS One 2014;9: e94109 doi: 10.1371/journal.pone.0094109 24847712
22. Becker MS, Watson FGR, Droge E, Leigh K, Carlson RS, Carlson AA. Estimating past and future male loss in three Zambian lion populations. J Wildl Manag 2013;77: 128–142.
23. Rosenblatt E, Creel S, Becker MS, Merkle J, Mwape H, Schuette P, Simpamba T. Effects of a protection gradient on carnivore density and survival: an example with leopards in the Luangwa valley, Zambia. Ecol Evol 2016;6: 3772–3785. doi: 10.1002/ece3.2155 27231529
24. Creel S, M'soka J, Dröge E, Rosenblatt E, Becker M, Matandiko W, Simpamba T. Assessing the sustainability of African lion trophy hunting, with recommendations for policy. Ecol Appl 2016;26: 2347–2357. doi: 10.1002/eap.1377 27755732
25. Mweetwa T, Christianson D, Becker MS, Creel S, Rosenblatt E, Merkle J, et al. Quantifying lion demographic responses during a three-year moratorium on trophy hunting. PLOS One 2018;13: e0197030. doi: 10.1371/journal.pone.0197030 29782514
26. Becker MS, McRobb R, Watson F, Droge E, Kanyembo B, Murdoch J, Kakumbi C. Evaluating wire-snare poaching trends and the impacts of by-catch on elephants and large carnivores. Biol Conserv 2013;158: 26–36.
27. Watson F, Becker MS, McRobb R, Kanyembo B. Spatial patterns of wire-snare poaching: implications for community conservation in buffer zones around National Parks. Biol Conserv 2013;168: 1–9.
28. White PA, Belant JL. Provisioning of game meat to rural communities as a benefit of sport hunting in Zambia. PloS one 2015;10: e0117237. doi: 10.1371/journal.pone.0117237 25693191
29. SLCS. South Luangwa Conservation Society Annual Report 2014. 2015. 22 pp.
30. Lewis DM, Phiri A. Wildlife snaring—an indicator of community response to a community-based conservation project. Oryx 1998;32: 111–121.
31. Astle WL, Webster R, Lawrance CJ. Land classification for management planning in the Luangwa Valley of Zambia. J Appl Ecol 1969;6: 143–169.
32. White F. The Vegetation of Africa: A Descriptive Memoir to Accompany the Unesco/AETFAT/UNSO Vegetation Map of Africa. Paris, France: United Nations Educational, Scientific and Cultural Organization; 1983.
33. Republic of Zambia. 2010 census of population and housing. Lusaka, Zambia; Central Statistical Office; 2011.
34. Maechler M, Rousseeuw P, Struyf A, Hubert M, Hornik K. cluster: Cluster analysis basics and extensions. R package version 2.0.4; 2016.
35. R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2019.
36. Pebesma EJ, Bivand RS. Classes and methods for spatial data in R. R News 5; 2005.
37. Bivand RS, Pebesma E, Gomez-Rubio V. Applied spatial data analysis with R. 2nd ed. New York: Springer; 2013.
38. Bivand RS, Keitt T, Rowlingson B. rgdal: Bindings for the geospatial data abstraction library. R package version 1.0–4; 2015.
39. Wickham H. The split-apply-combine strategy for data analysis. J Stat Softw 2011;40: 1–29.
40. Drӧge E, Creel S, Becker M, M'soka J. Spatial and temporal avoidance of risk within a large carnivore guild. Ecol Evol 2017;7: 189–199. doi: 10.1002/ece3.2616 28070283
41. Calenge C. The package adehabitat for the R software: a tool for the analysis of space and habitat use by animals. Ecol Model 2006;197: 516–519.
42. Hijmans RJ. raster: Geographic data analysis and modeling. R package version 2.4–15; (2015).
43. Anderson DR, Laake JL, Crain BR, Burnham KP. Guidelines for line transect sampling of biological populations. J Wildl Manag 1979;43: 70–78.
44. Wilsey BJ. Variation in use of green flushes following burns among African ungulate species: the importance of body size. Afr J Ecol 1996;34: 32–38.
45. Chandler RB, Royle JA, King DI. Inference about density and temporary emigration in unmarked populations. Ecol 2011;92: 1429–1435.
46. Zuur AF, Leno EN, Walker NJ, Savaliev AA, Smith GM. Mixed Effects Models and Extensions in Ecology with R. New York: Springer; 2009.
47. Yee TW. Vector Generalized Linear and Additive Models: With an Implementation in R. New York: Springer; 2015.
48. Frid A, Dill L. Human-caused disturbance stimuli as a form of predation risk. Conserv Ecol 2002;6: 11.
49. SLCS. Annual Report to Luangwa Conservation Community Fund. 2012. 22 pp.
50. SLCS. South Luangwa Conservation Society Annual Report 2013. 2014. 24 pp.
51. CSL. Conservation South Luangwa Annual Report 2015. 2015. 21 pp.
52. Nyirenda VR, Lindsey PA, Phiri E, Stevenson I, Chomba C, Namukonde N, et al. Trends in illegal killing of African elephants (Loxodonta africana) in the Luangwa and Zambezi ecosystems of Zambia. Environ Nat Resour Res 2015;5: 24–36.
53. Sankaran M, Augustine DJ, Ratnam J. Native ungulates of diverse body sizes collectively regulate long‐term woody plant demography and structure of a semi‐arid savanna. J Ecol 2013;101: 1389–1399.
54. DNPW. Report on the 2015 aerial census of elephants and other large mammals in Zambia: Volume II Population estimates for other large mammals and birds. Chilanga, Zambia: Department of National Parks and Wildlife; 2016.
55. Jachmann H. Estimating Abundance of African Wildlife: An Aid to Adaptive Management. Boston: Kluwer Academic Publishers; 2001.
56. Schuette P, Namukonde N, Becker MS, Watson F, Creel S, Chifunte C, et al. Boots on the ground: in defense of low‑tech, inexpensive, and robust survey methods for Africa’s under‑funded protected areas. Biodivers Conserv 2018;27: 2173–2191.
57. Marnewick K. Conservation biology of cheetahs Acinonyx jubatus (Schreber 1775) and African wild dogs Lycaon pictus (Temminck 1820) in South Africa. Doctoral Dissertation, University of Pretoria, South Africa. 2016.