The spatial distribution of African savannah herbivores: species associations and habitat occupancy in a landscape context

1. Introduction

Extinctions and population declines have disproportionately influenced large-bodied organisms [1]. For example, in Africa, 80% of the large-bodied taxa (which include members from the Primate, Perissodactyla and Carnivora orders) have been affected by harvesting [2]. Likewise, the consequences of losing species depend on their trophic position and functional redundancy. Thus, the loss of a top predator or dominant consumer often leads to a series of knock-on events that cause substantial changes in the distribution of biomass among trophic levels (e.g. [3,4]). These perturbations are especially important in ecosystems where predators or primary consumers lack functionally redundant trophic equivalents, i.e. systems in which dietary overlap is low among species with similar trophic positions [5–7].

Grass-dominated biomes throughout the tropics may be especially prone to species loss. Africa is home to more ungulates (hoofed mammals) than any other continent, with the greatest species richness (30 species) found in the grass-dominated savannahs of East Africa [8,9]. While the evolutionary timing and chronology of grassland-ungulate evolution is uncertain (e.g. [10–13]), molecular phylogenies, isotopic analysis and palaeo-reconstruction all confirm that a diverse East African ungulate community, with a diet dominated by C₄ savannah grasses, has existed for the last 8–10 Myr [11,14,15]. Additionally, data suggest that the carnivore community emerged more recently, likely within the last few million years [16]. We argue that future success in predicting the outcome of changes/losses to large mammal communities in tropical grass-dominated ecosystems is contingent upon a quantitative understanding of contemporary species interactions and associations in a fully functioning ecosystem.

Hierarchical spatial occupancy models provide probabilistic information on species occurrences in relation to key predictor variables (e.g. [17–19]). In conjunction with fine-scale time-series measurements, occupancy modelling can also address ecological questions involving population dynamics and species relationships, such as how species respond to resource availability and whether species associations (positive or negative) change in response to environmental variation. Camera trap studies have proved highly successful in characterizing ecological processes within ecosystems dominated by large-bodied, mobile mammals [20,21]. In this paper, we focus on the savannah habitats of East Africa because (i) this region harbours the greatest diversity of ungulate grazers on the Earth and (ii) numerous examples still exist of `naturally' functioning savannah ecosystems, devoid of major human disturbance.

Our study has two major aims: first, we introduce a simple conceptual model depicting the major interactions in savannahs among predators, herbivores, termites, fire, grasses and trees, which serves as a framework to understand herbivore habitat selection across African grass-dominated landscapes. Second, based on the framework from aim 1, we use a data-driven Bayesian model to quantify spatial habitat occupancy and species co-occurrence for a guild of eight herbivores and their most abundant predator (lions). The data originate from ‘Snapshot Serengeti’ [22], currently, the largest and longest running operational camera trap study.

2. Ecological interactions affecting herbivore habitat selection in African savannahs

African savannahs are typically represented by several guilds (figure 1): (i) a primary producer level consisting of sparse tree cover underlain by a highly flammable C₄ grassy understory, (ii) primary consumers consisting of (a) grazers and browsers that are subject to significant risks of predation by large carnivores and (b) mega-herbivores that are largely unaffected by predation risk [23] and (c) a diverse predator community, including large-bodied (typically lions) and meso-predators as dominant secondary consumers. This simplified schematic excludes small-bodied herbivores (dik-diks, rodents, lagomorphs, porcupines, etc.), scavengers (jackals, storks, vultures, etc.), insectivores (foxes, aardvarks, etc.), insects (e.g. grasshoppers, dung beetles, etc.) and ground-dwelling birds (bustards, cranes, secretary birds, etc.).

• Download figure
• Open in new tab
• Download PowerPoint

Fire distinguishes savannahs from other ecosystems by maintaining the tree–grass balance and ‘competing’ with herbivores for herbaceous vegetation [24]. Consequently, herbivores directly affect vegetation structure by consuming plants and indirectly by consuming herbaceous vegetation that fuels wildfire (figure 1). Owing to the regrowth of high-quality herbaceous forage following fire (e.g. high protein : lignin content or low nutrient : carbon ratio), herbivores often graze preferentially in recently burned areas [25]. Mega-herbivores, but elephants in particular, can facilitate forage availability for smaller bodied browsers by pushing over trees and encouraging re-sprouting, thus increasing availability for smaller-bodied herbviores [26]. On the other hand, competition for browse may have negative impacts, and it is unclear when positive or negative interactions will dominate (figure 1). Predators can have significant top-down influences both through direct consumption of herbivore prey and also through non-consumptive effects of the ‘landscape of fear’ that influence when and where herbivores concentrate their foraging [27,28]. Predator–prey relationships are body size-dependent, with a strong non-linear (sigmoidal) relationship between mass and per capita predation effects [29]. Above 150 kg, herbivores are relatively insensitive to predation and are regulated instead by food abundance, mid-sized herbivores are regulated by food quality, and smaller herbivores are largely regulated by predators [30]. As a consequence, herbivores are expected to respond to perceived predation risk in a way that depends strongly on body size, with small-bodied species avoiding habitats in which probability of predation is high. In African savannahs, areas with good cover for sit-and-wait predators, such as river confluences, woody vegetation and erosion terraces, have been identified as areas where prey ‘catchability’ increases [31]. Herbivores of large body size or those that move in large, aggregated groups can use risky habitats because, per capita, their predation risk is reduced, thus stabilizing predator–prey dynamics [32].

Ungulates also respond to landscape heterogeneity according to the spatial and temporal distribution of resources. A major axis of variation among mammalian herbivores is created by the foraging constraints which arise due to herbivore body size: the fermentation required to digest a low-quality, cellulose-rich diet requires long gut passage times and, therefore, a relatively large body cavity [33–35]. Thus, while large-bodied ruminant grazers, like wildebeest and buffalo, can tolerate low-quality forage because of their large gut capacity, smaller-bodied ruminants, such as gazelles, are constrained to eating relatively high-quality forage (plants with lower cellulose content and higher nutrient to carbon ratios). For example, islands of soil fertility and nutrient-rich forage, such as those created by termite activity [36], may be used differently based on dietary requirements and digestive physiology, with smaller-bodied ruminants more closely tied to nutrient hotspots than larger, hind gut fermenters requiring large volumes of lower quality forage (e.g. [33]). However, although termite mounds attract foraging grazers, browsers and mega-herbivores in some systems (e.g. [37–39]), they are actively avoided by browsers in others (e.g. [40,41]). The emerging picture is that vegetation composition and the quality of the forage in the matrix surrounding termite mounds determines foraging selectivity by ungulates [42]. These interactions provide a generalized framework for selecting environmental predictors, generating hypotheses and ultimately developing an integrated approach to study herbivore occupancy in a heterogeneous African savannah.

3. Modelling habitat occupancy of Serengeti herbivores

The topology of food web connections can stabilize communities with weak interactions and multi-linkage chains that are either nested [43] or asynchronous (e.g. [44]). Consumer–resource interactions (i.e. plant–herbivore and predator–prey) in the Serengeti ecosystem in northern Tanzania, have been described as (i) composed of fast and slow energy channels in grassland and woodland habitats, respectively [45] and (ii) having hierarchically nested consumer–resource interactions, with larger body-size consumers having generalist foraging patterns and smaller-bodied consumers specializing on a nested subset of prey [30,45]. Consequently, populations of savannah herbivores in Serengeti are expected to display stability under perturbation (i.e. absence of dynamic fluctuations), a prediction supported by long-term predator–prey population data under El Nino Southern Oscillation climatic variation [46]. Moreover, both theory (e.g. [47,48]) and analysis of Serengeti food web data [45,49] suggest that large mobile species link habitats through their movements and are key to understand large naturally functioning food webs. For example, the annual migration of approximately 1.3 million wildebeest and zebra physically connects otherwise disparate habitat patches, thus coupling productivity, energy and nutrients between the open Serengeti plains and the woodlands of Kenya's Masai Mara Reserve [45,49].

Despite recent theoretical development in stochastic food web theory, an empirical analysis of covariance dynamics (e.g. species associations in space and time) is still lacking between different species and subsystems across spatial scales [48], owing to an absence of spatially explicit time-series data for an entire guild of species. In the following section, we use data from the Snapshot Serengeti camera trap study to quantify associations among herbivores and between herbivores and their top predator across a large heterogeneous savannah landscape (i.e. aim 2). The multivariate model, which we develop in the following sections, provides quantitative estimates of herbivore associations with other aspects of the landscape, such as variation in forage quality, termite mound density, river distance, disturbance from fire and tree cover, as a means to explore variation in herbivore habitat preferences. We initially predicted that habitat preferences would depend strongly on herbivore body size, with larger species using habitats with abundant forage (e.g. high NDVI; e.g. [50]) and small-bodied herbivores avoiding areas with significant predation risk (woodlands and riverine areas; e.g. [31]) and occupying recently burned areas due to high forage quality (e.g. [25]). We also predicted that migrants, such as wildebeest, zebra and Thomson's gazelles, would show strong positive spatial associations, as would residents, such as topi and hartebeest. Finally, we predicted that migratory species would show indiscriminate use of the landscape, whereas resident species would be constrained to isolated, high-quality habitat patches.

4. Material and methods

(a) Study site

We deployed 225 cameras in a 1125 km² grid in the centre of Serengeti National Park (SNP), Tanzania (figure 2), a savannah ecosystem at the Kenya–Tanzania border in eastern Africa dominated by the annual migration of 1.6 million wildebeest and zebra. The migration follows a seasonal rainfall gradient from wetter northwest woodlands to drier, nutrient-rich southeast short-grass plains. The grid itself covers the intersection of open plains and woodland savannah and lies within an area of long-term (approx. 40 years) lion monitoring (Serengeti Lion Project, [22]). The camera trap survey was piloted in 2010 and has operated continuously since 2011. At each camera trap location, we collected metadata data on tree density, grass height and shade availability, and updated these values every 1–2 years.

• Download figure
• Open in new tab
• Download PowerPoint

(b) Cameras and layout

We used Scoutguard (SG565) incandescent cameras with passive infrared sensors that are triggered by a combination of heat and motion. We placed the cameras at the centre of each 5 km² grid cell on the nearest available tree or, lacking trees, on metal poles. Cameras were set approximately 50 cm above the ground to capture medium- to large-sized vertebrates and checked every two months. Low-hanging branches and tall grass were trimmed at regular intervals to minimize camera misfires and provide an unobstructed view.

(c) Data collection and processing

During the period from 2010 to 2013, the survey produced approximately 1.2 million image sets [22]. In partnership with Zooniverse, more than 70 000 members of the general public classified photographic images on our citizen science website, Snapshot Serengeti (www.snapshotserengeti.org). Users were shown photographs at random and were asked to use a series of web-based visual tools to classify all mammals larger than a mongoose [51]. The classifications of multiple users per image are aggregated into a ‘consensus dataset’ of final classifications. Classifications include the identification of all species in an image, the number of individuals of each species and the presence/absence of juveniles and specific behaviours. The consensus classifications have been validated against expert-identified image sets to be 96.6% accurate in terms of species identifications and 90% accurate for species counts. An exhaustive review of identification and misidentification rates by species can be found in [51]. Our analysis focuses on eight of the most common and readily identifiable herbivores in the study area, including elephant (Loxodonta africana), buffalo (Syncerus caffer), zebra (Equus quagga), wildebeest (Connochaetes taurinus), hartebeest (Alcelaphus buselaphus), topi (Damaliscus korrigum), impala (Aepyceros melampus) and Thomson's gazelle (Eudorcas thomsonii), and their most abundant predator, the African lion (Panthera leo). Body size for each species was estimated as the midpoint of the adult female body mass range as reported in Kingdon [52].

(d) Predictor variables

As predictors of ungulate occupancy, five variables were selected to represent forage quality/quantity, vegetation structure, disturbance (e.g. fire), soil fertility and probability of predator encounter. The Normalized Difference Vegetation Index (NDVI), a measure of vegetation greenness, was obtained from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) instrument operating on the Terra satellite platform. We used the MOD13Q1 product for the year 2012, which provides NDVI (and other vegetation indices) as 16-day composite images at a spatial resolution of 250 m for the entire grid area. Based on the previous work [50,53,54] and the observation that grazers respond to NDVI on relatively short timescales, our a priori assumption was that NDVI largely served as a proxy for changes in the herbaceous layer. As an indicator of woody vegetation structure (e.g. relative density of trees), we enumerated trees greater than 2 m in height in a 50 m radius (approx. 7854 m²) around each camera location and classified them into one of four categories: none (0), sparse (1–10), moderate (11–30) or dense (greater than 30) tree cover. The inclusion of both NDVI and a tree-density classification accounted for effects of both herbaceous and woody vegetation. Distance to nearest river was incorporated as an indicator of habitat selection based on species-specific requirements for water and on risk associated with increased predation near rivers and confluences [31]. Camera distance from the nearest permanent river was extracted from a GIS layer of river distances for the entire Serengeti available from Hopcraft et al. [31]; river distance was square-root transformed prior to inclusion in the model. The location of fires across the grid was determined using the MODIS burned area product (MCD45A1), which provides a burn date for pixels at a 500 m spatial resolution [55]. Fire = 1, if it occurred in the two previous sampling periods (i.e. 32 days), otherwise fire = 0. A seasonality (e.g. wet versus dry season) variable was included by determining rainfall with a Decagon Devices (Pullman, WA, USA) ERCN-100 tipping spoon rain gauge (1 tip = 0.2 mm) that were part of two weather stations within the study area, one in the northwest and the other in the southeast region of the grid. Cameras were assigned rainfall values of the nearest weather station (mean distance = 11.2 km, range = 0.64–23.5 km), and the wet season was defined as months when rainfall exceeded 100 mm. Termite mound density was enumerated within a 50 m radius of each trap as an indicator of soil nutrient concentration (e.g. [36]) and their potential as herbivory hotspots in the landscape (e.g. [56]). All mounds of approximately 1 m in height or greater were counted and an indicator variable was created to account for termite mounds in the model (0 = absent, 1 = present). These mounds are built by Odontotermes, Macrotermes or Trinervitermes genera [57] and contain locally enriched soils, stoloniforous grasses, forbs and shrubs (e.g. [38]). For temporally variable data (NDVI, fire and rainfall), we extracted values for each 16-day interval over which the animal's presence/absence was determined (22 time steps over the course of a year). In the case of static landscape variables (e.g. termite mound density), we assigned a single value at each trap and treated it as a site variable in which temporal variables were nested.

(e) Bayesian modelling

To quantify the spatial variability of ungulate occupancy over time, we modelled the simultaneous presence/absence of any two species using a bivariate Bayesian occupancy model that accounts for both within- and between-species spatial dependence as well as imperfect detection. To begin, we let the subscript s index each species, i index the camera trap location and t index time. The response variable Y_sit takes value 1 if species s is present at camera location i during time t, and otherwise takes value 0; thus, we model Y as a Bernoulli random variable with covariates. As this is a multivariate model, the dataset has a unique row for each combination of i and t, and all species show either a 1 or 0 value for the same record. Thus, one record contains all the presence/absence information for each species and all covariate data. For the bivariate species model (i.e. associations between any two species), we extracted the information for the two focal species plus all covariate data.

There is a well-documented issue of detectability with camera trap studies (e.g. [18,20]), wherein it is important to acknowledge that an observed value Y_sit = 0 may mean either that a species was not present, or that the species was in fact present but could not be detected. We call Z_sit the latent (unmeasured) variable for the true presence versus absence of species s at location i and time t. We assume that detection (Y_sit = 1) or non-detection (Y_sit = 0) depends on the latent presence (Z_sit = 1) or absence (Z_sit = 0) and detection probability p_sit. The observed presence/absence of each species is modelled as

Only when Z_sit = 1 does the model allow for the possibility that Y_sit = 1, and further, this detection occurs with probability p_sit. Therefore, p_sit is the probability of detecting species s at location i and time t. Moreover, the model assumes independence of Y_sit conditional on Z_sit and p_sit. We model the logit of this detection probability as where are observed landscape features and β_s is a vector of species-specific parameters for detection. We use latitude, longitude and an indicator for woodlands (greater than 10 trees within 50 m of the camera) for these landscape features.

Next, our model accounts for both within- and between-species spatial dependence using an extension of the univariate occupancy model proposed by Royle & Dorazio ([18], ch. 9). It is common to use an autologistic model for the presence/absence within a Bayesian occupancy model (e.g. [21]). An autologistic model considers occurrence at neighbouring locations when modelling the probability of occurrence for any particular location. In this paper, we consider both same- and different-species occurrence at neighbouring locations. Traditional autologistic models bias estimates of model parameters [58], but can be corrected by a centred parametrization of the autologistic model [59]. The presence/absence component of our model is an adaptation of the bivariate, centred autologistic model of Caragea & Berg [60]. We model the presence Z_sit conditional on all other species (anywhere and at any time) and also the same species s at all other locations and times. It is convenient to notate this set of the presence/absences as Z_-sit, which is all Z_s′i′t′ excluding Z_sit. Let where is a vector of landscape features and α_s is a vector of species-specific parameters which influence occupancy. Under the assumption of no spatial autocorrelation, μ_sit is interpreted as the probability of occurrence based solely on the landscape features. This term is needed in the centred autologistic model described above. The full hierarchical model is given by:

and where φ_sit is the probability of occurrence of species s at location i and time t, modelled as

In the above model, δ_s is a measure of the strength of the within-species spatial autocorrelation for species s, and δ is the strength of the cross-species spatial autocorrelation between species s and s′. The neighbourhood of location i at time t, N_t(i), consists of all camera trips within 2.2 km. The covariates used to model species presence are: NDVI, rainy season, square root of the distance to the nearest river (the square-root transformation was used purely for computational reasons), the presence of any termite mounds within 50 m, any incidence of fire within the previous 32 days and three indicators for tree cover (1–10, 11–30 and greater than 31 trees within 50 m). The covariates used to model the probability of detection are longitude, latitude and an indicator for greater than 10 trees within 50 m. The parameters for detection and presence/absence probability (β₁, β₂, α₁, α₂, δ₁, δ₂, δ) are assigned zero-mean normal prior distributions, with the autoregressive coefficients δ₁, δ₂ and δ each truncated to [−1, 1]. The posterior distribution is simulated using an adaptive Metropolis-within-Gibbs algorithm that used the pseudo-likelihood to update the occupancy parameters due to the intractability of the full conditional distributions [61]. The algorithm was coded in MATLAB and implemented on a single core of a research computing cluster. We report the mean coefficients ± standard deviations for each model parameter; however, those parameters with s.d. ≥ 50% of the mean are considered weak signals and are largely ignored.

Finally, we created maps of habitat occupancy probability for herbivores to explore how fine-scale (i.e. 250 m) occupancy varied as a function of body size. We restricted the mapping to the six grazing herbivores (buffalo, zebra, wildebeest, hartebeest, topi and Thomson's gazelle), because they span a large range of body sizes and their reliance on herbaceous forage means that they are more directly comparable for purposes of understanding resource selection. For each map, occupancy was scaled relative to the maximum probability of observing a species, so that cells vary in their occupancy from 1 (most probable) to 0 (least probable).

5. Results

(a) Spatial environmental variation

Environmental variables showed strong spatial gradients that varied with latitude and longitude across the camera grid (table 1; electronic supplementary material, figure S1). NDVI, the presence of termites and fire showed a strong spatial trend which increased from east to west, consistent with the greater rainfall, deeper soil and higher herbaceous biomass (and thus fuel) in that part of the study area. As latitude increased, dense tree cover, proximity to rivers (1/river distance), termite presence and fire increased. Other notable correlations included a positive association between sites that burned and dense tree cover and a negative association between moderate tree cover and river distance (table 1; electronic supplementary material, figure S1).

(b) Herbivore associations with landscape features

Species diverged in their associations with landscape features, often in a way that depended on herbivore body size (table 2 and figure 3). For example, large- and small-size herbivores had opposite relationships with NDVI: elephant, buffalo and wildebeest were positively associated while hartebeest, topi and Thomson's gazelle were negatively associated with NDVI (figure 3a). Large-bodied grazers, such as zebra and wildebeest avoided burned areas, while Thomson's gazelle and impala were associated with burned areas (figure 3b). Topi, in particular, showed strong positive associations with termite mounds, while impala and elephants showed weaker positive associations; hartebeest, a resident species similar in body size to topi, were negatively associated with termite mounds (figure 3c). Herbivore occupancy in relation to rivers and tree density was also associated with body size in that the largest herbivores, elephants and buffalo, were negatively associated with river distance and occupied areas of moderate and dense tree cover. The exception to this was that impala, a relatively small browser, preferred habitat of moderate and dense tree cover associated with rivers (table 2; figure 3d). Finally, woodland occupancy (in this case the measure of sparse tree cover) showed a striking linear relationship to herbivore log body size (figure 3e). To investigate this pattern further, we plotted the coefficient from the Bayesian occupancy model for sparse tree cover as a function of herbivore log body mass (kilograms). The relationship was positive and strongly linear across all eight species (figure 4): woodland coefficient = −2.79 + 0.58 × ln(mass (kg)).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

(c) Inter-specific spatial autocorrelation in occupancy

The tendency to observe the same species at neighbouring traps (approx. 2.2 km²) during the same sampling window, i.e. ‘intra-specific spatial dependency’, varied according to body size and migratory strategy: the largest (elephant and buffalo) and smallest species (topi–Thomson's) showed no pattern, while large-bodied migrants (zebra and wildebeest) and the largest resident species (hartebeest) showed positive spatial dependence (figure 3f).

(d) Pairwise associations between species

After accounting for the covariate information, our Bayesian analysis demonstrated surprisingly strong, positive pairwise species associations among herbivores across 16-day intervals (table 3 and figure 5). Surprisingly, impala showed positive spatial associations with all other herbivores in the dataset (n = 7), the strongest relations being with topi and hartebeest. The strong association between impala and other species may reflect their flexibility in foraging as both a browser and grazer. Remaining herbivores had either five or six positive spatial associations with other herbivores, except for elephants which, after impala, showed associations only with wildebeest and buffalo. Perhaps, not surprisingly, the strongest association involved the two most abundant herbivores in the system, wildebeest and zebra (0.86). Other strong spatial associations emerged between hartebeest and topi (0.51), impala and Thomson's (0.50), wildebeest and buffalo (0.49) and impala and topi (0.47). Lions had the strongest positive association with wildebeest, buffalo and zebra, putatively because lions track their large-bodied mobile prey, especially during the wet season [56]. Topi, a resident herbivore, was the only species that showed a negative association with lions (figure 4). Although this effect is relatively weak (table 3), we chose to include it, because it demonstrates an important contrast to the association between lions and their abundant, mobile prey. Important future work will investigate the degree to which this changes with the length of the temporal window over which species are aggregated (see Discussion).

• Download figure
• Open in new tab
• Download PowerPoint

Finally, the scaled probability of occupancy maps demonstrated considerable variation among the six grazing species at the extent of the entire landscape. The first striking pattern to emerge is mobility/occupancy in relation to body size (figure 6; electronic supplementary material, figure S2). Small and resident herbivores (Thomson's gazelle, topi and hartebeest) maintained local patches of preferred habitat, with most of the landscape showing low probability of occurrence (less than 0.30; figure 6a–c). By contrast, for the largest grazer, buffalo, much of the landscape is represented by a high-scaled probability of occurrence (greater than 0.60; figure 6f), whereas low probability of occurrence is only seen in areas to the southeast where standing biomass of herbaceous vegetation is low. Wildebeest and Zebra show intermediate usage of the landscape, with most of the landscape dominated by intermediate probability of occurrence (less than 0.70 and greater than 0.40). Thomson's gazelle occupancy decreased from southeast to northwest as grassland productivity and NDVI increased. However, their occupancy was concentrated in the portions of the landscape with the lowest NDVI, likely reflecting patches that were favourable because of low-growing grasses that offered both high-quality forage and safety from predators. Like Thomson's gazelle, topi are relatively patchy in their use of the landscape, occupying patches associated with termite mounds and away from river courses in the northwest (figure 6; electronic supplementary material, figure S1). Hartebeest were likewise patchy in their use of the landscape but a high probability swath was observed in the relatively open habitat in the eastern areas of the grid.

• Download figure
• Open in new tab
• Download PowerPoint

6. Discussion

Most camera trap studies model species occupancy in relation to habitat features, resources or enemies, without considering interactions among species on the same trophic level. Here, we model occupancy for a guild of African mammals over 16-day intervals at a relatively fine spatial resolution (2.2 km²) and over a relatively large landscape (1100 km²), while explicitly accounting for spatial associations among herbivores. The resultant analysis permits the evaluation of questions pertaining to food web behaviour, consumer–resource dynamics and the trade-off between foraging and risk.

(a) Species associations in a landscape context

The strong spatial associations which emerged from the analysis after accounting for the covariates was surprising, with likely explanations falling into three categories: facilitation of forage quality for smaller herbivores by larger herbivores, ideal free distributions based on an underlying resource gradient and diffusion of predation risk [62]. In the first category, large-bodied herbivores such as zebras and wildebeest may facilitate consumption by smaller bodied herbivores such as Thomson's gazelle, topi and hartebeest [45,63]. Others have rejected this hypothesis using population [64] and habitat overlap data [65]. In our study, the strong spatial association between species that are unlikely to elicit facilitation, such as between impala and topi, for example, also suggest that an explanation rooted in purely facilitation cannot account for our results. In the second category, ungulates may distribute themselves according to an ideal free distribution. Under this scenario, some ungulate species may simply occupy the most favourable patches in the landscape and distributions are driven by local resource availability rather than species interactions among herbivores. However, in the light of evidence from Serengeti that Thomson's gazelle populations are negatively influenced by inter-specific competition [64], we find this explanation unlikely. Third, group foraging may reduce per capita predation risk and reduce predator efficiency, a mechanism that has been proposed as fundamental to predator–prey dynamics in many ecosystems (e.g. [32]). Moreover, associations among conspecific individuals have positive and stabilizing effects on prey populations in Serengeti [32] and our data suggest that these associations may extend to the formation of mixed-species herds. The observation that elephants, not at significant risk of predation by lions [29], had relatively weak and few associations with other herbivores is consistent with the anti-predator hypothesis. However, as pointed out by others [62], these mechanisms need not be mutually exclusive and can accentuate herbivore aggregations.

(b) Landscape indicators of habitat Use

Several key landscape features emerged as indicators of habitat selection for herbivores, including NDVI, termite mound density, recent fire, distance to rivers and woody vegetation (figure 3). Our results highlight the well-known importance of body size as a key trait that influences herbivore habitat affinity (e.g. [30]). While the observational nature of our study did not physically separate the effects of forage resources from predation risk, the inclusion of both herbaceous biomass (NDVI) and woodland cover in our model provides compelling insight into patterns of herbivore occupancy. For example, the strong positive response of buffalo and negative response of Thomson's gazelle to NDVI is consistent with the trade-off experienced by herbivores of different body sizes: large herbivores attempt to maximize forage intake rates, while small herbivores attempt to maximize forage quality [30,33,34]. We acknowledge potential difficulties linking satellite-derived data to the forage quantity/quality trade-off, but note that our results are consistent with ground-truthed vegetation data that demonstrate the greater quality (low C : N ratios) of low biomass vegetation in regions of the Serengeti that map as low NDVI (e.g. [50,54]).

Smaller herbivores also associated with other habitat features, such as recent burns and termite mounds, which improve forage quality relative to the surrounding matrix. Our results confirm the strong allometric relationship between body size and the use of burned patches documented in savannah ecosystems [66,67]. Previous work identified strong links between nutritive quality of re-growing forage in recent burns and habitat preference by herbivores lasting up to three months [25]. This, together with evidence that vigilance by Thomson's gazelles does not differ in burned and unburned patches, suggests that herbivores are primarily using burn patches for their nutritional properties and not their decreased risk of predation (e.g. [68]).

Topi use termite mounds as a vantage point to signal territoriality (e.g. [69]) or carry out mating displays [70] and our data clearly support a strong association between topi and termite mounds [71]. However, our model also supports the use of areas of high termite mound density by impala and elephants and the avoidance of these same areas by hartebeest (figure 3). In Lake Mburo National Park, Uganda, both grazing and browsing herbivores used termitaria as preferential feed sites over sites without termite mounds [38]. In Kruger National Park, South Africa, the strong preference of herbivores for vegetation on and around termite mounds extends up to 20 m beyond the mound [72]. The selection of areas with greater termite mound density by elephant may be due to their preference for woody vegetation growing on mounds, which is enriched in nutrient content relative to the surrounding vegetation [37].

Smaller herbivores, such as Thomson's gazelle, topi and hartebeest, also avoided habitats used by lions, such as rivers and woodland cover (e.g. [31]), in ways suggesting that they were using habitat in a way that reduced risk of predation. The one exception was impala (e.g. table 2), which is a relatively small-bodied mixed feeder and known as an obligate woodland species. However, impala occur in large groups consisting of one male and up to 30 females, which may reduce the predation risk by providing greater vigilance by the group [73], and, indeed, lion predation on impala is extremely rare in the Serengeti: lions were only witnessed to kill three impala out of 1553 observed prey captures between 1966 and 2014 (C. Packer 2016, personal observation).

When viewed on the extent of the entire landscape, large-bodied species and migrants tend to use and move across the entire landscape while avoiding areas of low-standing biomass, such as areas of low NDVI, high termite mound density or recent fire. By contrast, smaller herbivores, such as the semi-migrant Thomson's gazelle and resident hartebeest, use a series of smaller, high-quality patches that are likely to offer high forage quality and improved visibility of plains predators such as cheetahs and wild dogs ([50,54]; figure 6). In fact, when the landscape occurrence probabilities of grazers are viewed as frequency distributions (electronic supplementary material, figure S2), three distinct groups appear: mobile mega-herbivores (buffalo); large-bodied migrants (zebra and wildebeest) and small-bodied (Thomson's gazelle) or resident species (topi and hartebeest). This suggests that, at the scale of this landscape (approx. 1100 km²), body size is also linked to movement patterns within the study area with larger species covering greater fractions of the landscape relative to smaller species.

7. Conclusion and further directions

As described in our conceptual framework savannahs are composed of complex, multivariate relationships among predators, herbivores and vegetation and disturbance (figure 1). Herbivore body size emerges as a major axis to understand within- and among-trophic level interactions in terms of landscape occupancy and spatial associations among pairs of mammalian herbivores. While our analysis identified a high degree of spatial association among herbivores across 16-day periods, the logical next step will be to vary the temporal window size and quantify the effects on species associations. If species that were positively associated at 16 days switch to negative associations at smaller temporal scales, then evidence exists for facilitation between large- and small-bodied herbivores (e.g. [62,63]). On the other hand, if species maintain their positive associations at smaller temporal windows, then evidence exists for inter-specific spatial associations to diffuse the predation risk (e.g. [64]).

Our study focused on data collected across one wet/dry season cycle, so it is unknown whether the habitat affinities (e.g. figure 3) and species linkages (e.g. figure 4) are stable through time. If perturbations, such as large changes in annual primary production due to the El Nino Southern Oscillation [46] were to alter species associations, putative effects from grazing facilitation could be more easily disentangled from anti-predator strategies. Years with greater than average primary productivity, and thus ample cover for ambush predators, are linked to increased rates of prey capture rates by lions and high recruitment rates of new prides [74]. Research on trophic cascades in savannahs (e.g. [75]) is beginning to explore how predator–prey relationships can influence other components of the systems, but the degree to which environmental perturbations drive predator-mediated behavioural shifts in herbivores is unknown. If predators show temporal variation in preference for residents and migrants (prey switching), we expect to see variation in the interaction network diagram (figure 4) over the years.

What has been more convincing is evidence that the elimination of large-bodied herbivores leads to complex and sometimes unpredictable trophic cascades in savannahs [4,76]. Notably, members of the large mammalian herbivore community appear to lack functional redundancy [77], suggesting that losing individual species will have unique and largely unpredictable effects on communities, capable of cascading from large mammals, to rodents, snakes and ultimately parasites that harbour wildlife disease (e.g. [78]). Consequently, more research is needed to understand how the loss of individual species may propagate through the Serengeti food web.

We conclude by acknowledging the correlational nature of the camera trap data and measures of species associations. Our species associations results (figure 4), while related to food web dynamics, are not quantitative measures of trophic interactions, but rather spatial associations which could arise via several different processes. Likewise, the strong relationships between herbivores and the environmental variables in our models, such as NDVI and fire, are suggestive, but they do not identify an actual underlying mechanism. Thus, only manipulative experiments involving fire, resources and predation risk can test explicit mechanisms identified in this study. A blend of large-scale observational studies with experimental approaches will elucidate how perturbations may propagate across landscapes and influence the future stability of the east African tropical grassy biome.

Data accessibility

Data, metadata and images can be accessed at: http://dx.doi.org/10.5061/dryad.5pt92.

Authors' Contributions

T.M.A. conceived of the study, conducted the data analysis and drafted the manuscript. S.W., B.D. and R.E. participated in data analysis and carried out the statistical analyses. M.P. and A.S. collected field data and helped draft the manuscript. M.K. managed and analyzed data. C.P. designed the study and helped draft the manuscript. All authors gave final approval for the publication.

Competing interests

We have no competing interests.

Funding

Funding for this project was provided by NSF DEB-1020479 to C.P., NSF DEB-1145861 to T.M.A., NSF BCS-1461728 to T.M.A. and separate explorer grants from National Geographic to C.P. and T.M.A. Snapshot Serengeti website development was funded by awards to Zooniverse from the Alfred P. Sloan Foundation.

Footnotes

One contribution of 15 to a theme issue ‘Tropical grassy biomes: linking ecology, human use and conservation’.